tshowhydrovel.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
---
tshowhydrovel.py (4991B)
---
1 #!/usr/bin/env python3
2
3 from numpy import *
4 from matplotlib.pyplot import *
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(
16 description='show quiver for the subglacial water velocity (or flux) field from a PISM file')
17 parser.add_argument('filename',
18 help='file from which to get V = bwatvel[2] (and W = bwat for flux)')
19 parser.add_argument('-b', type=float, default=-1.0,
20 help='upper bound on speed; -b 100 shows all speeds above 100 as 100')
21 parser.add_argument('-c', type=float, default=-1.0,
22 help='arrow crop size; -c 0.1 shortens arrows longer than 0.1 * speed.max()')
23 # e.g. ./showhydrovel.py -c 0.001 -q extras_nbreen_y0.25_250m_routing.nc -b 20
24 parser.add_argument('-d', type=int, default=-1,
25 help='index of frame (default: last frame which is D=-1)')
26 parser.add_argument('-q', action='store_true',
27 help='show advective flux q = V W instead of velocity')
28 parser.add_argument('-s', action='store_true',
29 help='show second figure with pcolor on components of velocity (or flux)')
30 parser.add_argument('-t', action='store_true',
31 help='transpose the x,y axes')
32 parser.add_argument('-x', action='store_true',
33 help='reverse order on x-axis')
34 parser.add_argument('-y', action='store_true',
35 help='reverse order on y-axis')
36
37 args = parser.parse_args()
38
39 try:
40 nc = NC(args.filename, 'r')
41 except:
42 print("ERROR: can't read from file ...")
43 sys.exit(1)
44
45 print(" reading x,y axes from file %s ..." % (args.filename))
46 if args.t:
47 xvar = nc.variables["y"] # note x-y swap
48 yvar = nc.variables["x"]
49 else:
50 xvar = nc.variables["x"]
51 yvar = nc.variables["y"]
52 x = asarray(squeeze(xvar[:]))
53 y = asarray(squeeze(yvar[:]))
54
55 print(" reading 'bwatvel[2]' field from file %s ..." % (args.filename))
56 try:
57 velx = nc.variables["bwatvel[0]"]
58 except:
59 print("ERROR: variable 'bwatvel[0]' not found ...")
60 sys.exit(2)
61 try:
62 vely = nc.variables["bwatvel[1]"]
63 except:
64 print("ERROR: variable 'bwatvel[1]' not found ...")
65 sys.exit(3)
66
67 if args.q:
68 try:
69 bwat = nc.variables["bwat"]
70 except:
71 print("ERROR: variable 'bwat' not found ...")
72 sys.exit(6)
73
74 if args.d >= 0:
75 if shape(velx)[0] <= args.d:
76 print("ERROR: frame %d not available in variable velx" % args.d)
77 sys.exit(3)
78 if shape(vely)[0] <= args.d:
79 print("ERROR: frame %d not available in variable vely" % args.d)
80 sys.exit(4)
81 print(" reading frame %d of %d frames" % (args.d, shape(velx)[0]))
82 else:
83 args.d = -1
84 print(" reading last frame of %d frames" % (shape(velx)[0]))
85
86 units = "m hr-1" # FIXME: make this merely the default scale?
87 scale = 3.1556926e7 / 3600.0
88 velx = asarray(squeeze(velx[args.d, :, :])).transpose() / scale
89 vely = asarray(squeeze(vely[args.d, :, :])).transpose() / scale
90
91 if args.q:
92 bwat = asarray(squeeze(bwat[args.d, :, :])).transpose()
93 velx = velx * bwat
94 vely = vely * bwat
95 units = "m2 hr-1" # FIXME: adjust units?
96
97 nc.close()
98
99 if args.t:
100 xytmp = velx.copy()
101 velx = vely.transpose()
102 vely = xytmp.transpose()
103 if args.x:
104 x = x[::-1]
105 velx = -velx
106 if args.y:
107 y = y[::-1]
108 vely = -vely
109
110 if args.s:
111 figure(2)
112 print(" generating pcolor() image of velocity (or flux) components in figure(2) ...")
113 for j in [1, 2]:
114 if j == 1:
115 data = velx
116 name = "velx"
117 else:
118 data = vely
119 name = "vely"
120 print(" %s stats:\n min = %9.3f %s, max = %9.3f %s, av = %8.3f %s" %
121 (name, data.min(), units, data.max(), units, data.sum() / (x.size * y.size), units))
122 subplot(1, 2, j)
123 pcolor(x / 1000.0, y / 1000.0, data, vmin=data.min(), vmax=data.max())
124 colorbar()
125 gca().set_aspect('equal')
126 gca().autoscale(tight=True)
127 xlabel('x (km)')
128 ylabel('y (km)')
129
130 speed = sqrt(velx * velx + vely * vely)
131
132 plotvelx = velx.copy()
133 plotvely = vely.copy()
134 if args.c > 0.0:
135 crop = (speed > args.c * speed.max())
136 plotvelx[crop] = args.c * speed.max() * velx[crop] / speed[crop]
137 plotvely[crop] = args.c * speed.max() * vely[crop] / speed[crop]
138
139 if args.b > 0.0:
140 speed[speed > args.b] = args.b
141
142 figure(1)
143 quiver(x / 1000.0, y / 1000.0, plotvelx, plotvely, speed)
144 colorbar()
145 gca().set_aspect('equal')
146 gca().autoscale(tight=True)
147 xlabel('x (km)')
148 ylabel('y (km)')
149
150 if args.q:
151 print(" maximum water flux magnitude = %8.3f %s" % (speed.max(), units))
152 titlestr = "water flux in %s" % units
153 else:
154 print(" maximum water speed = %8.3f %s = %6.3f %s" %
155 (speed.max(), units, speed.max() / 3600.0, 'm s-1')) # assumes units is m hr-1
156 titlestr = "water velocity in %s" % units
157 title(titlestr)
158
159 show()
160
161 print(" done.")