tplot.py - pism-exp-gsw - ice stream and sediment transport experiments | |
git clone git://src.adamsgaard.dk/pism-exp-gsw | |
Log | |
Files | |
Refs | |
README | |
LICENSE | |
--- | |
tplot.py (7433B) | |
--- | |
1 #!/usr/bin/env python3 | |
2 | |
3 import MISMIP | |
4 | |
5 from pylab import figure, subplot, plot, xlabel, ylabel, title, axis, vl… | |
6 from sys import exit | |
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 | |
19 def parse_filename(filename, opts): | |
20 "Get MISMIP info from a file name." | |
21 tokens = filename.split('_') | |
22 if tokens[0] == "ex": | |
23 tokens = tokens[1:] | |
24 | |
25 try: | |
26 model = tokens[0] | |
27 experiment = tokens[1] | |
28 mode = int(tokens[2][1]) | |
29 step = int(tokens[3][1]) | |
30 | |
31 except: | |
32 if opts.experiment is None: | |
33 print("Please specify the experiment name (e.g. '-e 1a').") | |
34 exit(0) | |
35 else: | |
36 experiment = opts.experiment | |
37 | |
38 if opts.step is None: | |
39 print("Please specify the step (e.g. '-s 1').") | |
40 exit(0) | |
41 else: | |
42 step = opts.step | |
43 | |
44 if opts.model is None: | |
45 print("Please specify the model name (e.g. '-m ABC1').") | |
46 exit(0) | |
47 else: | |
48 model = opts.model | |
49 | |
50 try: | |
51 nc = NC(filename) | |
52 x = nc.variables['x'][:] | |
53 N = (x.size - 1) / 2 | |
54 if N == 150: | |
55 mode = 1 | |
56 elif N == 1500: | |
57 mode = 2 | |
58 else: | |
59 mode = 3 | |
60 except: | |
61 mode = 3 | |
62 | |
63 return model, experiment, mode, step | |
64 | |
65 | |
66 def process_options(): | |
67 "Process command-line options and arguments." | |
68 parser = OptionParser() | |
69 parser.usage = "%prog <input files> [options]" | |
70 parser.description = "Plots the ice flux as a function of the distan… | |
71 parser.add_option("-o", "--output", dest="output", type="string", | |
72 help="Output image file name (e.g. -o foo.png)") | |
73 parser.add_option("-e", "--experiment", dest="experiment", type="str… | |
74 help="MISMIP experiment: 1a,1b,2a,2b,3a,3b (e.g. -… | |
75 parser.add_option("-s", "--step", dest="step", type="int", | |
76 help="MISMIP step: 1,2,3,... (e.g. -s 1)") | |
77 parser.add_option("-m", "--model", dest="model", type="string", | |
78 help="MISMIP model (e.g. -M ABC1)") | |
79 parser.add_option("-f", "--flux", dest="profile", action="store_fals… | |
80 help="Plot ice flux only") | |
81 parser.add_option("-p", "--profile", dest="flux", action="store_fals… | |
82 help="Plot geometry profile only") | |
83 | |
84 opts, args = parser.parse_args() | |
85 | |
86 if len(args) == 0: | |
87 print("ERROR: An input file is requied.") | |
88 exit(0) | |
89 | |
90 if len(args) > 1 and opts.output: | |
91 print("More than one input file given. Ignoring the -o option...… | |
92 opts.output = None | |
93 | |
94 if opts.output and opts.profile and opts.flux: | |
95 print("Please choose between flux (-f) and profile (-p) plots.") | |
96 exit(0) | |
97 | |
98 return args, opts.output, opts | |
99 | |
100 | |
101 def read(filename, name): | |
102 "Read a variable and extract the middle row." | |
103 nc = NC(filename) | |
104 | |
105 try: | |
106 var = nc.variables[name][:] | |
107 except: | |
108 print("ERROR: Variable '%s' not present in '%s'" % (name, filena… | |
109 exit(1) | |
110 | |
111 N = len(var.shape) | |
112 if N == 1: | |
113 return var[:] # a coordinate variable ('x') | |
114 elif N == 2: | |
115 return var[1] # get the middle row | |
116 elif N == 3: | |
117 return var[-1, 1] # get the middle row of the last re… | |
118 else: | |
119 raise Exception("Can't read %s. (It's %d-dimensional.)" % (name,… | |
120 | |
121 | |
122 def find_grounding_line(x, topg, thk, mask): | |
123 "Find the modeled grounding line position." | |
124 # "positive" parts of x, topg, thk, mask | |
125 topg = topg[x > 0] | |
126 thk = thk[x > 0] | |
127 mask = mask[x > 0] | |
128 x = x[x > 0] # this should go last | |
129 | |
130 def f(j): | |
131 "See equation (7) in Pattyn et al, 'Role of transition zones in … | |
132 z_sl = 0 | |
133 return (z_sl - topg[j]) * MISMIP.rho_w() / (MISMIP.rho_i() * thk… | |
134 | |
135 for j in range(x.size): | |
136 if mask[j] == 2 and mask[j + 1] == 3: # j is grounded, j+1 floa… | |
137 nabla_f = (f(j + 1) - f(j)) / (x[j + 1] - x[j]) | |
138 | |
139 # See equation (8) in Pattyn et al | |
140 return (1.0 - f(j) + nabla_f * x[j]) / nabla_f | |
141 | |
142 raise Exception("Can't find the grounding line") | |
143 | |
144 | |
145 def plot_profile(in_file, out_file): | |
146 print("Reading %s to plot geometry profile for model %s, experiment … | |
147 in_file, model, experiment, mode, step)) | |
148 | |
149 if out_file is None: | |
150 out_file = os.path.splitext(in_file)[0] + "-profile.pdf" | |
151 | |
152 mask = read(in_file, 'mask') | |
153 usurf = read(in_file, 'usurf') | |
154 topg = read(in_file, 'topg') | |
155 thk = read(in_file, 'thk') | |
156 x = read(in_file, 'x') | |
157 | |
158 # theoretical grounding line position | |
159 xg = MISMIP.x_g(experiment, step) | |
160 # modeled grounding line position | |
161 xg_PISM = find_grounding_line(x, topg, thk, mask) | |
162 | |
163 # mask out ice-free areas | |
164 usurf = np.ma.array(usurf, mask=mask == 4) | |
165 | |
166 # compute the lower surface elevation | |
167 lsrf = topg.copy() | |
168 lsrf[mask == 3] = -MISMIP.rho_i() / MISMIP.rho_w() * thk[mask == 3] | |
169 lsrf = np.ma.array(lsrf, mask=mask == 4) | |
170 | |
171 # convert x to kilometers | |
172 x /= 1e3 | |
173 | |
174 figure(1) | |
175 ax = subplot(111) | |
176 plot(x, np.zeros_like(x), ls='dotted', color='red') | |
177 plot(x, topg, color='black') | |
178 plot(x, usurf, 'o', color='blue', markersize=4) | |
179 plot(x, lsrf, 'o', color='blue', markersize=4) | |
180 xlabel('distance from the divide, km') | |
181 ylabel('elevation, m') | |
182 title("MISMIP experiment %s, step %d" % (experiment, step)) | |
183 text(0.6, 0.9, "$x_g$ (model) = %4.0f km" % (xg_PISM / 1e3), color='… | |
184 transform=ax.transAxes) | |
185 text(0.6, 0.85, "$x_g$ (theory) = %4.0f km" % (xg / 1e3), color='bla… | |
186 transform=ax.transAxes) | |
187 | |
188 #_, _, ymin, ymax = axis(xmin=0, xmax=x.max()) | |
189 _, _, ymin, ymax = axis(xmin=x.min(), xmax=x.max()) | |
190 | |
191 vlines(xg / 1e3, ymin, ymax, linestyles='dashed', color='black') | |
192 vlines(xg_PISM / 1e3, ymin, ymax, linestyles='dashed', color='red') | |
193 | |
194 print("Saving '%s'...\n" % out_file) | |
195 savefig(out_file) | |
196 | |
197 | |
198 def plot_flux(in_file, out_file): | |
199 print("Reading %s to plot ice flux for model %s, experiment %s, grid… | |
200 in_file, model, experiment, mode, step)) | |
201 | |
202 if out_file is None: | |
203 out_file = os.path.splitext(in_file)[0] + "-flux.pdf" | |
204 | |
205 x = read(in_file, 'x') | |
206 flux_mag = read(in_file, 'flux_mag') | |
207 | |
208 # plot positive xs only | |
209 flux_mag = flux_mag[x >= 0] | |
210 x = x[x >= 0] | |
211 | |
212 figure(2) | |
213 | |
214 plot(x / 1e3, flux_mag, 'k.-', markersize=10, linewidth=2) | |
215 plot(x / 1e3, x * MISMIP.a() * MISMIP.secpera(), 'r:', linewidth=1.5) | |
216 | |
217 title("MISMIP experiment %s, step %d" % (experiment, step)) | |
218 xlabel("x ($\mathrm{km}$)", size=14) | |
219 ylabel(r"flux ($\mathrm{m}^2\,\mathrm{a}^{-1}$)", size=14) | |
220 | |
221 tight_layout() | |
222 print("Saving '%s'...\n" % out_file) | |
223 savefig(out_file, dpi=300, facecolor='w', edgecolor='w') | |
224 | |
225 | |
226 if __name__ == "__main__": | |
227 args, out_file, opts = process_options() | |
228 | |
229 for in_file in args: | |
230 model, experiment, mode, step = parse_filename(in_file, opts) | |
231 | |
232 if opts.profile: | |
233 plot_profile(in_file, out_file) | |
234 | |
235 if opts.flux: | |
236 plot_flux(in_file, out_file) |