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()