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.")