tplot_flowline_results.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_flowline_results.py (3571B)
---
1 #!/usr/bin/env python3
2
3 try:
4 from netCDF4 import Dataset
5 except:
6 print("netCDF4 is not installed!")
7 sys.exit(1)
8
9 import numpy as np
10 import pylab as plt
11
12 from optparse import OptionParser
13
14 parser = OptionParser()
15 parser.usage = "usage: %prog [options] FILE"
16 parser.description = "A script to compare PISM flowline velocities with Stokes solution."
17
18 (options, args) = parser.parse_args()
19
20 plot_acab = True
21
22 if len(args) != 1:
23 print('wrong number of arguments, 1 expected')
24 exit(1)
25
26 try:
27 nc = Dataset(args[0], 'r')
28 except:
29 print("file %s not found ... ending ..." % args[0])
30 exit(2)
31
32
33 def permute(variable, output_order=('time', 'z', 'zb', 'y', 'x')):
34 """Permute dimensions of a NetCDF variable to match the output storage order."""
35 input_dimensions = variable.dimensions
36
37 # filter out irrelevant dimensions
38 dimensions = [x for x in output_order if x in input_dimensions]
39
40 # create the mapping
41 mapping = [dimensions.index(x) for x in input_dimensions]
42
43 if mapping:
44 return np.transpose(variable[:], mapping)
45 else:
46 return variable[:] # so that it does not break processing "mapping"
47
48
49 x = nc.variables["x"][:]
50 b = np.squeeze(nc.variables["topg"][:])
51 s = np.squeeze(nc.variables["usurf"][:])
52 h = np.squeeze(nc.variables["thk"][:])
53 z = nc.variables["z"][:]
54
55 mask = h <= 1
56
57 us = np.ma.array(data=np.squeeze(nc.variables["uvelsurf"][:]), mask=mask)
58 ub = np.ma.array(data=np.squeeze(nc.variables["uvelbase"][:]), mask=mask)
59 # stuff needed for contour plots
60 xx = (np.tile(x, [len(z), 1]))
61 zz = ((np.tile(z, [len(x), 1])).transpose() + b)
62
63 # ignore the first level
64 cts = np.squeeze(permute(nc.variables["cts"]))
65 liqfrac = np.squeeze(permute(nc.variables["liqfrac"]))
66 temppa = np.squeeze(permute(nc.variables["temp_pa"]))
67 mask2 = np.zeros_like(cts)
68 mask2[zz > s] = 1
69 cts = np.ma.array(data=cts, mask=mask2)
70 liqfrac = np.ma.array(data=liqfrac, mask=mask2)
71 temppa = np.ma.array(data=temppa, mask=mask2)
72
73 # Contour level of the CTS
74 cts_level = [1]
75 liqfrac_levels = np.arange(0, 2.5, .25)
76 temppa_levels = [-6, -5, -4, -3, -2, -1, -.0001]
77
78 fig = plt.figure(figsize=(6.4, 7.4))
79
80 axUpperLeft = plt.axes([0.1, 0.6, 0.8, 0.25])
81 axLower = plt.axes([0.1, 0.05, 0.8, 0.5])
82
83 axUpperLeft.plot(x, us, color='#377EB8', lw=1.5)
84 axUpperLeft.plot(x, ub, '--', color='#377EB8', lw=1.5)
85 axUpperLeft.axes.set_xlim(-250, 3500)
86 axUpperLeft.axes.set_ylabel("velocity [m a$^{-1}$]")
87 plt.setp(axUpperLeft, xticks=[])
88 if (plot_acab == True):
89 acab = np.squeeze(nc.variables["climatic_mass_balance"][:])
90 axUpperRight = axUpperLeft.twinx()
91 axUpperRight.plot(x, acab / 910.0, color='#984EA3', lw=1.5)
92 axUpperRight.axes.set_ylabel("mass balance [m a$^{-1}$]")
93 axUpperRight.axes.set_xlim(-250, 3500)
94 axLower.plot(x, b, color='black', lw=1.5)
95 axLower.plot(x, s, color='black', lw=1.5)
96 c1 = axLower.contourf(xx, zz, liqfrac * 100, liqfrac_levels, cmap=plt.cm.Reds)
97 plt.colorbar(mappable=c1, ax=axLower, orientation='horizontal', pad=0.05, shrink=0.75, extend="max")
98 c2 = axLower.contourf(xx, zz, temppa, temppa_levels, cmap=plt.cm.Blues_r)
99 plt.colorbar(mappable=c2, ax=axLower, orientation='horizontal',
100 ticks=[-6, -5, -4, -3, -2, -1, 0], pad=0.20, shrink=0.75)
101 axLower.contour(xx, zz, cts, cts_level, colors='black', linestyles='dashed')
102 axLower.axes.set_xlim(-250, 3500)
103 axLower.axes.set_ylim(1100, 1800)
104 axLower.axes.set_xlabel("distance from bergschrund [m]")
105 axLower.axes.set_ylabel("elevation [m a.s.l.]")
106
107
108 plt.savefig('sg_results.pdf', bbox_inches='tight', pad_inches=0.35)
109
110 nc.close()