tshowradius.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
       ---
       tshowradius.py (1950B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 # Copyright (C) 2013, 2014, 2016, 2017 the PISM Authors
            4 
            5 # This script creates a graph of the modeled margin radius time series
            6 # by using the ice_area_glacierized variable.  Compare Figure 4(a) in Sayag & Worster (2012).
            7 # Try something like:
            8 #    ./showradius.py -o foo.png ts_lab51.nc
            9 
           10 import numpy as np
           11 import pylab as plt
           12 from sys import exit
           13 try:
           14     import netCDF4 as netCDF
           15 except:
           16     print("netCDF4 is not installed!")
           17     sys.exit(1)
           18 
           19 import argparse
           20 
           21 parser = argparse.ArgumentParser(
           22     description='A script to show margin radius time series.  Requires one or more NetCDF files with "time" dimension and "time" and "ice_area_glacierized" variables.  Generates a .png image file with the graph.')
           23 parser.add_argument('-o', '--outfile', metavar='FILENAME',
           24                     help='output file name (PNG)', default='foo.png')
           25 parser.add_argument('-d', '--datafile', metavar='FILENAME',
           26                     help='data file name (ASCII with two columns:  t r_N)')
           27 parser.add_argument('infiles', metavar='FILENAME', nargs='+',
           28                     help='input file name (NetCDF)')
           29 args = parser.parse_args()
           30 
           31 plt.figure(figsize=(12, 6))
           32 
           33 for j in range(len(args.infiles)):
           34     nc = netCDF.Dataset(args.infiles[j], "r")
           35     t = nc.variables["time"][:]
           36     ice_area_glacierized = nc.variables["ice_area_glacierized"][:]
           37     nc.close()
           38     plt.loglog(t[t > 2], np.sqrt(ice_area_glacierized[t > 2] / np.pi) * 100.0, linewidth=2.0,  # after t=2s, and in cm
           39                label=args.infiles[j])
           40 
           41 if args.datafile != None:
           42     A = np.loadtxt(args.datafile)
           43     data_t = A[:, 0]
           44     data_rN = A[:, 1]
           45     plt.loglog(data_t, 100.0 * data_rN, 'ko', label='observed', ms=4)  # cm versus s
           46 
           47 plt.legend(loc='upper left')
           48 plt.xticks([1.0, 10.0, 100.0, 1000.0])
           49 plt.yticks([1.0, 10.0])
           50 plt.axis([1.0, 1000.0, 1.0, 30.0])
           51 plt.xlabel("t  (s)", size=14)
           52 plt.ylabel(r"$r_N$  (cm)", size=14)
           53 plt.grid(True)
           54 
           55 plt.savefig(args.outfile)