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
------------------------------------------------------------------------------
The ultimate all-in-one performance toolkit: Intel(R) Parallel Studio XE:
Pinpoint memory and threading errors before they happen.
Find and fix more than 250 security defects in the development cycle.
Locate bottlenecks in serial and parallel code that limit performance.
http://p.sf.net/sfu/intel-dev2devfeb
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
