GopherProxy

	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)