tvnreport.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
---
tvnreport.py (8968B)
---
1 #!/usr/bin/env python3
2
3 from pylab import close, figure, clf, hold, plot, xlabel, ylabel, xticks, yticks, axis, legend, title, grid, show, savefig
4 from numpy import array, polyfit, polyval, log10, floor, ceil, unique
5 import sys
6
7 try:
8 from netCDF4 import Dataset as NC
9 except:
10 print("netCDF4 is not installed!")
11 sys.exit(1)
12
13
14 class Plotter:
15
16 def __init__(self, save_figures, nc, file_format):
17 self.save_figures = save_figures
18 self.nc = nc
19 self.file_format = file_format
20
21 def plot(self, x, vars, testname, plot_title):
22 # This mask lets us choose data corresponding to a particular test:
23 test = array(list(map(chr, self.nc.variables['test'][:])))
24 mask = (test == testname)
25
26 # If we have less than 2 points to plot, then bail.
27 if (sum(mask) < 2):
28 print("Skipping Test %s %s (not enough data to plot)" % (testname, plot_title))
29 return
30
31 # Get the independent variable and transform it. Note that everywhere here
32 # I assume that neither dx (dy, dz) nor errors can be zero or negative.
33 dx = self.nc.variables[x][mask]
34 dim = log10(dx)
35
36 figure(figsize=(10, 6))
37 clf()
38 hold(True)
39
40 colors = ['red', 'blue', 'green', 'black', 'brown', 'cyan']
41 for (v, c) in zip(vars, colors):
42 # Get a particular variable, transform and fit a line through it:
43 data = log10(self.nc.variables[v][mask])
44 p = polyfit(dim, data, 1)
45
46 # Try to get the long_name, use short_name if it fails:
47 try:
48 name = self.nc.variables[v].long_name
49 except:
50 name = v
51
52 # Create a label for the independent variable:
53 if (x == "dx"):
54 dim_name = "\Delta x"
55 if (x == "dy"):
56 dim_name = "\Delta y"
57 if (x == "dz"):
58 dim_name = "\Delta z"
59 if (x == "dzb"):
60 dim_name = "\Delta z_{bed}"
61
62 # Variable label:
63 var_label = "%s, $O(%s^{%1.2f})$" % (name, dim_name, p[0])
64
65 print("Test {} {}: convergence rate: O(dx^{:1.4f})".format(testname, name, p[0]))
66
67 # Plot errors and the linear fit:
68 plot(dim, data, label=var_label, marker='o', color=c)
69 plot(dim, polyval(p, dim), ls="--", color=c)
70
71 # Shrink axes, then expand vertically to have integer powers of 10:
72 axis('tight')
73 _, _, ymin, ymax = axis()
74 axis(ymin=floor(ymin), ymax=ceil(ymax))
75
76 # Switch to km if dx (dy, dz) are big:
77 units = self.nc.variables[x].units
78 if (dx.min() > 1000.0 and (units == "meters")):
79 dx = dx / 1000.0
80 units = "km"
81 # Round grid spacing in x-ticks:
82 xticks(dim, ["%d" % x for x in dx])
83 xlabel("$%s$ (%s)" % (dim_name, units))
84
85 # Use default (figured out by matplotlib) locations, but change labels for y-ticks:
86 loc, _ = yticks()
87 yticks(loc, ["$10^{%1.1f}$" % x for x in loc])
88
89 # Make sure that all variables given have the same units:
90 try:
91 ylabels = array([self.nc.variables[x].units for x in vars])
92 if (any(ylabels != ylabels[0])):
93 print("Incompatible units!")
94 else:
95 ylabel(ylabels[0])
96 except:
97 pass
98
99 # Legend, grid and the title:
100 legend(loc='best', borderpad=1, labelspacing=0.5, handletextpad=0.75, handlelength=0.02)
101 # prop = FontProperties(size='smaller'),
102 grid(True)
103 title("Test %s %s (%s)" % (testname, plot_title, self.nc.source))
104
105 if self.save_figures:
106 filename = "%s_%s_%s.%s" % (self.nc.source.replace(" ", "_"),
107 testname.replace(" ", "_"),
108 plot_title.replace(" ", "_"),
109 self.file_format)
110 savefig(filename)
111
112 def plot_tests(self, list_of_tests):
113 for test_name in list_of_tests:
114 # thickness, volume and eta errors:
115 if test_name in ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'L']:
116 self.plot('dx', ["maximum_thickness", "average_thickness"], test_name, "ice thickness errors")
117 self.plot('dx', ["relative_volume"], test_name, "relative ice volume errors")
118 self.plot('dx', ["relative_max_eta"], test_name, r"relative max eta errors")
119
120 # errors that are reported for test E only:
121 if (test_name == 'E'):
122 self.plot('dx', ["maximum_basal_velocity", "average_basal_velocity"], 'E', r"basal velocity errors")
123 self.plot('dx', ["maximum_basal_u", "maximum_basal_v"], 'E', "basal velocity (ub and vb) errors")
124 self.plot('dx', ["relative_basal_velocity"], 'E', "relative basal velocity errors")
125
126 # F and G temperature, sigma and velocity errors:
127 if test_name in ['F', 'G']:
128 self.plot('dx', ["maximum_sigma", "average_sigma"],
129 test_name, "strain heating errors")
130 self.plot('dx', ["maximum_temperature", "average_temperature",
131 "maximum_basal_temperature", "average_basal_temperature"],
132 test_name, "ice temperature errors")
133
134 self.plot('dx', ["maximum_surface_velocity", "maximum_surface_w"],
135 test_name, "maximum ice surface velocity errors")
136 self.plot('dx', ["average_surface_velocity", "average_surface_w"],
137 test_name, "average ice surface velocity errors")
138
139 # test I: plot only the u component
140 if test_name == 'I':
141 self.plot('dy', ["relative_velocity"],
142 test_name, "relative velocity errors")
143 self.plot('dy', ["maximum_u", "average_u"],
144 test_name, "velocity errors")
145
146 # tests J and M:
147 if test_name in ['J', 'M']:
148 self.plot('dx', ["relative_velocity"],
149 test_name, "relative velocity errors")
150 self.plot('dx', ["max_velocity", "maximum_u", "average_u", "maximum_v", "average_v"],
151 test_name, "velocity errors")
152
153 # test K temperature errors:
154 if (test_name == 'K'):
155 self.plot('dz', ["maximum_temperature", "average_temperature",
156 "maximum_bedrock_temperature", "average_bedrock_temperature"],
157 'K', "temperature errors")
158
159 # test O temperature and basal melt rate errors:
160 if (test_name == 'O'):
161 self.plot('dz', ["maximum_temperature", "average_temperature",
162 "maximum_bedrock_temperature", "average_bedrock_temperature"],
163 'K', "temperature errors")
164 self.plot('dz', ["maximum_basal_melt_rate"],
165 'O', "basal melt rate errors")
166
167 # test V: plot only the u component
168 if test_name == 'V':
169 self.plot('dx', ["relative_velocity"],
170 test_name, "relative velocity errors")
171 self.plot('dx', ["maximum_u", "average_u"],
172 test_name, "velocity errors")
173
174
175 from argparse import ArgumentParser
176 parser = ArgumentParser()
177 parser.description = """Plot script for PISM verification results."""
178
179 parser.add_argument("filename",
180 help="The NetCDF error report file name, usually produces by running vfnow.py")
181 parser.add_argument("-t", nargs="+", dest="tests_to_plot", default=None,
182 help="Test results to plot (space-delimited list)")
183 parser.add_argument("--save_figures", dest="save_figures", action="store_true",
184 help="Save figures to .png files")
185 parser.add_argument("--file_format", dest="file_format", default="png",
186 help="File format for --save_figures (png, pdf, jpg, ...)")
187
188 options = parser.parse_args()
189
190 input_file = NC(options.filename, 'r')
191 available_tests = unique(array(list(map(chr, input_file.variables['test'][:]))))
192 tests_to_plot = options.tests_to_plot
193
194 if len(available_tests) == 1:
195 if tests_to_plot == None:
196 tests_to_plot = available_tests
197 else:
198 if (tests_to_plot == None):
199 print("""Please choose tests to plot using the -t option.
200 (Input file %s has reports for tests %s available.)""" % (input, str(available_tests)))
201 sys.exit(0)
202
203 if (tests_to_plot[0] == "all"):
204 tests_to_plot = available_tests
205
206 close('all')
207
208 p = Plotter(options.save_figures, input_file, options.file_format)
209
210 p.plot_tests(tests_to_plot)
211 try:
212 # show() will break if we didn't plot anything
213 if not options.save_figures:
214 show()
215 except:
216 pass