tshowPvsW.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
---
tshowPvsW.py (6200B)
---
1 #!/usr/bin/env python3
2
3 import numpy as np
4 import matplotlib.pyplot as plt
5
6 import sys
7 import argparse
8
9 try:
10 from netCDF4 import Dataset as NC
11 except:
12 print("netCDF4 is not installed!")
13 sys.exit(1)
14
15 parser = argparse.ArgumentParser(description='show scatter plot P versus W from a PISM run')
16
17 parser.add_argument('filename',
18 help='file from which to get P = bwprel and W = bwat')
19 parser.add_argument('-o', default=None,
20 help='output file for image, in a format matplotlib can write (e.g. .png, .pdf)')
21 parser.add_argument('-d', type=int, default=-1,
22 help='index of frame (default: last frame which is D=-1)')
23
24 parser.add_argument('-c', default=None,
25 help='name of variable to use to color points in scatter plot')
26 parser.add_argument('-cmin', type=float, default=None,
27 help='crop color values below at this value')
28 parser.add_argument('-cmax', type=float, default=None,
29 help='crop color values above at this value')
30 parser.add_argument('--colorbar', action='store_true',
31 help='whether to show colorbar() beside figure')
32
33 parser.add_argument('-s', default=None,
34 help='name of variable to use to select whether points appear in scatter plot')
35 parser.add_argument('-smin', type=float, default=None,
36 help='select minimum: if -c is used then below this value (of -s var) the points will not be plotted')
37 parser.add_argument('-smax', type=float, default=None,
38 help='select minimum: if -c is used then above this value (of -s var) the points will not be plotted')
39
40 parser.add_argument('-wmin', type=float, default=None,
41 help='lower limit on W axis')
42 parser.add_argument('-wmax', type=float, default=None,
43 help='upper limit on W axis')
44
45 args = parser.parse_args()
46
47 try:
48 nc = NC(args.filename, 'r')
49 except:
50 print("ERROR: can't read from file ...")
51 sys.exit(1)
52
53 print(" reading 'bwprel' field from file %s ..." % (args.filename))
54 try:
55 bwprel = nc.variables["bwprel"]
56 except:
57 print("ERROR: variable 'bwprel' not found ...")
58 sys.exit(2)
59
60 print(" reading 'bwat' field from file %s ..." % (args.filename))
61 try:
62 bwat = nc.variables["bwat"]
63 except:
64 print("ERROR: variable 'bwat' not found ...")
65 sys.exit(3)
66
67 if args.s != None:
68 print(" reading '%s' field, for point selection, from file %s ..." % (args.s, args.filename))
69 try:
70 ss = nc.variables[args.s]
71 except:
72 print("ERROR: variable '%s' not found ..." % (args.s))
73 sys.exit(5)
74 if args.smin != None:
75 print(" minimum value of '%s' for selection is %f" % (args.s, args.smin))
76 if args.smax != None:
77 print(" maximum value of '%s' for selection is %f" % (args.s, args.smax))
78
79 if args.c != None:
80 print(" reading '%s' field, for point color, from file %s ..." % (args.c, args.filename))
81 try:
82 cc = nc.variables[args.c]
83 except:
84 print("ERROR: variable '%s' not found ..." % (args.c))
85 sys.exit(4)
86 if args.cmin != None:
87 print(' color minimum value is %f' % args.cmin)
88 if args.cmax != None:
89 print(' color maximum value is %f' % args.cmax)
90
91 if args.d >= 0:
92 if np.shape(bwprel)[0] <= args.d:
93 print("ERROR: frame %d not available in variable bwprel" % args.d)
94 sys.exit(4)
95 if np.shape(bwat)[0] <= args.d:
96 print("ERROR: frame %d not available in variable bwat" % args.d)
97 sys.exit(5)
98 print(" using frame %d of %d frames" % (args.d, np.shape(bwat)[0]))
99 else:
100 args.d = -1
101 print(" reading last frame of %d frames" % (np.shape(bwat)[0]))
102
103 bwat = np.asarray(np.squeeze(bwat[args.d, :, :])).flatten()
104 bwprel = np.asarray(np.squeeze(bwprel[args.d, :, :])).flatten()
105
106 bwprel[bwprel > 1.0] = 1.0
107 bwprel[bwprel < 0.0] = 0.0
108
109 if args.c != None:
110 ccc = np.asarray(np.squeeze(cc[args.d, :, :])).flatten()
111 if args.cmin != None:
112 ccc[ccc < args.cmin] = args.cmin
113 if args.cmax != None:
114 ccc[ccc > args.cmax] = args.cmax
115
116 if args.s != None:
117 sss = np.asarray(np.squeeze(ss[args.d, :, :])).flatten()
118 if args.smin != None:
119 bwat = bwat[sss >= args.smin]
120 bwprel = bwprel[sss >= args.smin]
121 if args.c != None:
122 ccc = ccc[sss >= args.smin]
123 sss = sss[sss >= args.smin]
124 if args.smax != None:
125 bwat = bwat[sss <= args.smax]
126 bwprel = bwprel[sss <= args.smax]
127 if args.c != None:
128 ccc = ccc[sss <= args.smax]
129 sss = sss[sss <= args.smax]
130
131 nc.close()
132
133 # to reduce file size, remove zero water points
134 if args.c != None:
135 ccc = ccc[bwat > 0]
136 bwprel = bwprel[bwat > 0.0]
137 bwat = bwat[bwat > 0.0]
138
139 # to reduce file size, remove zero pressure points
140 if args.c != None:
141 ccc = ccc[bwprel > 0]
142 bwat = bwat[bwprel > 0.0]
143 bwprel = bwprel[bwprel > 0.0]
144
145 print(" scatter plot has %d points ...\n" % len(bwat))
146
147 # W axis limits
148 if args.wmin == None:
149 wwmin = min(bwat)
150 else:
151 wwmin = args.wmin
152 if args.wmax == None:
153 wwmax = max(bwat)
154 else:
155 wwmax = args.wmax
156
157 # color axis limits
158 if args.c != None:
159 if args.cmin == None:
160 ccmin = min(ccc)
161 else:
162 ccmin = args.cmin
163 if args.cmax == None:
164 ccmax = max(ccc)
165 else:
166 ccmax = args.cmax
167
168 import matplotlib
169 if args.colorbar:
170 matplotlib.rcParams['figure.figsize'] = (3.8, 3.0)
171 else:
172 matplotlib.rcParams['figure.figsize'] = (3.0, 3.0)
173
174 plt.figure(1)
175 if args.c != None:
176 plt.scatter(bwat, bwprel, c=ccc, vmin=ccmin, vmax=ccmax, linewidth=0.0, cmap='hsv')
177 else:
178 plt.scatter(bwat, bwprel, c='k')
179 plt.axis('tight')
180 plt.gca().set_xlim((wwmin, wwmax))
181 plt.gca().set_ylim((-0.02, 1.05))
182 plt.gca().set_xticks((0.0, 0.05, 0.10, 0.15))
183 plt.gca().set_xticklabels(('0', '.05', '.10', '.15'))
184 plt.text(0.07, 0.25, '%d' % int(args.smin) + r'$ < |\mathbf{v}_b| < $' + '%d' % int(args.smax))
185 if args.colorbar:
186 plt.colorbar()
187 plt.xlabel(r'$W$ (m)')
188 plt.ylabel(r'$P/P_0$')
189
190 if args.o == None:
191 plt.show()
192 else:
193 print(" saving scatter plot in %s ...\n" % args.o)
194 plt.savefig(args.o, dpi=300, bbox_inches='tight')