/usr/rmt_share/baldwin/EAS591/plotpcp.py
#!/usr/bin/env python
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
from optparse import OptionParser
from pylab import *
import Nio
import scipy.ndimage
import os
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import cm
from mpl_toolkits.axes_grid import make_axes_locatable
import matplotlib.axes as maxes
#########################################################################################
parser = OptionParser()
parser.add_option("-d", "--fcst_date", dest="date", default="20070501",
help="fcst date", metavar="YYYYMMDD")
(options, args) = parser.parse_args()
file = 'stage2.' + options.date + '.grb'
print file
#open file
gr = Nio.open_file(file,'r')
#read in variables
obs = gr.variables["A_PCP_240_SFC_acc24h"][:]
lats = gr.variables["gridlat_240"][:]
lons = gr.variables["gridlon_240"][:]
# plot precip
fig = plt.figure()
fig.set_size_inches(8,5)
ax = fig.add_axes([0.1,0.1,0.8,0.8])
m = Basemap(projection='lcc',lon_0=-95.,
llcrnrlat=20.,urcrnrlat=50., llcrnrlon=-120.,urcrnrlon=-60.,
lat_1=25.,lat_2=46.,resolution='i',area_thresh=1000.,ax=ax)
m.drawcoastlines(linewidth=0.5)
m.drawstates()
m.drawcountries()
parallels = np.arange(0.,90,5.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
meridians = np.arange(180.,360.,5.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
x, y = m(lons,lats)
intervals = [0.1,2,5,10,15,20,25,35,50,75,100,125,150,200,300,400,500]
obsobj = m.contourf(x,y,obs,intervals,cmap=cm.s3pcpn)
colbar = m.colorbar(obsobj,location="right",size="5%",pad=0.05,extend='both',ticks=intervals)
colbar.set_label('mm',rotation='horizontal',fontsize='medium')
titlestr = 'NCEP STAGE2 Accum precip '
plt.title(titlestr, fontsize='medium')
fname = 'stage2.'+options.date+'.png'
plt.savefig(fname)
#close file
gr.close()
# now do the wrf precip
file = 'wrfnssl.' + options.date + '.grb'
print file
#open file
gr = Nio.open_file(file,'r')
#read in variables
fcst = gr.variables["A_PCP_240_SFC_acc24h"][:]
lats = gr.variables["gridlat_240"][:]
lons = gr.variables["gridlon_240"][:]
# plot precip
fig = plt.figure()
fig.set_size_inches(8,5)
ax = fig.add_axes([0.1,0.1,0.8,0.8])
m = Basemap(projection='lcc',lon_0=-95.,
llcrnrlat=20.,urcrnrlat=50., llcrnrlon=-120.,urcrnrlon=-60.,
lat_1=25.,lat_2=46.,resolution='i',area_thresh=1000.,ax=ax)
m.drawcoastlines(linewidth=0.5)
m.drawstates()
m.drawcountries()
parallels = np.arange(0.,90,5.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
meridians = np.arange(180.,360.,5.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
x, y = m(lons,lats)
intervals = [0.1,2,5,10,15,20,25,35,50,75,100,125,150,200,300,400,500]
fcstobj = m.contourf(x,y,fcst,intervals,cmap=cm.s3pcpn)
colbar = m.colorbar(fcstobj,location="right",size="5%",pad=0.05,extend='both',ticks=intervals)
colbar.set_label('mm',rotation='horizontal',fontsize='medium')
titlestr = 'WRFNSSL Accum precip '
plt.title(titlestr, fontsize='medium')
fname = 'wrfnssl.'+options.date+'.png'
plt.savefig(fname)
#close file
gr.close()
# plot difference
fig = plt.figure()
fig.set_size_inches(8,5)
ax = fig.add_axes([0.1,0.1,0.8,0.8])
m = Basemap(projection='lcc',lon_0=-95.,
llcrnrlat=20.,urcrnrlat=50., llcrnrlon=-120.,urcrnrlon=-60.,
lat_1=25.,lat_2=46.,resolution='i',area_thresh=1000.,ax=ax)
m.drawcoastlines(linewidth=0.5)
m.drawstates()
m.drawcountries()
parallels = np.arange(0.,90,5.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
meridians = np.arange(180.,360.,5.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
intervals = np.arange(-300.,300.,30.)
diffobj = m.contourf(x,y,fcst-obs,intervals,cmap=plt.cm.spectral)
colbar = m.colorbar(diffobj,location="right",size="5%",pad=0.05,extend='both',ticks=intervals)
colbar.set_label('mm',rotation='horizontal',fontsize='medium')
titlestr = 'WRFNSSL-STAGE2 precip '
plt.title(titlestr, fontsize='medium')
fname = 'diff.'+options.date+'.png'
plt.savefig(fname)
# plot a,b,c,d
thresh=1.
a=np.logical_and(fcst>thresh,obs>thresh)
b=np.logical_and(fcst>thresh,obs<=thresh)
c=np.logical_and(fcst<=thresh,obs>thresh)
d=np.logical_and(fcst<=thresh,obs<=thresh)
fig = plt.figure()
fig.set_size_inches(8,5)
ax = fig.add_axes([0.1,0.1,0.8,0.8])
m = Basemap(projection='lcc',lon_0=-95.,
llcrnrlat=20.,urcrnrlat=50., llcrnrlon=-120.,urcrnrlon=-60.,
lat_1=25.,lat_2=46.,resolution='i',area_thresh=1000.,ax=ax)
m.drawcoastlines(linewidth=0.5)
m.drawstates()
m.drawcountries()
parallels = np.arange(0.,90,5.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
meridians = np.arange(180.,360.,5.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
intervals = np.arange(-1.,2.,1.)
aobj = m.contourf(x,y,a,intervals,cmap=plt.cm.spectral)
colbar = m.colorbar(aobj,location="right",size="5%",pad=0.05,extend='both',ticks=intervals)
titlestr = 'WRFNSSL-STAGE2 a '
plt.title(titlestr, fontsize='medium')
fname = 'a.'+options.date+'.png'
plt.savefig(fname)
quit()