plot.py (7433B)
1 #!/usr/bin/env python3 2 3 import MISMIP 4 5 from pylab import figure, subplot, plot, xlabel, ylabel, title, axis, vlines, savefig, text, tight_layout 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 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='r', 184 transform=ax.transAxes) 185 text(0.6, 0.85, "$x_g$ (theory) = %4.0f km" % (xg / 1e3), color='black', 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 mode %s, step %s" % ( 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)