Introduction
Introduction Statistics Contact Development Disclaimer Help
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)
You are viewing proxied material from mx1.adamsgaard.dk. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.