In [1]:
import sys
sys.path.append('../scripts')
import constants as c
import qcio
import qcutils
import netCDF4
In [2]:
def nc_read_series_3d(ncFullName):
''' Read a netCDF file and put the data and meta-data into a DataStructure'''
print ' Reading netCDF file '+ncFullName
netCDF4.default_encoding = 'latin-1'
ds = qcio.DataStructure()
# check to see if the requested file exists, return empty ds if it doesn't
if not qcutils.file_exists(ncFullName,mode="quiet"):
print ' netCDF file '+ncFullName+' not found'
raise Exception("GapFillFromACCESS: ACCESS file not found")
# file probably exists, so let's read it
ncFile = netCDF4.Dataset(ncFullName,'r')
# now deal with the global attributes
gattrlist = ncFile.ncattrs()
if len(gattrlist)!=0:
for gattr in gattrlist:
ds.globalattributes[gattr] = getattr(ncFile,gattr)
if "time_step" in ds.globalattributes: c.ts = ds.globalattributes["time_step"]
# get a list of the variables in the netCDF file (not their QC flags)
varlist = [x for x in ncFile.variables.keys() if "_QCFlag" not in x]
for ThisOne in varlist:
# skip variables that do not have time as a dimension
dimlist = [x.lower() for x in ncFile.variables[ThisOne].dimensions]
if "time" not in dimlist: continue
# create the series in the data structure
ds.series[unicode(ThisOne)] = {}
# get the data and the QC flag
nDims = len(ncFile.variables[ThisOne].shape)
if nDims==1:
# single dimension
ds.series[ThisOne]["Data"] = ncFile.variables[ThisOne][:]
# netCDF4 returns a masked array if the "missing_variable" attribute has been set
# for the variable, here we trap this and force the array in ds.series to be ndarray
if numpy.ma.isMA(ds.series[ThisOne]["Data"]):
ds.series[ThisOne]["Data"],dummy = qcutils.MAtoSeries(ds.series[ThisOne]["Data"])
# check for a QC flag
if ThisOne+'_QCFlag' in ncFile.variables.keys():
# load it from the netCDF file
ds.series[ThisOne]["Flag"] = ncFile.variables[ThisOne+'_QCFlag'][:]
else:
# create an empty flag series if it does not exist
nRecs = numpy.size(ds.series[ThisOne]["Data"])
ds.series[ThisOne]["Flag"] = numpy.zeros(nRecs,dtype=numpy.int32)
elif nDims==3:
# 3 dimensions
ds.series[ThisOne]["Data"] = ncFile.variables[ThisOne][:,0,0]
# netCDF4 returns a masked array if the "missing_variable" attribute has been set
# for the variable, here we trap this and force the array in ds.series to be ndarray
if numpy.ma.isMA(ds.series[ThisOne]["Data"]):
ds.series[ThisOne]["Data"],dummy = qcutils.MAtoSeries(ds.series[ThisOne]["Data"])
# check for a QC flag
if ThisOne+'_QCFlag' in ncFile.variables.keys():
# load it from the netCDF file
ds.series[ThisOne]["Flag"] = ncFile.variables[ThisOne+'_QCFlag'][:,0,0]
else:
# create an empty flag series if it does not exist
nRecs = numpy.size(ds.series[ThisOne]["Data"])
ds.series[ThisOne]["Flag"] = numpy.zeros(nRecs,dtype=numpy.int32)
# get the variable attributes
vattrlist = ncFile.variables[ThisOne].ncattrs()
ds.series[ThisOne]["Attr"] = {}
if len(vattrlist)!=0:
for vattr in vattrlist:
ds.series[ThisOne]["Attr"][vattr] = getattr(ncFile.variables[ThisOne],vattr)
ncFile.close()
# make sure all values of -9999 have non-zero QC flag
#qcutils.CheckQCFlags(ds)
# get a series of Python datetime objects
#qcutils.get_datetimefromymdhms(ds)
# get series of UTC datetime
#qcutils.get_UTCfromlocaltime(ds)
return ds
In [3]:
def nc_read_series_1d(ncFullName):
''' Read a netCDF file and put the data and meta-data into a DataStructure'''
print ' Reading netCDF file '+ncFullName
netCDF4.default_encoding = 'latin-1'
ds = qcio.DataStructure()
# check to see if the requested file exists, return empty ds if it doesn't
if not qcutils.file_exists(ncFullName,mode="quiet"):
print ' netCDF file '+ncFullName+' not found'
raise Exception("GapFillFromACCESS: ACCESS file not found")
# file probably exists, so let's read it
ncFile = netCDF4.Dataset(ncFullName,'r')
# now deal with the global attributes
gattrlist = ncFile.ncattrs()
if len(gattrlist)!=0:
for gattr in gattrlist:
ds.globalattributes[gattr] = getattr(ncFile,gattr)
if 'time_step' in ds.globalattributes: c.ts = ds.globalattributes['time_step']
for ThisOne in ncFile.variables.keys():
if '_QCFlag' not in ThisOne:
# create the series in the data structure
ds.series[unicode(ThisOne)] = {}
# get the data variable object
ds.series[ThisOne]['Data'] = ncFile.variables[ThisOne][:]
# netCDF4 returns a masked array if the "missing_variable" attribute has been set
# for the variable, here we trap this and force the array in ds.series to be ndarray
if numpy.ma.isMA(ds.series[ThisOne]["Data"]):
ds.series[ThisOne]["Data"],dummy = qcutils.MAtoSeries(ds.series[ThisOne]["Data"])
# check for a QC flag and if it exists, load it
if ThisOne+'_QCFlag' in ncFile.variables.keys():
ds.series[ThisOne]['Flag'] = ncFile.variables[ThisOne+'_QCFlag'][:]
else:
nRecs = numpy.size(ds.series[ThisOne]['Data'])
ds.series[ThisOne]['Flag'] = numpy.zeros(nRecs,dtype=numpy.int32)
# get the variable attributes
vattrlist = ncFile.variables[ThisOne].ncattrs()
if len(vattrlist)!=0:
ds.series[ThisOne]['Attr'] = {}
for vattr in vattrlist:
ds.series[ThisOne]['Attr'][vattr] = getattr(ncFile.variables[ThisOne],vattr)
ncFile.close()
# get a series of Python datetime objects
qcutils.get_datetimefromymdhms(ds)
# get series of UTC datetime
qcutils.get_UTCfromlocaltime(ds)
return ds
In [37]:
def CheckQCFlags(ds):
for ThisOne in ds.series.keys():
data = numpy.ma.masked_values(ds.series[ThisOne]["Data"],-9999)
flag = numpy.ma.masked_equal(ds.series[ThisOne]["Flag"],0)
mask = data.mask&flag.mask
index = numpy.ma.where(mask==True)[0]
ds.series[ThisOne]["Flag"][index] = numpy.int32(8)
#index = numpy.where((abs(data-float(-9999))<c.eps)&(abs(flag-float(0))<c.eps))[0]
#if len(index)!=0:
# ds.series[ThisOne]["Flag"][index] = numpy.int32(8)
In [5]:
ncname_3d="/home/peter/OzFlux/Sites/HowardSprings/Data/Processed/2012/HowardSprings_2012_L3.nc"
In [6]:
ncname_1d="/home/peter/OzFlux/Sites/Calperum/Data/Processed/2012/Calperum_2012_L3.nc"
In [38]:
ds_3d=nc_read_series_3d(ncname_3d)
In [8]:
ds_1d=nc_read_series_1d(ncname_1d)
In [39]:
ds_3d.series["Fsd"]["Data"][5]=numpy.float64(-9999)
In [41]:
%timeit CheckQCFlags(ds_3d)
In [42]:
print ds_3d.series["Fsd"]["Data"][0:10],ds_3d.series["Fsd"]["Flag"][0:10]
In [ ]: