>________________________________ > From: Jeff Whitaker <jeffrey.s.whita...@noaa.gov> >To: Benjamin Root <ben.r...@ou.edu>; Michael Rawlins <rawl...@geo.umass.edu>; >Matplotlib Users <matplotlib-users@lists.sourceforge.net> >Sent: Wednesday, September 26, 2012 5:10 PM >Subject: Re: [Matplotlib-users] error reading netcdf file > > >Michael: You can't change an attribute in a read-only dataset. It should >never have worked. > >-Jeff > >
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. Mike #!/usr/bin/env python # v0.5 19 June 2010 # General purpose plotter of 2-D gridded data from NetCDF files, # plotted with map boundaries. # NetCDF file should have data in either 2-D, 3-D or 4-D arrays. # Works with the netCDF files in the tutorial, and also the # files available for download at: # http://www.cdc.noaa.gov/cdc/data.ncep.reanalysis.html # Adapted from the basemap example plotmap_pcolor.py, # Some of the variable names from that example are retained. # # 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() ------------------------------------------------------------------------------ How fast is your code? 3 out of 4 devs don\\\'t know how their code performs in production. Find out how slow your code is with AppDynamics Lite. http://ad.doubleclick.net/clk;262219672;13503038;z? http://info.appdynamics.com/FreeJavaPerformanceDownload.html _______________________________________________ Matplotlib-users mailing list Matplotlib-users@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/matplotlib-users
