plot-till-evolution.py (3787B)
1 #!/usr/bin/env python3 2 3 from pylab import figure, subplots, plot, xlabel, ylabel, title, axis, vlines, savefig, text, tight_layout, cm, legend 4 from sys import exit 5 6 import MISMIP 7 8 import numpy as np 9 from optparse import OptionParser 10 import os.path 11 12 try: 13 from netCDF4 import Dataset as NC 14 except: 15 print("netCDF4 is not installed!") 16 sys.exit(1) 17 18 def process_options(): 19 "Process command-line options and arguments." 20 parser = OptionParser() 21 parser.usage = "%prog <input files> [options]" 22 parser.description = "Plots the ice flux as a function of the distance from the divide." 23 parser.add_option("-o", "--output", dest="output", type="string", 24 help="Output image file name (e.g. -o foo.png)") 25 26 opts, args = parser.parse_args() 27 28 if len(args) == 0: 29 print("ERROR: An input file is requied.") 30 exit(0) 31 32 if len(args) > 1 and opts.output: 33 print("More than one input file given. Ignoring the -o option...\n") 34 opts.output = None 35 36 return args, opts.output, opts 37 38 39 def read(filename, name): 40 "Read a variable and extract the middle row." 41 nc = NC(filename) 42 43 try: 44 var = nc.variables[name][:] 45 except: 46 print("ERROR: Variable '%s' not present in '%s'" % (name, filename)) 47 exit(1) 48 49 return var 50 51 52 def find_grounding_line(x, topg, thk, mask): 53 "Find the modeled grounding line position." 54 # "positive" parts of x, topg, thk, mask 55 topg = topg[x > 0] 56 thk = thk[x > 0] 57 mask = mask[x > 0] 58 x = x[x > 0] # this should go last 59 60 def f(j): 61 "See equation (7) in Pattyn et al, 'Role of transition zones in marine ice sheet dynamics', 2005." 62 z_sl = 0 63 return (z_sl - topg[j]) * MISMIP.rho_w() / (MISMIP.rho_i() * thk[j]) 64 65 for j in range(x.size): 66 if mask[j] == 2 and mask[j + 1] == 3: # j is grounded, j+1 floating 67 nabla_f = (f(j + 1) - f(j)) / (x[j + 1] - x[j]) 68 69 # See equation (8) in Pattyn et al 70 return (1.0 - f(j) + nabla_f * x[j]) / nabla_f 71 72 raise Exception("Can't find the grounding line") 73 74 75 def plot_profile(in_file, out_file): 76 77 if out_file is None: 78 out_file = os.path.splitext(in_file)[0] + "-till-evol.pdf" 79 80 steps = read(in_file, 'thk').shape[0] 81 82 fig, ax = subplots(1, 1, sharex=True, figsize=[6, 3]) 83 for i in range(1, steps): 84 85 mask = read(in_file, 'mask')[i] 86 uvelbase = read(in_file, 'uvelbase')[i] 87 utillflux = read(in_file, 'utillflux')[i] 88 till_deposit = read(in_file, 'tilldeposit')[i] 89 x = read(in_file, 'x') 90 91 # convert x to kilometers 92 x /= 1e3 93 94 # modeled grounding line position 95 #xg_PISM = find_grounding_line(x, lsrf, thk, mask) 96 #plot(x, np.zeros_like(x), ls='dotted', color='red') 97 icecolor = cm.cool(i / steps) 98 #ax[0].plot(x, uvelbase, color=icecolor) 99 #ax[1].plot(x, utillflux, color=icecolor) #, label='{}'.format(i)) 100 ax.plot(x, till_deposit, color=icecolor) 101 102 ax.set_xlabel('distance from the divide, km') 103 #ax[0].set_ylabel('$v_{SSA}$, m/a') 104 #ax[1].set_ylabel('$q_{t,x}$, m$^2$/a') 105 ax.set_ylabel('$\Delta b$, m') 106 107 #_, _, ymin, ymax = axis(xmin=0, xmax=x.max()) 108 _, _, ymin, ymax = axis(xmin=950, xmax=1150) 109 #_, _, ymin, ymax = axis(xmin=x.min(), xmax=x.max()) 110 111 #vlines(xg / 1e3, ymin, ymax, linestyles='dashed', color='black') 112 #vlines(xg_PISM / 1e3, ymin, ymax, linestyles='dashed', color='red') 113 114 #legend() 115 fig.tight_layout() 116 print("Saving '%s'...\n" % out_file) 117 savefig(out_file) 118 119 if __name__ == "__main__": 120 args, out_file, opts = process_options() 121 122 for in_file in args: 123 plot_profile(in_file, out_file)