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