tplot-time-series.py - pism-exp-gsw - ice stream and sediment transport experiments
(HTM) git clone git://src.adamsgaard.dk/pism-exp-gsw
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tplot-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