plot grib data in python

install Basemap on coates


#untar the basemap tarball

cd $HOME
mkdir $HOME/GEOS
mkdir -p $HOME/local/lib/python2.7/site-packages
tar -zxvf /usr/rmt_share/baldwin/basemap-1.0.5.tar.gz
cd basemap-1.0.5

#install GEOS
cd geos-3.3.3
export GEOS_DIR=$HOME/GEOS
./configure --prefix=$GEOS_DIR
make; make install

cd $HOME/basemap-1.0.5
module load python/2.7
module load ncl
python setup.py install --prefix=$HOME/local

# add these lines to your .bash_profile
export NETCDF_HOME=/apps/rhel5/graphics/x86_64/gnu/gcc-4.1.2/netcdf-4.1.1
export PYTHONPATH=${PYTHONPATH}:~/local/lib/python2.7/site-packages
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${NETCDF_HOME}/lib
export PATH=${PATH}:${NETCDF_HOME}/bin

python script to print out variable names in GRIB file

/usr/rmt_share/baldwin/EAS591/gribvars.py

#!/usr/bin/env python

import matplotlib.pyplot as plt
import math
import numpy as np
from pylab import *
import Nio
import scipy.ndimage
from optparse import OptionParser
import os
import datetime
import time
import sys
import textwrap
import subprocess

#########################################################################################


file='stage2.20070501.grb'
#open file(s)
gr = Nio.open_file(file,'r',format='grb')
print gr
gr = Nio.open_file(file,'r',format='grb')
print gr
names =  gr.variables.keys()
print "Variable names:", names

#dimensions
dims = gr.dimensions
print "dimensions:",dims

#global attributes
attr = gr.attributes.keys()
print "global attributes",attr

data = gr.variables['gridlat_240']
attrib = data.attributes.keys()
print attrib
dx = getattr(data,'Dx')
dy = getattr(data,'Dy')
print dx, dy

python script to read in and plot QPF and observed precip

/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()


Comments