python to read/plot grib

log onto conte, if you do not have an account, let me know

ssh -Y conte.rcac.purdue.edu

load a specific python module

module load python/2.7.2

download the grib file of interest

wget ftp://ldas.ncep.noaa.gov/nldas2/nco_nldas/2017/nldas.20170104/noah.t12z.grbf09

to get an inventory of the grib file, run the following script: edit a file called "gribvars.py" and paste the following code in, then do "python gribvars.py"
#!/usr/bin/env python

import Nio

file='noah.t12z.grbf09'

#open file(s)
gr = Nio.open_file(file,'r',format='grb2')
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['lat_0']
attrib = data.attributes.keys()
print attrib
#'GridType', 'Dj', 'Di', 'Lo2', 'long_name', 'Lo1', 'units', 'La1', 'La2'
dj = getattr(data,'Dj')
di = getattr(data,'Di')
print dj, di

to print values of specific variables at the nearest grid point to a given lat/lon
to plot a 2D map of a variable
run the following script: 
edit a file called "gribplot.py" and paste the following code in, then do "python gribplot.py"
grib file name and lat (-y) / lon (-x) values can be specified as options in the call:
python gribplot.py -f filename.grib -y 40.418 -x -86.933

#!/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

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

#constants


dx=0.125


#########################################################################################
parser = OptionParser()
parser.add_option("-f", "--file_name", dest="name", default="noah.t12z.grbf09",
                  help="input file name", metavar="file name")
parser.add_option("-y", "--latitude", dest="lat_location", type="float", default=40.4148,
                  help="north latitude location", metavar="loc")
parser.add_option("-x", "--longitude", dest="lon_location", type="float", default=-86.9333,
                  help="east long location", metavar="loc")


(options, args) = parser.parse_args()

file = options.name
print file
gr = Nio.open_file(file,'r',format='grb2')
#read in model variables
tsfc = gr.variables["AVSFT_P0_L1_GLL0"][:]
snod = gr.variables["SNOD_P0_L1_GLL0"][:]
swrad = gr.variables["NSWRF_P8_L1_GLL0_avg"][:]
lwrad = gr.variables["NLWRF_P8_L1_GLL0_avg"][:]
lhflux = gr.variables["LHTFL_P8_L1_GLL0_avg"][:]
shflux = gr.variables["SHTFL_P8_L1_GLL0_avg"][:]
gflux = gr.variables["GFLUX_P8_L1_GLL0_avg"][:]
snow = gr.variables["TSNOWP_P8_L1_GLL0_acc"][:]
rain = gr.variables["TWATP_P8_L1_GLL0_acc"][:]
lat = gr.variables["lat_0"][:]
lon = gr.variables["lon_0"][:]
lons, lats = np.meshgrid(lon,lat)
lons = lons-360.

#close file
gr.close()

[imax,jmax]=shape(tsfc)

#print values at a point
lon1=options.lon_location
lat1=options.lat_location
deltax=((lon1-lons)*np.cos(lats*np.pi/180))**2
deltay=(lat1-lats)**2
dist=np.sqrt(deltax+deltay)
loc=np.unravel_index(np.argmin(dist),dist.shape)
row1=loc[0]
col1=loc[1]
print 'lat lon requested = ', lat1,lon1
print 'lat lon of nearest grid point = ', lats[row1,col1],lons[row1,col1]
print 'surface temp at nearest grid point = ', tsfc[row1,col1]
print 'sw radiation at nearest grid point = ', swrad[row1,col1]
print 'sens heat flx at nearest grid point = ', shflux[row1,col1]

# set up sub domain for plot
lat1=37.2
lon1=-90.3
lat2=42.4
lon2=-82.2
#lat1=lats[0,0]
#lon1=lons[0,0]
#lat2=lats[imax-1,jmax-1]
#lon2=lons[imax-1,jmax-1]
m = Basemap(projection='cyl',
          llcrnrlat=lats[0,0],urcrnrlat=lats[imax-1,jmax-1], llcrnrlon=lons[0,0],urcrnrlon=lons[imax-1,jmax-1],
          resolution='i',area_thresh=1000.)
xpt, ypt = m(lon1,lat1)
jll=floor(xpt/dx)
ill=floor(ypt/dx)
xpt, ypt = m(lon2,lat2)
jur=floor(xpt/dx)
iur=floor(ypt/dx)
lonll,latll=m(jll*dx,ill*dx,inverse=True)
lonur,latur=m(jur*dx,iur*dx,inverse=True)
print latll,lonll,latur,lonur
print imax,jmax
print ill,jll,iur,jur




fig = plt.figure()
fig.set_size_inches(8,6)
matplotlib.rcParams['contour.negative_linestyle'] = 'solid'
ax = fig.add_axes([0.1,0.1,0.8,0.8])
m = Basemap(projection='lcc',lon_0=-90,
          llcrnrlat=latll,urcrnrlat=latur, llcrnrlon=lonll,urcrnrlon=lonur,
          lat_1=45.,lat_2=28.,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)
apcpobj = m.contourf(x,y,tsfc)
colbar = m.colorbar(apcpobj,location="right",size="5%",pad=0.05,extend='both')
colbar.set_label('K',rotation='horizontal',fontsize='medium')
titlestr = 'surface temp'
plt.title(titlestr, fontsize='medium')
fname = 'tsfc.png'
plt.savefig(fname)
plt.close()


Comments