Re: FillValue error in PyNio

Hi Gary,

Can you provide your input NetCDF file or let me know where I could
get a similar file that would demonstrate the problem?
To send us the file do this:
ftp ftp.cgd.ucar.edu
login as anonymous
password your email address
cd incoming
bin
put <the_file>

Let me know the name of the file.
Thanks,
  -dave

On Nov 30, 2009, at 11:37 AM, Gary Bates wrote:

> hi,
> I have a python code (attached) which reads 1 year of daily data
> from a netcdf file and writes out a subset of this data. The code
> successfully reads and writes out the 1st day's data (1 Jan), but I
> get the following error when trying to write data for 2 Jan:
>
> Traceback (most recent call last):
> File "test.py", line 154, in <module>
> apcpdat[nr] = precip.astype('f')
> File "/opt/local/Library/Frameworks/Python.framework/Versions/2.5/
> lib/python2.5/site-packages/PyNIO/Nio.py", line 316, in __setitem__
> fval = self.__dict__['_FillValue'][0]
> IndexError: 0-d arrays can't be indexed
>
> It appears that somehow _FillValue was an array for 1 Jan but
> becomes a scalar for 2 Jan. Why? I never set _FillValue in the code.
> -Gary
>
>
> p.s. On the line before the error, print precip.__dict__ gives:
>
> {'_baseclass': <type 'numpy.ndarray'>, '_hardmask': False,
> '_optinfo': {}, '_basedict': {}, '_sharedmask': False, '_isfield':
> False, '_mask': False, '_fill_value': None}
> import os, sys, math
> import dateutils
> import numpy
> import Nio
> import proj
>
> def round(input):
> "take a float and return the nearest integer"
> remainder = math.ceil(input) - input
> if remainder > 0.5:
> result = math.floor(input)
> else:
> result = math.ceil(input)
> return int(result)
>
> field = 'apcp'
>
> # open file containing NARR data, get lons and lats.
> filename = '/Datasets/NARR/monolevel/apcp.2001.nc'
> ncfile = Nio.open_file(filename,'r')
> latsout = numpy.array(ncfile.variables['lat'][:,:])
> lonsout = numpy.array(ncfile.variables['lon'][:,:])
> xgrid = numpy.array(ncfile.variables['x'][:])
> ygrid = numpy.array(ncfile.variables['y'][:])
>
> nx = len(xgrid)
> dx = xgrid[-1]/(nx-1)
> ny = len(ygrid)
> dy = ygrid[-1]/(ny-1)
>
> # make lonsout go from 0 to 360
> lonsout = numpy.where(lonsout < 0, lonsout+360, lonsout)
> latmin = ncfile.latcorners[0]
> latmax = ncfile.latcorners[2]
> lonmin = ncfile.loncorners[0]
> lonmax = ncfile.loncorners[2]
> print latmin, latmax, lonmin, lonmax
>
> # define awips grid 221 projection.
> projparams = {}
> projparams['proj'] = 'lcc'
> projparams['R'] = 6371200
> projparams['lat_1'] = ncfile.standardpar1[0]
> projparams['lat_2'] = ncfile.standardpar1[0]
> projparams['lon_0'] = ncfile.centerlon[0]
> awips221 =
> proj
> .Proj
> (projparams
> ,llcrnrlon
> =
> lonmin
> ,llcrnrlat=latmin,urcrnrx=lonmax,urcrnry=latmax,urcrnrislatlon=True)
> ncfile.close()
>
> # read in cpc .25 degree us mask.
> file="/Datasets/cpc_us_precip/lsmask.nc"
> ncfile = Nio.open_file(file,'r')
> lats = numpy.array(ncfile.variables['lat'][:])
> lons = numpy.array(ncfile.variables['lon'][:]) # lons go from 0 to
> 360.
> mask = numpy.array(ncfile.variables['lsmask'][0,:,:])
> mask[0:16,200:320]=-1 # get rid of cuba.
> # interpolate mask to eta 221 grid.
> maskout =
> proj
> .interp
> (mask
> .astype
> ('f'),lons,lats,lonsout,latsout,order=1,mode='constant',cval=-1.0)
> ncfile.close()
>
> # define boundaries of us subset region.
> latminout = 22.0
> lonminout = -125.
> xmin,ymin = awips221(lonminout,latminout)
> xmax = xmin + 21*8*dx
> ymax = ymin + 15*8*dy
> lonmaxout,latmaxout = awips221(xmax,ymax,inverse=True)
> print lonminout,latminout,lonmaxout,latmaxout
>
> # define indices of slice for subset region.
> ixmin = round((len(xgrid)-1)*(xmin-xgrid[0])/(xgrid[-1]-xgrid[0]))
> iymin = round((len(ygrid)-1)*(ymin-ygrid[0])/(ygrid[-1]-ygrid[0]))
> ixmax = round((len(xgrid)-1)*(xmax-xgrid[0])/(xgrid[-1]-xgrid[0]))
> iymax = round((len(ygrid)-1)*(ymax-ygrid[0])/(ygrid[-1]-ygrid[0]))
> print 'ixmin,ixmax,iymin,iymax:', ixmin, ixmax, iymin, iymax
>
> # conus where usmask=1, otherwise usmask=0
> usmask = (numpy.where(maskout[iymin:iymax+1,ixmin:ixmax+1] < 0, 0,
> 1)).astype('h')
>
> x = xgrid[ixmin:ixmax+1]-xmin
> y = ygrid[iymin:iymax+1]-ymin
> longitudes = lonsout[iymin:iymax+1,ixmin:ixmax+1]
> latitudes = latsout[iymin:iymax+1,ixmin:ixmax+1]
>
>
> year = '1979'
> filename = '/nas/gbates/narr/'+field+'usdaily_'+year+'.nc'
> try:
> os.remove(filename)
> except:
> print 'Not found: ', filename
> ncfile = Nio.open_file(filename,'c')
>
> ncfile.create_dimension('time', None)
> ncfile.create_dimension('ydist', iymax-iymin+1)
> ncfile.create_dimension('xdist', ixmax-ixmin+1)
>
> xdat = ncfile.create_variable('xdist', 'f', ('xdist',))
> ydat = ncfile.create_variable('ydist', 'f', ('ydist',))
> latsdat = ncfile.create_variable('lat', 'f', ('ydist','xdist',))
> lonsdat = ncfile.create_variable('lon', 'f', ('ydist','xdist',))
> maskdat = ncfile.create_variable('conusmask', 'h', ('ydist','xdist',))
> timesdat = ncfile.create_variable('time', 'l', ('time',))
> timeunits='hours since 1800-1-1 00:00:0.0'
>
> xdat[:] = x.astype('f')
> ydat[:] = y.astype('f')
> latsdat[:] = latitudes.astype('f')
> lonsdat[:] = longitudes.astype('f')
> maskdat[:,:] = usmask.astype('h')
>
> setattr(latsdat, 'long_name','latitude')
> setattr(latsdat, 'units', 'degrees_north')
> setattr(lonsdat, 'units', 'degrees_east')
> setattr(lonsdat, 'long_name','longitude')
> setattr(xdat, 'long_name','eastward distance from southwest corner
> of domain in projection coordinates')
> setattr(ydat, 'long_name','northward distance from southwest corner
> of domain in projection coordinates')
> setattr(maskdat, 'long_name','continental united states mask')
> setattr(xdat,'units','m')
> setattr(ydat,'units','m')
> setattr(timesdat, 'units', timeunits)
> setattr(timesdat, 'long_name','analysis time')
> setattr(ncfile, 'standardpar1', 50.)
> setattr(ncfile, 'standardpar2', 50.)
> setattr(ncfile, 'centerlon', -107.)
> setattr(ncfile, 'centerlat', 50.)
> setattr(ncfile, 'title','NCEP North American Regional Reanalysis
> Data')
>
> apcpdat = ncfile.create_variable(field, 'f',
> ('time','ydist','xdist',))
> setattr(apcpdat, 'units', 'mm')
> setattr(apcpdat, 'long_name', 'Daily Accumulated Precipitation')
> setattr(apcpdat, 'least_significant_digit', 1)
>
> filename = '/Datasets/NARR/monolevel/'+field+'.'+year+'.nc'
> narrfile = Nio.open_file(filename,'r')
> times = narrfile.variables['time']
> print times
> iprecip = narrfile.variables[field]
> print iprecip, iprecip.shape, iprecip.scale_factor, iprecip.add_offset
>
> startdate = year + '010100'
> enddate = year + '123100'
> print startdate, enddate
> nr=0
> for nt,date in enumerate(dateutils.daterange(startdate,enddate,24)):
> nt1 = nt*8
> precip8 = iprecip.scale_factor*iprecip[nt1:nt1+8,iymin:iymax
> +1,ixmin:ixmax+1] + iprecip.add_offset
> precip = numpy.sum(precip8,axis=0)
> print date,nt1,times[nt1],numpy.amin(precip),numpy.amax(precip)
>
> if field == 'apcp' or field == 'prate': precip =
> numpy.where(precip > 0., precip, 0.)
>
> print apcpdat.shape
> print precip, precip.shape
> print precip.__dict__
> apcpdat[nr] = precip.astype('f')
> timesdat[nr] = times[nt1]
> nr += 1
>
> narrfile.close()
>
> ncfile.close()
> _______________________________________________
> pyngl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/pyngl-talk

_______________________________________________
pyngl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/pyngl-talk
Received on Mon Nov 30 13:00:34 2009