m.pcolor: shade a grid based on value

I'm trying to shade grids below a certain threshold a certain color. Using m.pcolor. Also would like to know more about how the routine I've been given, and am using using, (below) works. I'm fairly new to matplotlib so thanks for the help.

Looks like array is filled in these lines:

thevar='mean'
data=ncfile.variables[thevar]
dataa=data[:]
theshape=dataa.shape

Missing values I've chosen set here:

data.missing_value=0
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)

Map shading here, I believe:
p = m.pcolor(x,y,maskdat,shading='flat',cmap=mycolormap)

Questions. 1) How does the function know of the data array? Then, how would I shade grids, with values that are less than zero, purple?

Thanks!
Mike

···

-------------------------------------------------

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='Mean Snow Mass, January-March 1980s, HADCM3_HRM3'
#thetitle='Mean Snow Depth, January-March 1980s, CGCM3_CRCM'
#thetitle='Mean Snow Depth, January-March 2060s, CGCM3_CRCM'
#thetitle='Mean Snow Depth, January-March 1980s, CCSM_WRFG'
#thetitle='Mean Snow Depth, January-March 2060s, CCSM_WRFG'
#thetitle='Mean Snow Depth, January-March 1980s, HADCM3_HRM3'
#thetitle='Mean Snow Depth, January-March 2060s, HADCM3_HRM3'

#thetitle='Mean Snow Depth, January-March 1980s, MultiModel Mean'
thetitle='Mean Snow Depth, January-March 2060s, MultiModel Mean'

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
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=[-10,90] # for diff in days
colorbounds=[0,100]
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
colordict=cm.jet._segmentdata.copy() # dictionary ('blue', 'green', 'red') of nested tuples
for k in colordict.keys():
  colordict[k]=[list(q) for q in colordict[k]] #convert nested tuples to nested list
  for a in colordict[k]:
    a[0]=1.-a[0] #in inner list, change normalized value to 1 - value.
  colordict[k].reverse() #reverse order of outer list
jetrev = cm.colors.LinearSegmentedColormap("jetrev", colordict)
# jetrev can now be used in place of cm.jet
#
# setup of basemap ('lcc' = lambert conformal conic).
# use major and minor sphere radii from WGS84 ellipsoid.
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.)
  xtxt=200000. #offset for text
  ytxt=200000.
  parallels = arange(20.,50.,10.)
# meridians = arange(10.,360.,10.)
  meridians = arange(-130.,-50.,20.)
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.,urcrnrlat=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__
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')
# make a pcolor plot.
x, y = m(lons, lats)
mycolormap=cm.jet
if usejetrev:
  mycolormap=jetrev
  print "using reverse of cm.jet"
p = m.pcolor(x,y,maskdat,shading='flat',cmap=mycolormap)
clim(*colorbounds)
cax = axes([0.85, 0.1, 0.05, 0.7]) # setup colorbar axes
#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='%d', cax=cax)
text(0.7,1.02,'mm')

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')

# draw coastlines and political boundaries.
m.drawcoastlines()
m.drawcountries()
m.drawstates()
# draw parallels and meridians.
# label on left, right and bottom of map.
m.drawparallels(parallels,labels=[1,1,0,0])
m.drawmeridians(meridians,labels=[1,1,0,1])
if not thetitle:
  title(thevar+extratext)
else:
  title(thetitle)
  
#if plotfile:
# savefig(plotfile, dpi=72, facecolor='w', bbox_inches='tight', edgecolor='w', orientation='portrait')
#else:
# show()

plt.savefig('map.png')
# comment show to mass produce

____________________________________________________________________________________
Now that's room service! Choose from over 150,000 hotels
in 45,000 destinations on Yahoo! Travel to find your fit.
http://farechase.yahoo.com/promo-generic-14795097