tbasemapfigs.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
---
tbasemapfigs.py (6460B)
---
1 #!/usr/bin/env python3
2
3 # generate figures in Getting Started section of User's Manual
4
5 # usage:
6 # $ python basemapfigs.py FILEROOT [FIELD] [DPI]
7 # where
8 # FILEROOT root of NetCDF filename and output .png figures
9 # FIELD optional: one of {velbase_mag, [velsurf_mag], mask, usurf} (edit script to add more)
10 # DPI optional: resolution in dots per inch [200]
11 #
12 # equivalent usages:
13 # $ python basemapfigs.py g20km_10ka_hy velsurf_mag 200
14 # $ python basemapfigs.py g20km_10ka_hy velsurf_mag
15 # $ python basemapfigs.py g20km_10ka_hy
16 #
17 # generate figs like those in Getting Started section of User's Manual:
18 # $ for FLD in velsurf_mag usurf velbase_mag mask; do ./basemapfigs.py g20km_10ka_hy ${FLD}; done
19 #
20 # crop out western Greenland with command like this (uses ImageMagick):
21 # $ ./basemapfigs.py g20km_10ka_hy velsurf_mag 500
22 # $ convert -crop 600x800+400+800 +repage g20km_10ka_hy-velsurf_mag.png g20km-detail.png
23 #
24 # batch generate figures from a parameter study like this:
25 # $ for QQ in 0.1 0.5 1.0; do for EE in 1 3 6; do ../basemapfigs.py p10km_q${QQ}_e${EE} velsurf_mag 100; done; done
26 # $ for QQ in 0.1 0.5 1.0; do for EE in 1 3 6; do convert -crop 274x486+50+6 +repage p10km_q${QQ}_e${EE}-velsurf_mag.png p10km-${QQ}-${EE}-csurf.png; done; done
27
28 from mpl_toolkits.basemap import Basemap
29
30 try:
31 from netCDF4 import Dataset as NC
32 except:
33 print("netCDF4 is not installed!")
34 sys.exit(1)
35
36 import numpy as np
37 import matplotlib.pyplot as plt
38 from matplotlib import colors
39 import sys
40
41 if len(sys.argv) < 2:
42 print("ERROR: first argument must be root of filename ...")
43 sys.exit(1)
44 rootname = sys.argv[1]
45 try:
46 nc = NC(rootname + '.nc', 'r')
47 except:
48 print("ERROR: can't read from file %s.nc ..." % rootname)
49 sys.exit(2)
50
51 if len(sys.argv) >= 3:
52 field = sys.argv[2]
53 else:
54 field = 'velsurf_mag'
55
56 if len(sys.argv) >= 4:
57 mydpi = float(sys.argv[3])
58 else:
59 mydpi = 200
60
61 bluemarble = False # if True, use Blue Marble background
62
63 if (field == 'velsurf_mag') | (field == 'velbase_mag'):
64 fill = nc.variables[field]._FillValue
65 logscale = True
66 contour100 = True
67 myvmin = 1.0
68 myvmax = 6.0e3
69 ticklist = [2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000]
70 elif field == 'surfvelmag':
71 fill = 0.0
72 logscale = True
73 contour100 = True
74 myvmin = 1.0
75 myvmax = 6.0e3
76 ticklist = [2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000]
77 elif field == 'usurf':
78 fill = 0.0
79 logscale = False
80 contour100 = False
81 myvmin = 1.0
82 myvmax = 3500.0
83 ticklist = [100, 500, 1000, 1500, 2000, 2500, 3000, 3500]
84 elif field == 'mask':
85 fill = -1.0
86 logscale = False
87 contour100 = False
88 myvmin = 0.0
89 myvmax = 4.0
90 ticklist = [0, 1, 2, 3, 4]
91 elif field == 'basal_melt_rate_grounded':
92 fill = -2.0e+09
93 logscale = True
94 contour100 = False
95 myvmin = 0.9e-4
96 myvmax = 1.1
97 ticklist = [0.0001, 0.001, 0.01, 0.1, 1.0]
98 elif field == 'tillwat':
99 fill = -2.0e+09
100 logscale = False
101 contour100 = False
102 myvmin = 0.0
103 myvmax = 2.0
104 ticklist = [0.0, 0.5, 1.0, 1.5, 2.0]
105 elif field == 'bwat':
106 fill = -2.0e+09
107 logscale = True
108 contour100 = False
109 myvmin = 0.9e-4
110 myvmax = 1.1
111 ticklist = [0.0001, 0.001, 0.01, 0.1, 1.0]
112 elif field == 'bwprel':
113 fill = -2.0e+09
114 logscale = False
115 contour100 = False
116 myvmin = 0.0
117 myvmax = 1.0
118 ticklist = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
119 else:
120 print('invalid choice for FIELD option')
121 sys.exit(3)
122
123 # we need to know longitudes and latitudes corresponding to grid
124 lon = nc.variables['lon'][:]
125 lat = nc.variables['lat'][:]
126 if field == 'surfvelmag':
127 lon = np.squeeze(lon).transpose()
128 lat = np.squeeze(lat).transpose()
129
130 # x and y *in the dataset* are only used to determine plotting domain
131 # dimensions
132 if field == 'surfvelmag':
133 x = nc.variables['x1'][:]
134 y = nc.variables['y1'][:]
135 else:
136 x = nc.variables['x'][:]
137 y = nc.variables['y'][:]
138 width = x.max() - x.min()
139 height = y.max() - y.min()
140
141 # load data
142 if field == 'bwprel':
143 thkvar = np.squeeze(nc.variables['thk'][:])
144 myvar = np.squeeze(nc.variables['bwp'][:])
145 myvar = np.ma.array(myvar, mask=(thkvar == 0.0))
146 thkvar = np.ma.array(thkvar, mask=(thkvar == 0.0))
147 myvar = myvar / (910.0 * 9.81 * thkvar)
148 else:
149 myvar = np.squeeze(nc.variables[field][:])
150
151 # mask out ice free etc.; note 'mask' does not get masked
152 if (field == 'surfvelmag'):
153 myvar = myvar.transpose()
154 thkvar = np.squeeze(nc.variables['thk'][:]).transpose()
155 myvar = np.ma.array(myvar, mask=(thkvar == 0.0))
156 elif (field != 'mask'):
157 maskvar = np.squeeze(nc.variables['mask'][:])
158 if (field == 'basal_melt_rate_grounded') | (field == 'bwat'):
159 myvar[myvar < myvmin] = myvmin
160 if (field == 'usurf'):
161 myvar = np.ma.array(myvar, mask=(maskvar == 4))
162 else:
163 myvar = np.ma.array(myvar, mask=(maskvar != 2))
164
165 m = Basemap(width=1.1 * width, # width in projection coordinates, in meters
166 height=1.05 * height, # height
167 resolution='l', # coastline resolution, can be 'l' (low), 'h'
168 # (high) and 'f' (full)
169 projection='stere', # stereographic projection
170 lat_ts=71, # latitude of true scale
171 lon_0=-41, # longitude of the plotting domain center
172 lat_0=72) # latitude of the plotting domain center
173
174 # m.drawcoastlines()
175
176 # draw the Blue Marble background (requires PIL, the Python Imaging Library)
177 if bluemarble: # seems to reverse N and S
178 m.bluemarble()
179
180 # convert longitudes and latitudes to x and y:
181 xx, yy = m(lon, lat)
182
183 if contour100:
184 # mark 100 m/a contour in black:
185 m.contour(xx, yy, myvar, [100], colors="black")
186
187 # plot log color scale or not
188 if logscale:
189 m.pcolormesh(xx, yy, myvar,
190 norm=colors.LogNorm(vmin=myvmin, vmax=myvmax))
191 else:
192 m.pcolormesh(xx, yy, myvar, vmin=myvmin, vmax=myvmax)
193
194 # add a colorbar:
195 plt.colorbar(extend='both',
196 ticks=ticklist,
197 format="%d")
198
199 # draw parallels and meridians
200 # labels kwarg is where to draw ticks: [left, right, top, bottom]
201 m.drawparallels(np.arange(-55., 90., 5.), labels=[1, 0, 0, 0])
202 m.drawmeridians(np.arange(-120., 30., 10.), labels=[0, 0, 0, 1])
203
204 outname = rootname + '-' + field + '.png'
205 print("saving image to file %s ..." % outname)
206 plt.savefig(outname, dpi=mydpi, bbox_inches='tight')