tshowPvsW.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 --- tshowPvsW.py (6200B) --- 1 #!/usr/bin/env python3 2 3 import numpy as np 4 import matplotlib.pyplot as plt 5 6 import sys 7 import argparse 8 9 try: 10 from netCDF4 import Dataset as NC 11 except: 12 print("netCDF4 is not installed!") 13 sys.exit(1) 14 15 parser = argparse.ArgumentParser(description='show scatter plot P versus W from a PISM run') 16 17 parser.add_argument('filename', 18 help='file from which to get P = bwprel and W = bwat') 19 parser.add_argument('-o', default=None, 20 help='output file for image, in a format matplotlib can write (e.g. .png, .pdf)') 21 parser.add_argument('-d', type=int, default=-1, 22 help='index of frame (default: last frame which is D=-1)') 23 24 parser.add_argument('-c', default=None, 25 help='name of variable to use to color points in scatter plot') 26 parser.add_argument('-cmin', type=float, default=None, 27 help='crop color values below at this value') 28 parser.add_argument('-cmax', type=float, default=None, 29 help='crop color values above at this value') 30 parser.add_argument('--colorbar', action='store_true', 31 help='whether to show colorbar() beside figure') 32 33 parser.add_argument('-s', default=None, 34 help='name of variable to use to select whether points appear in scatter plot') 35 parser.add_argument('-smin', type=float, default=None, 36 help='select minimum: if -c is used then below this value (of -s var) the points will not be plotted') 37 parser.add_argument('-smax', type=float, default=None, 38 help='select minimum: if -c is used then above this value (of -s var) the points will not be plotted') 39 40 parser.add_argument('-wmin', type=float, default=None, 41 help='lower limit on W axis') 42 parser.add_argument('-wmax', type=float, default=None, 43 help='upper limit on W axis') 44 45 args = parser.parse_args() 46 47 try: 48 nc = NC(args.filename, 'r') 49 except: 50 print("ERROR: can't read from file ...") 51 sys.exit(1) 52 53 print(" reading 'bwprel' field from file %s ..." % (args.filename)) 54 try: 55 bwprel = nc.variables["bwprel"] 56 except: 57 print("ERROR: variable 'bwprel' not found ...") 58 sys.exit(2) 59 60 print(" reading 'bwat' field from file %s ..." % (args.filename)) 61 try: 62 bwat = nc.variables["bwat"] 63 except: 64 print("ERROR: variable 'bwat' not found ...") 65 sys.exit(3) 66 67 if args.s != None: 68 print(" reading '%s' field, for point selection, from file %s ..." % (args.s, args.filename)) 69 try: 70 ss = nc.variables[args.s] 71 except: 72 print("ERROR: variable '%s' not found ..." % (args.s)) 73 sys.exit(5) 74 if args.smin != None: 75 print(" minimum value of '%s' for selection is %f" % (args.s, args.smin)) 76 if args.smax != None: 77 print(" maximum value of '%s' for selection is %f" % (args.s, args.smax)) 78 79 if args.c != None: 80 print(" reading '%s' field, for point color, from file %s ..." % (args.c, args.filename)) 81 try: 82 cc = nc.variables[args.c] 83 except: 84 print("ERROR: variable '%s' not found ..." % (args.c)) 85 sys.exit(4) 86 if args.cmin != None: 87 print(' color minimum value is %f' % args.cmin) 88 if args.cmax != None: 89 print(' color maximum value is %f' % args.cmax) 90 91 if args.d >= 0: 92 if np.shape(bwprel)[0] <= args.d: 93 print("ERROR: frame %d not available in variable bwprel" % args.d) 94 sys.exit(4) 95 if np.shape(bwat)[0] <= args.d: 96 print("ERROR: frame %d not available in variable bwat" % args.d) 97 sys.exit(5) 98 print(" using frame %d of %d frames" % (args.d, np.shape(bwat)[0])) 99 else: 100 args.d = -1 101 print(" reading last frame of %d frames" % (np.shape(bwat)[0])) 102 103 bwat = np.asarray(np.squeeze(bwat[args.d, :, :])).flatten() 104 bwprel = np.asarray(np.squeeze(bwprel[args.d, :, :])).flatten() 105 106 bwprel[bwprel > 1.0] = 1.0 107 bwprel[bwprel < 0.0] = 0.0 108 109 if args.c != None: 110 ccc = np.asarray(np.squeeze(cc[args.d, :, :])).flatten() 111 if args.cmin != None: 112 ccc[ccc < args.cmin] = args.cmin 113 if args.cmax != None: 114 ccc[ccc > args.cmax] = args.cmax 115 116 if args.s != None: 117 sss = np.asarray(np.squeeze(ss[args.d, :, :])).flatten() 118 if args.smin != None: 119 bwat = bwat[sss >= args.smin] 120 bwprel = bwprel[sss >= args.smin] 121 if args.c != None: 122 ccc = ccc[sss >= args.smin] 123 sss = sss[sss >= args.smin] 124 if args.smax != None: 125 bwat = bwat[sss <= args.smax] 126 bwprel = bwprel[sss <= args.smax] 127 if args.c != None: 128 ccc = ccc[sss <= args.smax] 129 sss = sss[sss <= args.smax] 130 131 nc.close() 132 133 # to reduce file size, remove zero water points 134 if args.c != None: 135 ccc = ccc[bwat > 0] 136 bwprel = bwprel[bwat > 0.0] 137 bwat = bwat[bwat > 0.0] 138 139 # to reduce file size, remove zero pressure points 140 if args.c != None: 141 ccc = ccc[bwprel > 0] 142 bwat = bwat[bwprel > 0.0] 143 bwprel = bwprel[bwprel > 0.0] 144 145 print(" scatter plot has %d points ...\n" % len(bwat)) 146 147 # W axis limits 148 if args.wmin == None: 149 wwmin = min(bwat) 150 else: 151 wwmin = args.wmin 152 if args.wmax == None: 153 wwmax = max(bwat) 154 else: 155 wwmax = args.wmax 156 157 # color axis limits 158 if args.c != None: 159 if args.cmin == None: 160 ccmin = min(ccc) 161 else: 162 ccmin = args.cmin 163 if args.cmax == None: 164 ccmax = max(ccc) 165 else: 166 ccmax = args.cmax 167 168 import matplotlib 169 if args.colorbar: 170 matplotlib.rcParams['figure.figsize'] = (3.8, 3.0) 171 else: 172 matplotlib.rcParams['figure.figsize'] = (3.0, 3.0) 173 174 plt.figure(1) 175 if args.c != None: 176 plt.scatter(bwat, bwprel, c=ccc, vmin=ccmin, vmax=ccmax, linewidth=0.0, cmap='hsv') 177 else: 178 plt.scatter(bwat, bwprel, c='k') 179 plt.axis('tight') 180 plt.gca().set_xlim((wwmin, wwmax)) 181 plt.gca().set_ylim((-0.02, 1.05)) 182 plt.gca().set_xticks((0.0, 0.05, 0.10, 0.15)) 183 plt.gca().set_xticklabels(('0', '.05', '.10', '.15')) 184 plt.text(0.07, 0.25, '%d' % int(args.smin) + r'$ < |\mathbf{v}_b| < $' + '%d' % int(args.smax)) 185 if args.colorbar: 186 plt.colorbar() 187 plt.xlabel(r'$W$ (m)') 188 plt.ylabel(r'$P/P_0$') 189 190 if args.o == None: 191 plt.show() 192 else: 193 print(" saving scatter plot in %s ...\n" % args.o) 194 plt.savefig(args.o, dpi=300, bbox_inches='tight')