Thanks Jeff and Ben.
Below is a complete version of one such program I use. I was give the program by a colleague and I do not know who wrote that code. My previous post showed only as far as the missing data statement.
I'm relatively new to this software. It looks like the missing value goes into the maskdat array and maskdat is an input to pcolor function. Perhaps the best way is to get it working (shading a map) with just the array from the netCDF file read. Essentailly the missing value is to shade some grids outside the region white.
···
#
# Uses basemap's pcolor function. Pcolor accepts arrays
# of the longitude and latitude points of the vertices on the pixels,
# as well as an array with the numerical value in the pixel.
verbose=0 #verbose=2 says a bit more
import sys,getopt
from mpl_toolkits.basemap import Basemap, shiftgrid, cm
#from netCDF3 import Dataset as NetCDFFile
from mpl_toolkits.basemap import NetCDFFile
from pylab import *
alloptions, otherargs= getopt.getopt(sys.argv[1:],'ro:p:X:Y:v:t:l:u:n:') # note
the : after o and p
proj='lam'
#plotfile=None
#plotfile='testmap2.png'
usejetrev=False
colorbounds=[None,None]
extratext=""
xvar=None
yvar=None
thevar=None
thetitle='Bias - Winter Air Temperature'
#thetitle='Bias - Spring Air Temperature'
#thetitle='Bias - Summer Air Temperature'
#thetitle='Bias - Autumn Air Temperature'
therec=None
thelev=None
cbot=None
ctop=None
startlon=-180 #default assumption for starting longitude
for theopt,thearg in alloptions:
print theopt,thearg
if theopt=='-o': # -o needs filename after it, which is now thearg
plotfile=thearg
elif theopt=='-p':
proj=thearg
elif theopt=='-X':
xvar=thearg
elif theopt=='-Y':
yvar=thearg
elif theopt=='-v':
thevar=thearg
elif theopt=='-t':
thetitle=thearg
elif theopt=='-l':
cbot=thearg
elif theopt=='-u':
ctop=thearg
elif theopt=='-n':
therec=thearg
elif theopt=='-m':
thelev=thearg
elif theopt=='-r':
usejetrev=True
else: #something went wrong
print "hmm, what are these??? ", theopt, thearg
sys.exit()
print otherargs
try:
ncname=otherargs[0]
ncfile = NetCDFFile(ncname, 'r')
except:
# ncname = raw_input("\nenter NetCDF file name =>")
# ncfile = NetCDFFile(ncname, 'r')
ncfile = NetCDFFile('simple_xy.nc', 'r') # Here's filename
if verbose>0: #examine the NetCDF file
print "GLOBAL ATTRIBUTES:"
allattr=dir(ncfile)
normalattr=dir(NetCDFFile('/tmp/throwaway.nc','w'))
for item in allattr:
if item not in normalattr: print item,': ',getattr(ncfile,item)
# for name in ncfile.ncattrs():
# print name,'=', getattr(ncfile,name)
# if name=='Conventions' and getattr(ncfile,name)=='COARDS
':
# nmcrean=True
# startlon=0.
# print "guessing this is an NMC reananalysis file
"
# else:
# nmcrean=False
print "\n\nDIMENSIONS:"
for name in ncfile.dimensions.keys():
print "\nDIMENSION:",
try:
print name, ncfile.variables[name][:]
except:
print name, ncfile.dimensions[name]
try:
print ' ',ncfile.variables[name].units
except:
pass
try:
print ' ',ncfile.variables[name].gridtype
except:
pass
print "\n'*****************'\nVARIABLES:"
for name in ncfile.variables.keys():
if name in ncfile.dimensions.keys(): continue
print "\nVARIABLE:",
print name, ncfile.variables[name].dimensions,
try:
print ' ',ncfile.variables[name].units
except:
pass
try:
print " missing value=",ncfile.variables[name].missing
_value
except:
pass
#
print "********************\n"
#if not xvar: xvar=raw_input("\nenter X variable=>")
xvar='rlon'
print "shape of "+xvar+': ',ncfile.variables[xvar].shape
if ncfile.variables[xvar][:].ndim ==1:
print "X is independent of Y"
lon1d=ncfile.variables[xvar][:]
else:
lon1d=False
if ncfile.variables[xvar][:].ndim ==2: lon2d=ncfile.variables[xvar][:]
if ncfile.variables[xvar][:].ndim ==3: lon2d=ncfile.variables[xvar][0,:
,:] #WRF
print "shape of lond2d:", lon2d.shape
#print type(lon1d),lon1d,type(lon1d[1]),lon1d[1]
#
#if not yvar: yvar=raw_input("\nenter Y variable=>")
yvar='rlat'
print "shape of "+yvar+': ',ncfile.variables[yvar].shape
if ncfile.variables[yvar][:].ndim ==1:
print "Y is independent of X"
lat1d=ncfile.variables[yvar][:]
else:
lat1d=False
if ncfile.variables[yvar][:].ndim ==2: lat2d=ncfile.variables[yvar][:]
if ncfile.variables[yvar][:].ndim ==3: lat2d=ncfile.variables[yvar][0,:
,:] #WRF
print "shape of lat2d:", lat2d.shape
#
#if not thevar: thevar=raw_input("\nenter variable to plot=>")
#thevar='nseasondays'
thevar='diff'
#thevar='mean'
data=ncfile.variables[thevar]
dataa=data[:]
theshape=dataa.shape
try:
add_offset=ncfile.variables[thevar].add_offset
print "will use add_offset=",add_offset
except:
add_offset=0.
print "shape of "+thevar+': ',theshape
if len(theshape)>2:
print "there are apparently",theshape[0],"records"
if not therec: therec=raw_input("enter record number to plot=>")
therec=int(therec)
extratext=" rec=%d" %therec
if len(theshape)>3:
print "there are apparently",theshape[1],"levels"
if not thelev: thelev=raw_input("enter level number to plot=>")
thelev=int(thelev)
extratext=extratext+" lev=%d" %thelev
if len(theshape)>3:
slice=dataa[therec,thelev,:,:]+add_offset
elif len(theshape)>2:
slice=dataa[therec,:,:]+add_offset
else:
slice=dataa+add_offset
#data.missing_value=0
#data.missing_value=-999
data.missing_value=-9.99
try:
mval=data.missing_value
print " using specified missing value =",mval
except:
mval=1.e30
print " faking missing value=",mval
# make a masked array: using missing value(s) in mval
maskdat=ma.masked_values(slice,mval)
datamax=max(maskdat.compressed().flat)
datamin=min(maskdat.compressed().flat)
print "\n"+thevar+" has maximum ", datamax," and minimum:", datamin
#colorbounds=[None,None]
#if not cbot: cbot=raw_input("enter lower bound on color scale =>")
#if cbot: colorbounds[0]=float(cbot)
#if not ctop: ctop=raw_input("enter upper bound on color scale =>")
#if ctop: colorbounds[1]=float(ctop)
colorbounds=[-3,3] # for diff in C
print "using clim=",colorbounds
# shift lons and lats by 1/2 grid increment (so values represent the vertices
# of the grid box surrounding the data value, as pcolor expects).
if type(lon1d)!=type(False): #lon1d does exist
delon = lon1d[1]-lon1d[0]
delat = lat1d[1]-lat1d[0]
lons = zeros(len(lon1d)+1,'d')
lats = zeros(len(lat1d)+1,'d')
lons[0:len(lon1d)] = lon1d-0.5*delon
lons[-1] = lon1d[-1]+0.5*delon
lats[0:len(lat1d)] = lat1d-0.5*delon
lats[-1] = lat1d[-1]+0.5*delon
lons, lats = meshgrid(lons, lats)
else:
xdim,ydim=lon2d.shape
lons=zeros([xdim+1,ydim+1],'d')
lats=zeros([xdim+1,ydim+1],'d')
for i in range(1,xdim):
for j in range(1,ydim):
lats[i,j]=.5*lat2d[i,j]+.5*lat2d[i-1,j-1]
lons[i,j]=.5*lon2d[i,j]+.5*lon2d[i-1,j-1]
for i in range(1,xdim):
lons[i,-1]=-lons[i,-3]+2*lons[i,-2]
lats[i,-1]=-lats[i,-3]+2*lats[i,-2]
lons[i,0]=-lons[i,2]+2*lons[i,1]
lats[i,0]=-lats[i,2]+2*lats[i,1]
for j in range(ydim+1):
lons[-1,j]=-lons[-3,j]+2*lons[-2,j]
lats[-1,j]=-lats[-3,j]+2*lats[-2,j]
lons[0,j]=2*lons[1,j]-lons[2,j]
lats[0,j]=2*lats[1,j]-lats[2,j]
#
# alter a matplotlib color table,
# cm.jet is very useful scheme, but reversed colors are better for drought
cmap = cm.get_cmap('jet', 10) # 10 discrete colors
#cmap = cm.get_cmap('summer_r', 10) # 10 discrete colors
print "\nPlotting, please wait...maybe more than 10 seconds"
if proj=='lam': #Lambert Conformal
# m = Basemap(llcrnrlon=-118.0,llcrnrlat=24.,urcrnrlon=-60.0,urcrnrlat=48.
0,\
# resolution='l',area_thresh=1000.,projection='lcc',\
# lat_1=70.,lon_0=-98.)
m = Basemap(llcrnrlon=-80.6,llcrnrlat=38.4,urcrnrlon=-66.0,urcrnrlat=47.
7,\
resolution='l',area_thresh=1000.,projection='lcc',\
lat_1=65.,lon_0=-73.3)
xtxt=200000. #offset for text
ytxt=200000.
parallels = arange(38.,48.,2.)
meridians = arange(-80.,-66.,2.)
else: #cylindrical is default
# m = Basemap(llcrnrlon=-180.,llcrnrlat=-90,urcrnrlon=180.,urcrnrlat=90.,\
# resolution='c',area_thresh=10000.,projection='cyl')
m = Basemap(llcrnrlon=startlon,llcrnrlat=-90,urcrnrlon=startlon+360.,urc
rnrlat=90.,\
resolution='c',area_thresh=10000.,projection='cyl')
xtxt=1.
ytxt=0.
parallels = arange(-90.,90.,30.)
if startlon==-180:
meridians = arange(-180.,180.,60.)
else:
meridians = arange(0.,360.,60.)
if verbose>1: print m.__doc__
plt.rcParams['font.size'] = 20
xsize = rcParams['figure.figsize'][0]
fig=figure(figsize=(xsize,m.aspect*xsize))
#ax = fig.add_axes([0.08,0.1,0.7,0.7],axisbg='white')
ax = fig.add_axes([0.07,0.00,0.8,1.0],axisbg='white')
# make a pcolor plot.
x, y = m(lons, lats)
p = m.pcolor(x,y,maskdat,shading='flat',cmap=cmap)
clim(*colorbounds)
# axes units units are left, bottom, width, height
#cax = axes([0.85, 0.1, 0.05, 0.7]) # colorbar axes for map w/ no graticule
cax = axes([0.88, 0.1, 0.06, 0.81]) # colorbar axes for map w/ graticule
#if datamax>1000.:
# colorbar(format='%7.1e', cax=cax) # draw colorbar
#elif datamax>5.:
# colorbar(format='%d', cax=cax) # draw colorbar
#else:
# colorbar(format='%5.2f', cax=cax) # draw colorbar
colorbar(format='%3.1f', cax=cax)
#text(0.7,1.02,'mm')
text(0.3,1.02,'$^{o}$C',fontsize=20)
axes(ax) # make the original axes current again
# plot blue dot on Norman OK and label it as such.
#if startlon==-180:
# xpt,ypt = m(-97.4836,35.2556)
#else:
# xpt,ypt = m(-97.4836+360.,35.2556)
#m.plot([xpt],[ypt],'bo')
#text(xpt+xtxt,ypt+ytxt,'Norman')
xpt,ypt = m(-70.7500,41.7500)
text(xpt,ypt,'N')
# draw coastlines and political boundaries.
m.drawcoastlines(linewidth=2.0)
m.drawcountries(linewidth=2.0)
m.drawstates(linewidth=2.0)
# draw parallels and meridians.
# label on left, right and bottom of map.
#m.drawparallels(parallels,labels=[1,0,0,0])
#m.drawmeridians(meridians,labels=[1,1,0,1])
if not thetitle:
title(thevar+extratext,fontsize=20)
else:
title(thetitle,fontsize=20)
#if plotfile:
# savefig(plotfile, dpi=72, facecolor='w', bbox_inches='tight', edgecolor=
'w', orientation='portrait')
#else:
# show()
plt.savefig('map.png')
plt.savefig('map.eps')
show()