plot-time-series.py (4012B)
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] + "-profile.pdf" 79 80 steps = read(in_file, 'thk').shape[0] 81 82 fig, ax = subplots(1, 1, figsize=[4, 3]) 83 for i in range(1, steps): 84 85 mask = read(in_file, 'mask')[i] 86 usurf = read(in_file, 'usurf')[i] 87 topg = read(in_file, 'topg')[i] 88 lsrf = topg.copy() 89 thk = read(in_file, 'thk')[i] 90 sea_lvl = read(in_file, 'sea_level')[i] 91 #till_deposit = read(in_file, 'tilldeposit')[i] 92 x = read(in_file, 'x') 93 94 # mask out ice-free areas 95 usurf = np.ma.array(usurf, mask=mask == 4) 96 97 # convert x to kilometers 98 x /= 1e3 99 100 # compute the lower surface elevation 101 lsrf[mask == 3] = usurf[mask == 3] - thk[mask == 3] 102 lsrf = np.ma.array(lsrf, mask=mask == 4) 103 # modeled grounding line position 104 #xg_PISM = find_grounding_line(x, lsrf, thk, mask) 105 #plot(x, np.zeros_like(x), ls='dotted', color='red') 106 icecolor = cm.cool(i / steps) 107 ax.plot(x, topg, color='black') 108 ax.plot(x, usurf, color=icecolor, label='{}'.format(i)) 109 ax.plot(x, lsrf, color=icecolor) 110 111 ax.set_xlabel('distance from the divide, km') 112 ax.set_ylabel('elevation, m') 113 114 #_, _, ymin, ymax = axis(xmin=0, xmax=x.max()) 115 _, _, ymin, ymax = axis(xmin=950, xmax=1150, ymin=-500, ymax=1150) 116 #_, _, ymin, ymax = axis(xmin=x.min(), xmax=x.max()) 117 118 #vlines(xg / 1e3, ymin, ymax, linestyles='dashed', color='black') 119 #vlines(xg_PISM / 1e3, ymin, ymax, linestyles='dashed', color='red') 120 121 #legend() 122 fig.tight_layout() 123 print("Saving '%s'...\n" % out_file) 124 savefig(out_file) 125 126 if __name__ == "__main__": 127 args, out_file, opts = process_options() 128 129 for in_file in args: 130 plot_profile(in_file, out_file)