Re: [Matplotlib-users] error reading netcdf file

[Jeff W. <jef...@no...>]

On 9/26/12 1:41 PM, Benjamin Root wrote:
>
>
> On Wed, Sep 26, 2012 at 2:07 PM, Michael Rawlins <raw...@ya... 
> <mailto:raw...@ya...>> wrote:
>
>
>
>         ------------------------------------------------------------------------
>         *From:* Benjamin Root <ben...@ou... <mailto:ben...@ou...>>
>         *To:* Michael Rawlins <raw...@ya...
>         <mailto:raw...@ya...>>
>         *Cc:* "mat...@li...
>         <mailto:mat...@li...>"
>         <mat...@li...
>         <mailto:mat...@li...>>
>         *Sent:* Wednesday, September 26, 2012 10:33 AM
>         *Subject:* Re: [Matplotlib-users] error reading netcdf file
>
>
>
>         On Wed, Sep 26, 2012 at 10:27 AM, Michael Rawlins
>         <raw...@ya... <mailto:raw...@ya...>> wrote:
>
>             Recently built and installed netCDF4-1.0. I'm running a
>             script that has worked on two other linux OS systems. Error:
>
>             File "test.py <http://test.py/>", line 96, in <module>
>                 data.missing_value=-9.99
>               File "netCDF4.pyx", line 2570, in
>             netCDF4.Variable.__setattr__ (netCDF4.c:28242)
>               File "netCDF4.pyx", line 2392, in
>             netCDF4.Variable.setncattr (netCDF4.c:26309)
>               File "netCDF4.pyx", line 1013, in netCDF4._set_att
>             (netCDF4.c:12699)
>             AttributeError: NetCDF: Write to read only
>
>
>             The statement in the code triggers the error is:
>
>             data.missing_value=-9.99 .
>
>             MR
>
>
>
>
>         This typically happens when one opens a netCDF4 Dataset object
>         in 'r' mode, and/or if the file permissions for the file was
>         set to read-only.  When modifying an attribute, it is
>         technically trying to write to the file.
>
>         Ben Root
>
>
>
>         Here is the file open statement:
>
>         ncfile = NetCDFFile('statsPrcp_Uncertainty2_winter.nc', 'r')
>
>         This worked fine in earlier versions.  No where in the code do
>         I try to write to a netCDF file. So I'm confused as to why
>         specifying the read of the netCDF file would result in an
>         error when designating the missing data value.
>
>         Here's the code:
>
>         # 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 netCDF4 import Dataset as NetCDFFile
>         #from mpl_toolkits.basemap import NetCDFFile
>         from pylab import *
>         import  matplotlib.pyplot as plt
>
>         alloptions, otherargs=
>         getopt.getopt(sys.argv[1:],'ro:p:X:Y:v:t:l:u:n:') # note the :
>         after o and p
>         proj='lam'
>
>         cmap = cm.get_cmap('jet_r', 10)    # 10 discrete colors
>
>         print "\nPlotting, please wait...maybe more than 10 seconds"
>         if proj=='lam': #Lambert Conformal
>             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.,-64.,2.)
>
>         xsize = rcParams['figure.figsize'][0]
>         fig=figure(figsize=(xsize,m.aspect*xsize))
>         #cax = axes([0.88, 0.1, 0.06, 0.81])  # colorbar axes for map
>         w/ graticule
>
>         ############################################################################################
>         subplots_adjust(left=None, bottom=None, right=0.9, top=None,
>         wspace=0.05, hspace=0.03)
>         # Make the first map at upper left
>         plt.subplot(2,2,1)
>
>         ncfile = NetCDFFile('statsPrcp_Uncertainty2_winter.nc', 'r')
>
>         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,:,:]
>             print "shape of lond2d:", lon2d.shape
>
>         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,:,:]
>             print "shape of lat2d:", lat2d.shape
>
>         thevar='diff'
>         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=-9.99
>
>
> I think I know why this used to work in the past.  I think netCDF4 
> used to recognize the "missing_value" attribute as a special value.  
> It is now ._FillValue.  But I am still not sure you can save a value 
> to that attribute in a read-only Dataset. You might have to perform 
> the mask after extracting the data.
>
> Note, setting the _FillValue *after* extracting the data to the 
> "dataa" variable won't change "dataa".  I don't see how setting 
> data.missing_value is supposed to actually help you.
>
> Ben Root

Michael:  You can't change an attribute in a read-only dataset.  It 
should never have worked.

-Jeff

-- 
Jeffrey S. Whitaker         Phone  : (303)497-6313
Meteorologist               FAX    : (303)497-6449
NOAA/OAR/PSD  R/PSD1        Email  : Jef...@no...
325 Broadway                Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web    : http://tinyurl.com/5telg