tplot.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
       ---
       tplot.py (2854B)
       ---
            1 #!/usr/bin/env python3
            2 import numpy as np
            3 import matplotlib.pyplot as plt
            4 from matplotlib import colors
            5 from argparse import ArgumentParser
            6 
            7 # Import all necessary modules here so that if it fails, it fails early.
            8 try:
            9     import netCDF4 as NC
           10 except:
           11     print("netCDF4 is not installed!")
           12     sys.exit(1)
           13 
           14 # Set up the option parser
           15 parser = ArgumentParser()
           16 parser.description = "A script to plot results of the ROSS example."
           17 parser.add_argument("FILE", nargs='*')
           18 
           19 options = parser.parse_args()
           20 args = options.FILE
           21 
           22 if len(args) == 1:
           23     pism_output = args[0]
           24     plotname = pism_output.split('.')[0]
           25 else:
           26     print("wrong number of arguments, 1 expected, %i given" % int(len(args)))
           27     import sys
           28     exit(1)
           29 
           30 try:
           31     nc = NC.Dataset(pism_output, 'r')
           32 except:
           33     print("file %s not found" % pism_output)
           34     import sys
           35     exit(1)
           36 
           37 
           38 def get(name):
           39     global nc
           40     return np.squeeze(nc.variables[name][:])
           41 
           42 
           43 def floating(a):
           44     global mask
           45     return np.ma.array(a, mask=mask != 3)
           46 
           47 
           48 # load data and restrict velocities to floating areas
           49 seconds_per_year = 3.1556926e7
           50 x = get('x')
           51 y = get('y')
           52 velsurf_mag = get('velsurf_mag')
           53 mask = get('mask')
           54 u = floating(get('u_ssa'))
           55 v = floating(get('v_ssa'))
           56 # B.C.s are observations, so a PISM output file contains everything we need
           57 u_bc = floating(get('u_ssa_bc'))
           58 v_bc = floating(get('v_ssa_bc'))
           59 
           60 plt.clf()
           61 
           62 f, (a0, a1) = plt.subplots(1, 2,
           63                            gridspec_kw={'width_ratios': [1.2, 1]},
           64                            figsize=(16, 8))
           65 
           66 # mark the grounding line
           67 a0.contour(x, y, mask, [2.5], colors="black", lw=2)
           68 
           69 # plot velsurf_mag using log color scale
           70 p = a0.pcolormesh(x, y, velsurf_mag, norm=colors.LogNorm(vmin=1, vmax=1.5e3))
           71 # add a colorbar:
           72 f.colorbar(p, ax=a0, extend='both', ticks=[1, 10, 100, 500, 1000], format="%d")
           73 
           74 # quiver plot of velocities
           75 s = 10                                  # stride in grid points
           76 a0.quiver(x[::s], y[::s], u_bc[::s, ::s], v_bc[::s, ::s], color='white')
           77 a0.quiver(x[::s], y[::s], u[::s, ::s], v[::s, ::s], color='black')
           78 a0.set_xticks([])
           79 a0.set_yticks([])
           80 a0.set_title("Ross ice velocity (m/year)\nwhite=observed, black=model")
           81 
           82 # do the scatter plot
           83 magnitude = np.sqrt(np.abs(u[::s, ::s]) ** 2 + np.abs(v[::s, ::s]) ** 2)
           84 bc_magnitude = np.sqrt(np.abs(u_bc[::s, ::s]) ** 2 + np.abs(v_bc[::s, ::s]) ** 2)
           85 
           86 max_velocity = np.maximum(magnitude.max(), bc_magnitude.max())
           87 
           88 a1.scatter(magnitude, bc_magnitude, marker=".")
           89 a1.plot([0, max_velocity], [0, max_velocity], color='black', ls='--')
           90 a1.axis(xmin=0, xmax=max_velocity, ymin=0, ymax=max_velocity)
           91 a1.set_xlabel('modeled speed')
           92 a1.set_ylabel('observed speed')
           93 a1.set_title("Observed versus modeled speed (m/year)\nat points in quiver plot")
           94 
           95 output_filename = plotname + '.png'
           96 
           97 f.tight_layout()
           98 
           99 f.savefig(output_filename, dpi=300)
          100 
          101 print("saving figure %s" % output_filename)