GMTB Workflow Documentation
diff_GMTB_graphics.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 
3 import matplotlib as mpl
4 from array import *
5 mpl.use('Agg')
6 import pygrib
7 from mpl_toolkits.basemap import Basemap
8 import numpy as np
9 import matplotlib.pyplot as plt
10 import datetime
11 import os,sys,math
12 import maps, regrib
13 import get_grib_field as ggf
14 
15 def set_def_maps():
16 
17  map_confs = {
18  '200_temp': {'lev': '200', 'pvar': 't', 'pfill': 's', 'barbs': True, 'hgt': False, 'vc': 'isobaricInhPa'},
19  '500_hgt': {'lev': '500', 'pvar': 'gh', 'pfill': 's', 'barbs': True, 'hgt': True, 'vc': 'isobaricInhPa'},
20  '500_temp': {'lev': '500', 'pvar': 't', 'pfill': 's', 'barbs': True, 'hgt': False, 'vc': 'isobaricInhPa'},
21  '500_absv': {'lev': '500', 'pvar': 'absv', 'pfill': 's', 'barbs': False, 'hgt': False, 'vc': 'isobaricInhPa'},
22  '700_temp': {'lev': '700', 'pvar': 't', 'pfill': 's', 'barbs': True, 'hgt': False, 'vc': 'isobaricInhPa'},
23 # '700_vvel': {'lev': '700', 'pvar': 'w', 'pfill': 's', 'barbs': False, 'hgt': False, 'vc': 'isobaricInhPa'},
24  '850_hgt': {'lev': '850', 'pvar': 'gh', 'pfill': 's', 'barbs': True, 'hgt': True, 'vc': 'isobaricInhPa'},
25  '850_temp': {'lev': '850', 'pvar': 't', 'pfill': 's', 'barbs': True, 'hgt': False, 'vc': 'isobaricInhPa'},
26  '850_rh': {'lev': '850', 'pvar': 'r', 'pfill': 's', 'barbs': True, 'hgt': False, 'vc': 'isobaricInhPa'},
27  '2m_tmp': {'lev': '2', 'pvar': '2t', 'pfill': 's', 'barbs': False, 'hgt': False, 'vc': 'heightAboveGround'},
28  '2m_dpt': {'lev': '2', 'pvar': '2d', 'pfill': 's', 'barbs': False, 'hgt': False, 'vc': 'heightAboveGround'},
29  '10m_u': {'lev': '10', 'pvar': '10u', 'pfill': 's', 'barbs': False, 'hgt': False, 'vc': 'heightAboveGround'},
30  '10m_v': {'lev': '10', 'pvar': '10v', 'pfill': 's', 'barbs': False, 'hgt': False, 'vc': 'heightAboveGround'},
31  'prmsl': {'lev': '0', 'pvar': 'prmsl', 'pfill': 's', 'pline': True, 'barbs': True, 'hgt': False, 'vc': 'meanSea'},
32  'pwat': {'lev': '0', 'pvar': 'pwat', 'pfill': 's', 'pline': False, 'barbs': False, 'hgt': False, 'vc': 'unknown'},
33  'prate': {'lev': '0', 'pvar': 'prate','pfill': 's', 'pline': False, 'barbs': False, 'hgt': False, 'vc': 'surface', 'acc_field': True},
34  'cprat': {'lev': '0', 'pvar': 'cprat','pfill': 's', 'pline': False, 'barbs': False, 'hgt': False, 'vc': 'surface', 'acc_field': True},
35  '6h_acpcp': {'lev': '0', 'pvar': 'acpcp', 'pfill': 's', 'barbs': False, 'hgt': False, 'vc': 'surface', 'acc_field': True},
36  '6h_apcp': {'lev': '0', 'pvar': 'tp', 'pfill': 's', 'barbs': False, 'hgt': False, 'vc': 'surface', 'acc_field': True}
37 # '6h_ncpcp': {'lev': '0', 'pvar': 'ncpcp', 'pfill': 's', 'barbs': False, 'hgt': False, 'vc': 'surface'}
38  }
39 
40  return map_confs
41 
42 
43 def main():
44  # Get the cycle information from env vars
45  moad_dataroot=os.getenv('MOAD_DATAROOT')
46  rotdir=os.getenv('ROTDIR')
47  expdir=os.getenv('EXPDIR')
48  rotnm=os.getenv('ROTNM')
49  expnm=os.getenv('EXPNM')
50  print
51  print 'plotting the difference between:',expnm,'and',rotnm
52  diffstr=''.join(['(',expnm,'-',rotnm,')'])
53  print diffstr
54  print
55  date=os.getenv('START_TIME') # @Y@m@d@H
56  fhr=os.getenv('FCST_TIME') # 3 digit fcst hour
57  grid_list=os.getenv('GRID_VX_LIST')
58  vx_grid=grid_list.split(" ")
59  filename=''.join(['pgrbq', fhr, '.gfs.', date, '_', grid_list, '.grib2'])
60  rotf='/'.join([rotdir, filename])
61  expf='/'.join([expdir, filename])
62  print ' rotdir = ', rotdir
63  print ' expdir = ', expdir
64  print
65  print '... plotting the difference fields ...'
66  user_def = os.getenv('USER_MAPS',False)
67 
68  # Set the output file location (template set in maps.py).
69  #outdir='/'.join([moad_dataroot, 'figprd'])
70  outdir='/'.join(['./', 'figprd'])
71 
72  for grid in vx_grid:
73 
74  print 'grid = ', grid
75  outmap = grid
76  plot_wind=True
77 
78  def_maps=set_def_maps()
79  for fig in def_maps:
80  plt.figure(figsize=(12,8))
81  # Optional if statement here can control evaluating one type of figure at a time.
82  # if fig == '6h_apcp' or fig == '6h_acpcp':
83  acc_field=def_maps[fig].get('acc_field',False)
84  # avg_field=def_maps[fig].get('avg_field',False)
85 
86  if (acc_field and fhr != '000') or (not acc_field):
87  level = def_maps[fig].get('lev')
88  var = def_maps[fig].get('pvar')
89  hgt = def_maps[fig].get('hgt')
90  plot_wind = def_maps[fig].get('barbs',False)
91  plot_line = def_maps[fig].get('pline',False)
92  coord = def_maps[fig].get('vc')
93  #print level,var #, hgt,plot_wind, hgt
94 
95  # Grab winds from file if plotting the barbs
96  if plot_wind and var == 'wind':
97  u = ggf.get_ua_field(rotf,level,coord,shortName='u')
98  v = ggf.get_ua_field(rotf,level,coord,shortName='v')
99  winds=[u,v]
100  elif plot_wind and var == '10wind':
101  u = ggf.get_ua_field(rotf,level,coord,shortName='10u')
102  v = ggf.get_ua_field(rotf,level,coord,shortName='10v')
103  winds=[u,v]
104  else:
105  winds=None
106 
107  height = None
108  if hgt:
109  height = ggf.get_ua_field(rotf,level,coord,shortName='gh')
110 
111  # Get the variable from different file
112  fd_rot=ggf.get_ua_field(rotf,level,coord,shortName=var)
113  fd_exp=ggf.get_ua_field(expf,level,coord,shortName=var)
114  atime = fd_rot.analDate.strftime('Analysis: %Y%m%d %H UTC')
115  altime= fd_rot.analDate.strftime('%Y%m%d%H')
116  vtime = fd_rot.validDate.strftime('Valid: %Y%m%d %H UTC')
117 
118  fdrot = fd_rot.values
119  fdexp = fd_exp.values
120 
121 
122  data = fdexp - fdrot
123  lat,lon = fd_rot.latlons()
124 
125  if grid == 'glob' or grid == 'G3':
126  m =Basemap(projection='cyl',lat_ts=10,llcrnrlon=lon.min(), \
127  urcrnrlon=lon.max(),llcrnrlat=lat.min(),urcrnrlat=lat.max(), \
128  resolution='c')
129  elif grid == 'nh_mill':
130  m =Basemap(projection='mill',lat_ts=10,llcrnrlon=lon.min(), \
131  urcrnrlon=lon.max(),llcrnrlat=0,urcrnrlat=90, \
132  resolution='c')
133  elif grid == 'sh_mill':
134  m =Basemap(projection='mill',lat_ts=10,llcrnrlon=lon.min(), \
135  urcrnrlon=lon.max(),llcrnrlat=-90,urcrnrlat=0, \
136  resolution='c')
137  elif grid == 'trop_mill':
138  m =Basemap(projection='mill',lat_ts=10,llcrnrlon=lon.min(), \
139  urcrnrlon=lon.max(),llcrnrlat=-23,urcrnrlat=23, \
140  resolution='c')
141  elif grid == 'nps':
142  m =Basemap(projection='npstere',boundinglat=-0.268, \
143  lon_0=-105,resolution='c')
144  elif grid == 'sps':
145  m =Basemap(projection='npstere',boundinglat=-60, \
146  lon_0=105,resolution='c')
147  elif grid == 'G218':
148  m =Basemap(projection='lcc',lat_ts=10,llcrnrlon=226.514, \
149  urcrnrlon=-49.420,llcrnrlat=12.190,urcrnrlat=57.328, \
150  resolution='c',lat_1=25,lat_2=25,lon_0=265,rsphere=6371200)
151  elif grid == 'G104':
152  lon_w,lon_e,lat_s,lat_n =[-147, -64, 10.8, 89.9]
153  lon_0 = lon_w + (lon_e - lon_w) / 2.
154  ref = lat_s if abs(lat_s) > abs(lat_n) else lat_n
155  lat_0 = math.copysign(90., ref)
156  proj = 'npstere' if lat_0 > 0 else 'spstere'
157  prj = Basemap(projection=proj, lon_0=lon_0, lat_0=lat_0,
158  boundinglat=0, resolution='c')
159  lons = [lon_w, lon_e, lon_w, lon_e, lon_0, lon_0]
160  lats = [lat_s, lat_s, lat_n, lat_n, lat_s, lat_n]
161  x, y = prj(lons, lats)
162  ll_lon, ll_lat = prj(min(x), min(y), inverse=True)
163  ur_lon, ur_lat = prj(max(x), max(y), inverse=True)
164  m =Basemap(projection='stere', lat_0=lat_0, lon_0=lon_0,
165  llcrnrlon=ll_lon, llcrnrlat=ll_lat,
166  urcrnrlon=ur_lon, urcrnrlat=ur_lat, resolution='c')
167 
168  x, y = m(lon,lat)
169 
170  # Draw political boundaries
171  m.drawcoastlines()
172  m.drawstates()
173  m.drawcountries()
174  m.drawmapboundary()
175  m.drawlsmask(ocean_color='white',land_color='white')
176 
177  # Draw lat/lon grid
178  parallels = np.arange(-90.,90,10.)
179  m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
180  meridians = np.arange(0.,360.,20.)
181  m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
182 
183  #clevs=np.arange(-7,7,1)
184  # adjust the values based on exp. performance
185  if var == 'pwat':
186  cbd= 50
187  itv= 5
188  elif var == 'prate':
189  cbd=0.004
190  itv=0.001
191  elif var == 'absv':
192  data=data*10000 #scale
193  cbd=8
194  itv=1
195  elif var == '2d':
196  cbd= 40
197  itv= 5
198  elif var == 'r':
199  cbd=120
200  itv=20
201  elif var == 't' or var=='2t':
202  cbd=40
203  itv=5
204  elif var == 'acpcp':
205  cbd=25
206  itv=5
207  elif var == '10u' or var== '10v':
208  cbd=50
209  itv=10
210  elif var == 'cprat':
211  cbd=0.002
212  itv=0.0004
213  elif var == 'gh':
214  cbd=400
215  itv=50
216  elif var == 'tp':
217  cbd=80
218  itv=10
219  elif var =='prmsl':
220  cbd=5000
221  itv=500
222 
223  clevs=np.arange(-cbd,cbd,itv)
224 
225  #cs = m.contourf(x,y,data,cmap=plt.cm.RdBu_r)
226  cs = m.contourf(x,y,data,clevs,cmap=plt.cm.RdBu_r)
227  if grid == 'G218':
228  csl = m.contour(x,y,data,clevs,colors='k')
229 
230  # Determine how to label the vertical coordinate based on the contents of the grib file, or other
231  # information.'''
232  if coord == 'isobaricInhPa':
233  vert_unit = 'hPa'
234  elif coord == 'hybrid':
235  vert_unit = 'sigma'
236  elif coord == 'heightAboveGround':
237  vert_unit = 'm'
238  else:
239  vert_unit ='m'
240 
241  print 'diff ', var,',', level, vert_unit, ': max/min = ',data.max(), data.min(),'; clevs = ',clevs
242 
243 
244  #Add a colorbar and title, and then show the plot.
245  plt.colorbar(cs,orientation='horizontal',extend='both')
246 
247  #
248  maptype = 'shaded'
249  myvar = var
250  title_str=str('%s %s' %(myvar,diffstr))
251  plt.title('%s \nFcst Hr: %s' % (atime, fhr) , loc='left')
252  plt.title('%s %s' % (myvar,diffstr), loc = 'right')
253  plt.title('%s %s' % (level,vert_unit))
254  plt.xlabel('%s' % (vtime), fontsize=18, labelpad=40)
255 
256  # Output file name definition
257  filename=''.join(['%s' % (myvar),'%s' % (level),vert_unit,'_',grid,'_',altime,'_f',fhr,'_diff.png'])
258 
259  # Check to see if output location exists, if not, make it.
260  if not os.path.exists(outdir):
261  os.makedirs(outdir)
262 
263  plt.savefig('/'.join([outdir, filename]))
264  plt.close()
265 
266 main()
diff_GMTB_graphics.main
def main()
Definition: diff_GMTB_graphics.py:43
diff_GMTB_graphics.set_def_maps
def set_def_maps()
Definition: diff_GMTB_graphics.py:15