netcdf4-python build

After the build, I determined that ‘install’ was also needed.

python setup.py install

completed with no errors. OK, finally built and installed. But now my matplotlib script gives this error:

Traceback (most recent call last):
File “map_PrcpBias_Northeast.py”, line 21, in
from netCDF4 import Dataset as NetCDFFile
ImportError: /usr/local/lib/python2.7/dist-packages/netCDF4.so: undefined symbol: nc_inq_var_endian

So, checking shared library dependencies:

ldd
/usr/local/lib/python2.7/dist-packages/netCDF4.so

linux-gate.so.1 =>  (0xb776a000)
libnetcdf.so.7 => /usr/local/lib/libnetcdf.so.7 (0xb7604000)
libpthread.so.0 => /lib/i386-linux-gnu/libpthread.so.0 (0xb75d3000)
libc.so.6 => /lib/i386-linux-gnu/libc.so.6 (0xb742d000)
libm.so.6 => /lib/i386-linux-gnu/libm.so.6 (0xb7401000)
/lib/ld-linux.so.2 (0xb776b000)

and

ldd /usr/local/lib/libnetcdf.so.7

linux-gate.so.1 =>  (0xb7765000)
libm.so.6 => /lib/i386-linux-gnu/libm.so.6 (0xb7663000)
libc.so.6 => /lib/i386-linux-gnu/libc.so.6 (0xb74be000)
/lib/ld-linux.so.2 (0xb7766000)

no libhdf5 there. Can this be
fixed?

MR

From what I remember dealing with the netcdf c library, you have to explicitly set a compile flag to enable hdf5 support. That was a while ago, though. I’m not sure if things have changed.

Hope this helps.

···

On Friday, September 21, 2012, Michael Rawlins wrote:

After the build, I determined that ‘install’ was also needed.

python setup.py install

completed with no errors. OK, finally built and installed. But now my matplotlib script gives this error:

Traceback (most recent call last):
File “map_PrcpBias_Northeast.py”, line 21, in
from netCDF4 import Dataset as NetCDFFile
ImportError: /usr/local/lib/python2.7/dist-packages/netCDF4.so: undefined symbol: nc_inq_var_endian

So, checking shared library dependencies:

ldd
/usr/local/lib/python2.7/dist-packages/netCDF4.so

linux-gate.so.1 =>  (0xb776a000)
libnetcdf.so.7 => /usr/local/lib/libnetcdf.so.7 (0xb7604000)
libpthread.so.0 => /lib/i386-linux-gnu/libpthread.so.0 (0xb75d3000)

libc.so.6 => /lib/i386-linux-gnu/libc.so.6 (0xb742d000)
libm.so.6 => /lib/i386-linux-gnu/libm.so.6 (0xb7401000)
/lib/ld-linux.so.2 (0xb776b000)

and

ldd /usr/local/lib/libnetcdf.so.7

linux-gate.so.1 =>  (0xb7765000)
libm.so.6 => /lib/i386-linux-gnu/libm.so.6 (0xb7663000)
libc.so.6 => /lib/i386-linux-gnu/libc.so.6 (0xb74be000)
/lib/ld-linux.so.2 (0xb7766000)

no libhdf5 there. Can this be
fixed?

MR


Damon McDougall
http://www.damon-is-a-geek.com
B2.39
Mathematics Institute
University of Warwick

Coventry
West Midlands
CV4 7AL
United Kingdom


From: Damon McDougall <damon.mcdougall@…287…>
To: Michael Rawlins <rawlins02@…9…>
Cc: “matplotlib-users@…564…net”
matplotlib-users@lists.sourceforge.net
Sent: Friday, September 21, 2012 2:50 PM
Subject: Re: [Matplotlib-users] netcdf4-python build

From what I remember dealing with the netcdf c library, you have to explicitly set a compile flag to enable hdf5 support. That was a while ago, though. I’m not sure if things have changed.

Hope this helps.

I did not do a source compile.

I’ve reinstalled libhdf5-serial-1.8.4 package and now have the right HDF5 library files. Successfully built and installed netCDF4-1.0. My script loads in module but a read/write issue is present. This is code that worked on previous system. Error:

File “test.py”, line 96, in
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.

