tplot.py - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tplot.py (7378B)
---
1 #!/usr/bin/env python3
2
3 import MISMIP
4
5 from pylab import figure, subplot, hold, plot, xlabel, ylabel, title, axis, vlines, savefig, text
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 distance from the divide."
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="string",
74 help="MISMIP experiment: 1a,1b,2a,2b,3a,3b (e.g. -e 1a)")
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_false", default=True,
80 help="Plot ice flux only")
81 parser.add_option("-p", "--profile", dest="flux", action="store_false", default=True,
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...\n")
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, filename))
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 record
118 else:
119 raise Exception("Can't read %s. (It's %d-dimensional.)" % (name, N))
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 marine ice sheet dynamics', 2005."
132 z_sl = 0
133 return (z_sl - topg[j]) * MISMIP.rho_w() / (MISMIP.rho_i() * thk[j])
134
135 for j in range(x.size):
136 if mask[j] == 2 and mask[j + 1] == 3: # j is grounded, j+1 floating
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 %s, grid mode %s, step %s" % (
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 hold(True)
177 plot(x, np.zeros_like(x), ls='dotted', color='red')
178 plot(x, topg, color='black')
179 plot(x, usurf, 'o', color='blue', markersize=4)
180 plot(x, lsrf, 'o', color='blue', markersize=4)
181 xlabel('distance from the divide, km')
182 ylabel('elevation, m')
183 title("MISMIP experiment %s, step %d" % (experiment, step))
184 text(0.6, 0.9, "$x_g$ (model) = %4.0f km" % (xg_PISM / 1e3), color='r',
185 transform=ax.transAxes)
186 text(0.6, 0.85, "$x_g$ (theory) = %4.0f km" % (xg / 1e3), color='black',
187 transform=ax.transAxes)
188
189 _, _, ymin, ymax = axis(xmin=0, xmax=x.max())
190 vlines(xg / 1e3, ymin, ymax, linestyles='dashed', color='black')
191 vlines(xg_PISM / 1e3, ymin, ymax, linestyles='dashed', color='red')
192
193 print("Saving '%s'...\n" % out_file)
194 savefig(out_file)
195
196
197 def plot_flux(in_file, out_file):
198 print("Reading %s to plot ice flux for model %s, experiment %s, grid mode %s, step %s" % (
199 in_file, model, experiment, mode, step))
200
201 if out_file is None:
202 out_file = os.path.splitext(in_file)[0] + "-flux.pdf"
203
204 x = read(in_file, 'x')
205 flux_mag = read(in_file, 'flux_mag')
206
207 # plot positive xs only
208 flux_mag = flux_mag[x >= 0]
209 x = x[x >= 0]
210
211 figure(2)
212 hold(True)
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 print("Saving '%s'...\n" % out_file)
222 savefig(out_file, dpi=300, facecolor='w', edgecolor='w')
223
224
225 if __name__ == "__main__":
226 args, out_file, opts = process_options()
227
228 for in_file in args:
229 model, experiment, mode, step = parse_filename(in_file, opts)
230
231 if opts.profile:
232 plot_profile(in_file, out_file)
233
234 if opts.flux:
235 plot_flux(in_file, out_file)