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)