data.missing_value=-9.99 .

MR

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”, line 96, in
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

···

On Wed, Sep 26, 2012 at 10:27 AM, Michael Rawlins <rawlins02@…9…> 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”, line 96, in

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


From: Benjamin Root <ben.root@…1304…>
To: Michael Rawlins <rawlins02@…9…>
Cc:matplotlib-users@lists.sourceforge.net
matplotlib-users@lists.sourceforge.net
Sent: Wednesday, September 26, 2012 10:33 AM
Subject: Re: [Matplotlib-users] error reading netcdf file

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”, line 96, in

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

···

On Wed, Sep 26, 2012 at 10:27 AM, Michael Rawlins <rawlins02@…9…> wrote:

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

···

On Wed, Sep 26, 2012 at 2:07 PM, Michael Rawlins <rawlins02@…9…> wrote:


From: Benjamin Root <ben.root@…1304…>
To: Michael Rawlins <rawlins02@…9…>
Cc: "matplotlib-users@lists.sourceforge.net "
matplotlib-users@lists.sourceforge.net
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 <rawlins02@…9…> 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”, line 96, in

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

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

···

On 9/26/12 1:41 PM, Benjamin Root
wrote:

    On Wed, Sep 26, 2012 at 2:07 PM, Michael

Rawlins <rawlins02@…9…>
wrote:


From:
Benjamin Root <ben.root@…1304…>
To:
Michael Rawlins <rawlins02@…9… >
Cc:
"matplotlib-users@lists.sourceforge.net "
<matplotlib-users@lists.sourceforge.net >
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 <rawlins02@…9…>
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

                                    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](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
-- Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : 325 Broadway Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web :

Jeffrey.S.Whitaker@…259…http://tinyurl.com/5telg

________________________________
From: Jeff Whitaker <jeffrey.s.whitaker@...259...>
To: Benjamin Root <ben.root@...1304...>; Michael Rawlins <rawlins@...4218.....>; 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()

<snip>

I've gotten my script to work by setting the missing value another way. I'm showing lines commented out as well.

#data.missing_value=-9.99
#try:
# mval=data.missing_value
# print " using specified missing value =",mval
#except:
mval=-9.99
# print " faking missing value=",mval
# make a masked array: using missing value(s) in mval
maskdat=ma.masked_values(slice,mval)

If I understand correctly, in previous versions the value of -9.99 was passed through the data.missing_value into mval. A mistake in coding that did not result in an error. I will implement the change. Tomorrow I'll check that this solution is backward compatible.

Mike

···

----- Original Message -----

From: Michael Rawlins <rawlins02@...9...>
To: Jeff Whitaker <jeffrey.s.whitaker@...259...>; Benjamin Root <ben.root@...1304...>; Matplotlib Users <matplotlib-users@lists.sourceforge.net>
Cc:
Sent: Wednesday, September 26, 2012 6:23 PM
Subject: Re: [Matplotlib-users] error reading netcdf file

________________________________
From: Jeff Whitaker <jeffrey.s.whitaker@...259...>
To: Benjamin Root <ben.root@...1304...>; Michael Rawlins

<rawlins@...3986...>; 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.

Michael,

The .missing_value attribute is not used anymore (It is ._FillValue now). Anyway, if your data had any value that matched ._FillValue, then, by default, netCDF4 will give you a masked array anyway. You will only need to set the mask if the fill value doesn’t exist or if it is different from what you were expecting.

Ben Root

________________________________
From: Benjamin Root <ben.root@...1304...>
To: Michael Rawlins <rawlins02@...9...>
Cc: Matplotlib Users <matplotlib-users@lists.sourceforge.net>
Sent: Thursday, September 27, 2012 2:29 PM
Subject: Re: [Matplotlib-users] error reading netcdf file

Michael,

The .missing_value attribute is not used anymore (It is ._FillValue now). Anyway, if your data had any value that matched ._FillValue, then, by default, netCDF4 will give you a masked array anyway. You will only need to set the mask if the fill value doesn't exist or if it is different from what you were expecting.

Ben Root

Thanks Ben. I'll play around with ._FillValue and masked arrays. Glad to have everything working again with this new OS.

Mike