tshowhydro.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 --- tshowhydro.py (2227B) --- 1 #!/usr/bin/env python3 2 3 from numpy import * 4 from matplotlib.pyplot import * 5 from sys import exit 6 7 try: 8 from netCDF4 import Dataset as NC 9 except: 10 print("netCDF4 is not installed!") 11 sys.exit(1) 12 13 nameroot = "routing" 14 15 for dx in ("100", "50", "25", "15", "10", "5"): 16 basename = nameroot + dx + "km" 17 filename = basename + ".nc" 18 print("%s: looking for file ..." % filename) 19 try: 20 nc = NC(filename, 'r') 21 except: 22 print(" can't read from file ...") 23 continue 24 25 xvar = nc.variables["x"] 26 yvar = nc.variables["y"] 27 x = asarray(squeeze(xvar[:])) 28 y = asarray(squeeze(yvar[:])) 29 30 for varname in ("bwat", "bwp", "psi"): # psi must go after bwat, bwp 31 print(" %s: generating pcolor() image ..." % varname) 32 try: 33 if varname == "psi": 34 var = nc.variables["topg"] 35 else: 36 var = nc.variables[varname] 37 except: 38 print("variable '%s' not found ... continuing ..." % varname) 39 continue 40 41 data = asarray(squeeze(var[:])).transpose() 42 43 if varname == "bwat": 44 bwatdata = data.copy() 45 if varname == "bwp": 46 bwpdata = data.copy() 47 48 if varname == "psi": 49 # psi = bwp + rho_w g (topg + bwat) 50 data = bwpdata + 1000.0 * 9.81 * (data + bwatdata) 51 52 if varname == "bwat": 53 units = "m" 54 barmin = 0.0 55 barmax = 650.0 56 scale = 1.0 57 else: 58 units = "bar" 59 barmin = -20.0 60 barmax = 360.0 61 scale = 1.0e5 62 63 print(" [stats: max = %9.3f %s, av = %8.3f %s]" % 64 (data.max() / scale, units, data.sum() / (scale * x.size * y.size), units)) 65 pcolor(x / 1000.0, y / 1000.0, data / scale, vmin=barmin, vmax=barmax) 66 colorbar() 67 gca().set_aspect('equal') 68 gca().autoscale(tight=True) 69 xlabel('x (km)') 70 ylabel('y (km)') 71 dxpad = "%03d" % int(dx) 72 pngfilename = varname + "_" + dxpad + "km" + ".png" 73 print(" saving figure in %s ..." % pngfilename) 74 savefig(pngfilename, dpi=300, bbox_inches='tight') 75 close() 76 77 nc.close()