ttauc_compare.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
       ---
       ttauc_compare.py (2434B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 try:
            4     import netCDF4 as netCDF
            5 except:
            6     print("netCDF4 is not installed!")
            7     sys.exit(1)
            8 from matplotlib import pyplot as pp
            9 from matplotlib import colors as mc
           10 from optparse import OptionParser
           11 from siple.reporting import endpause
           12 import numpy as np
           13 
           14 usage = """Usage: %prog [options]
           15 
           16 Example: %prog -N 100 -n 0.1"""
           17 
           18 parser = OptionParser(usage=usage)
           19 parser.add_option("-i", "--input_file", type='string',
           20                   help='input file')
           21 parser.add_option("-c", "--tauc_cap", type='float', default=200000,
           22                   help='maximum tauc value to display')
           23 
           24 parser.add_option("-e", "--tauc_error_cap", type='float', default=0.2,
           25                   help='maximum relative error to display')
           26 
           27 
           28 (options, args) = parser.parse_args()
           29 
           30 try:
           31     ds = netCDF.Dataset(options.input_file)
           32 except:
           33     print('ERROR: option -i is required')
           34     parser.print_help()
           35     exit(0)
           36 
           37 secpera = 3.15569259747e7
           38 tauc = ds.variables['tauc'][...].squeeze()
           39 tauc_true = ds.variables['tauc_true'][...].squeeze()
           40 
           41 tauc_diff = tauc - tauc_true
           42 
           43 not_ice = abs(ds.variables['mask'][...].squeeze() - 2) > 0.01
           44 tauc[not_ice] = 0
           45 tauc_true[not_ice] = 0
           46 tauc_diff[not_ice] = 0.
           47 
           48 
           49 u_computed = ds.variables['u_computed'][...].squeeze() * secpera
           50 v_computed = ds.variables['v_computed'][...].squeeze() * secpera
           51 velbase_mag_computed = np.sqrt(u_computed * u_computed + v_computed * v_computed)
           52 
           53 not_sliding = np.logical_and((abs(u_computed) < 10.), (abs(v_computed) < 10.))
           54 tauc[not_ice] = 0
           55 tauc_true[not_ice] = 0
           56 tauc_diff[not_sliding] = 0.
           57 
           58 
           59 # difference figure
           60 pp.clf()
           61 pp.imshow(tauc_diff.transpose() / tauc_true.transpose(), origin='lower',
           62           vmin=-options.tauc_error_cap, vmax=options.tauc_error_cap)
           63 pp.title(r'$(\tau_c$ - true) / true')
           64 pp.colorbar()
           65 
           66 # side-by-side comparison
           67 pp.figure()
           68 pp.subplot(1, 2, 1)
           69 pp.imshow(tauc.transpose(), origin='lower', vmin=0.0, vmax=options.tauc_cap)
           70 pp.title(r'$\tau_c$  [from inversion]')
           71 pp.colorbar()
           72 pp.subplot(1, 2, 2)
           73 pp.imshow(tauc_true.transpose(), origin='lower', vmin=0.0, vmax=options.tauc_cap)
           74 pp.title(r'true $\tau_c$   [prior]')
           75 pp.colorbar()
           76 
           77 # show computed sliding speed
           78 pp.figure()
           79 im = pp.imshow(velbase_mag_computed.transpose(), origin='lower',
           80                norm=mc.LogNorm(vmin=0.1, vmax=1000.0))
           81 pp.title('computed sliding speed')
           82 t = [0.1, 1.0, 10.0, 100.0, 1000.0]
           83 pp.colorbar(im, ticks=t, format='$%.1f$')
           84 
           85 # pp.ion()
           86 pp.show()
           87 # endpause()