Re: FillValue error in PyNio

From: Gary Bates <gary.bates_at_nyahnyahspammersnyahnyah>
Date: Mon Nov 30 2009 - 15:25:34 MST

hi Dave,
Thanks for your help.

I've ftp'd the following 2 input files that you'll need:
apcp.2007.nc
lsmask.nc

Also, I've attached a slightly modified version of my python code.
-gary

David Brown wrote:
> 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 15:25:50 2009

This archive was generated by hypermail 2.1.8 : Tue Dec 01 2009 - 09:32:36 MST