From a92089475a99a6e7832747bab9fea8b7e8fef574 Mon Sep 17 00:00:00 2001
From: Wouter Peters <wouter.peters@wur.nl>
Date: Thu, 26 Apr 2012 10:01:16 +0000
Subject: [PATCH] time averaging now included for flux files
---
da/analysis/expand_fluxes.py | 223 ++++++++++-------------------------
1 file changed, 62 insertions(+), 161 deletions(-)
diff --git a/da/analysis/expand_fluxes.py b/da/analysis/expand_fluxes.py
index 4eb200cc..418abb16 100755
--- a/da/analysis/expand_fluxes.py
+++ b/da/analysis/expand_fluxes.py
@@ -666,64 +666,7 @@ def SaveWeeklyAvgExtTCData(DaCycle):
return saveas
-def SaveEcoDataExt(rundat):
- """ Function SaveEcoDataExt saves surface flux data to NetCDF files for extended ecoregions
-
- *** Inputs ***
- rundat : a RunInfo object
-
- *** Outputs ***
- NetCDF file containing n-hourly global surface fluxes per TransCom region
-
- *** Example ***
- ./expand_savestate project=enkf_release sd=20000101 ed=20010101 """
- from ecotools import xte
- from da.analysis.tools_transcom import LookUpName
- import pycdf as CDF
-
- infile=os.path.join(rundat.outputdir,'data_eco_weekly','ecofluxes.nc')
- if not os.path.exists(infile):
- print "Needed input file (%s) does not exist yet, please create weekly ecoflux files first, returning..."%infile
- return None
-
- # Create NetCDF output file
- #
- saveas=os.path.join(rundat.outputdir,'data_eco_weekly','eco_extfluxes.nc')
- ncf=CT_CDF(saveas,'create')
- dimdate=ncf.AddDateDim()
- dimidateformat=ncf.AddDateDimFormat()
- dimregs = ncf.AddRegionDim(type='tc_ext')
-
-
- ncf_in=CDF.CDF(infile)
- vardict=ncf_in.variables()
- for vname, vprop in vardict.iteritems():
- dims=()
- data=array(ncf_in.var(vname)[:])
- atts=ncf_in.var(vname).attributes()
- if vname not in ['date','idate']:
- if 'cov' in vname:
- data=xte(data[:,0:rundat.n_land,0:rundat.n_land],cov=True)
- dims=dimdate+dimregs+dimregs
- else:
- data=xte(data[:,0:rundat.n_land])
- dims=dimdate+dimregs
- elif vname=='date': dims=dimdate
- elif vname=='idate': dims=dimdate+dimidateformat
- else:
- print 'Dataset with unknown dimensions encountered in file: %s'%vname
-
- savedict=ncf.StandardVar(varname=vname)
- savedict['units']=atts['units']
- savedict['values']=data.tolist()
- savedict['dims']=dims
- ncf.AddData(savedict,nsets=data.shape[0])
-
- ncf.close()
-
- return saveas
-
-def SaveTimeAvgData(rundat,infile,avg='monthly'):
+def SaveTimeAvgData(DaCycle,infile,avg='monthly'):
""" Function saves time mean surface flux data to NetCDF files
*** Inputs ***
@@ -738,14 +681,18 @@ def SaveTimeAvgData(rundat,infile,avg='monthly'):
import da.analysis.tools_time as timetools
import da.tools.io4 as io
+ if 'weekly' in infile: intime = 'weekly'
+ if 'monthly' in infile: intime = 'monthly'
+ if 'yearly' in infile: intime = 'yearly'
+
dirname,filename = os.path.split(infile)
- dirname = CreateDirs(os.path.join(rundat.outputdir,dirname.replace('weekly',avg) ) )
+ outdir = CreateDirs(os.path.join(DaCycle['dir.analysis'],dirname.replace(intime,avg)))
dectime0 = date2num(datetime(2000,1,1))
# Create NetCDF output file
#
- saveas = os.path.join(dirname,filename)
+ saveas = os.path.join(outdir,filename)
ncf = io.CT_CDF(saveas,'create')
dimdate = ncf.AddDateDim()
#
@@ -763,68 +710,49 @@ def SaveTimeAvgData(rundat,infile,avg='monthly'):
date = file.GetVariable('date')
globatts = file.ncattrs()
- time = [datetime(2000,1,1)+timedelta(days=d) for d in date]
+ for att in globatts:
+ attval = file.getncattr(att)
+ if not att in ncf.ncattrs():
+ ncf.setncattr(att,attval)
-#
-# Add global attributes, sorted
-#
- keys = globatts
- keys.sort()
- for k in keys:
- setattr(ncf,k,k)
+ time = [datetime(2000,1,1)+timedelta(days=d) for d in date]
# loop over datasets in infile, skip idate and date as we will make new time axis for the averaged data
- for sds in datasets:
- if sds in ['idate','date'] : continue
- print '['
+ for sds in ['date']+datasets:
# get original data
- data = array(file.GetVariable(sds))
+ data = file.GetVariable(sds)
varatts = file.variables[sds].ncattrs()
vardims = file.variables[sds].dimensions
#
# Depending on dims of input dataset, create dims for output dataset. Note that we add the new dimdate now.
#
- dims=()
for d in vardims:
- if 'latitude' in d:
- dimgrid = ncf.AddLatLonDim()
- dims += (dimgrid[0],)
- if 'longitude' in d:
- dimgrid = ncf.AddLatLonDim()
- dims += (dimgrid[1],)
- for type in ['eco','eco_ext','tc','tc_ext','olson']:
- if 'regions_%s' % type == d:
- dimregs = ncf.AddRegionDim(type)
- dims += dimregs
-#
-# If the variable name from the infile is not defined in the standard, simply copy attributes from the old to the new sds
-# Do not copy over the grid information
-#
-
- if ncf.StandardVar(sds)['name'] == 'unknown':
-
- savedict = ncf.StandardVar(sds)
- savedict['name'] = sds
- savedict['values'] = data.tolist()
- savedict['dims'] = dims
- savedict['units'] = file.variables[sds].units
- savedict['long_name'] = file.variables[sds].long_name
- savedict['comment'] = file.variables[sds].comment
- savedict['standard_name'] = file.variables[sds].standard_name
- savedict['count'] =0
- ncf.AddData(savedict)
- sys.stdout.write('.')
- sys.stdout.flush()
-
+ if 'date' in d:
+ continue
+ if d in ncf.dimensions.keys():
+ pass
+ else:
+ dim = ncf.createDimension(d,size=len(file.dimensions[d]))
+
+ savedict = ncf.StandardVar(sds)
+ savedict['name'] = sds
+ savedict['dims'] = vardims
+ savedict['units'] = file.variables[sds].units
+ savedict['long_name'] = file.variables[sds].long_name
+ savedict['comment'] = file.variables[sds].comment
+ savedict['standard_name'] = file.variables[sds].standard_name
+ savedict['count'] = 0
+
+ if not 'date' in vardims:
+ savedict['values'] = data
+ ncf.AddData(savedict)
else:
- if sds in ['latitude','longitude']: continue
-
if avg == 'monthly':
time_avg, data_avg = timetools.monthly_avg(time,data)
elif avg == 'seasonal':
@@ -841,34 +769,27 @@ def SaveTimeAvgData(rundat,infile,avg='monthly'):
count=-1
for dd,data in zip(time_avg,data_avg):
count=count+1
-
- if not ncf.has_date(date2num(dd)-dectime0):
-
- savedict=ncf.StandardVar('date')
- savedict['values']=date2num(dd)-dectime0
- savedict['dims']=dimdate
- savedict['count']=count
- ncf.AddData(savedict)
-
- savedict=ncf.StandardVar(sds)
- savedict['values']=data.tolist()
- savedict['units']=file.variables[sds].units
- if 'cov' in sds:
- savedict['dims']=dimdate+dims+dims
+ if sds == 'date':
+ savedict['values'] = date2num(dd)-dectime0
else:
- savedict['dims']=dimdate+dims
+ savedict['values']=data
savedict['count']=count
- ncf.AddData(savedict)
+ ncf.AddData(savedict,silent=True)
sys.stdout.write('.')
- sys.stdout.flush()
- print ']'
+ sys.stdout.write('\n')
+ sys.stdout.flush()
# end NetCDF file access
file.close()
ncf.close()
+ msg = "------------------- Finished time averaging---------------------------------" ; logging.info(msg)
+
+
+ return saveas
+
if __name__ == "__main__":
import logging
@@ -892,43 +813,23 @@ if __name__ == "__main__":
StateVector = CtGriddedStateVector(DaCycle)
dummy = StateVector.Initialize()
+ while DaCycle['time.start'] < DaCycle['time.finish']:
+
+ savedas_1x1=SaveWeeklyAvg1x1Data(DaCycle, StateVector)
+ savedas_state=SaveWeeklyAvgStateData(DaCycle, StateVector)
+ savedas_tc=SaveWeeklyAvgTCData(DaCycle, StateVector)
+ savedas_tcext=SaveWeeklyAvgExtTCData(DaCycle)
+
+ DaCycle.AdvanceCycleTimes()
+
+ StateVector = None # free memory
+
+ for avg in ['monthly','yearly','longterm']:
+
+ savedas_1x1 = SaveTimeAvgData(DaCycle,savedas_1x1,avg)
+ savedas_state = SaveTimeAvgData(DaCycle,savedas_state,avg)
+ savedas_tc = SaveTimeAvgData(DaCycle,savedas_tc,avg)
+ savedas_tcext = SaveTimeAvgData(DaCycle,savedas_tcext,avg)
- no, yes, all = range(3)
-
- proceed = no
-
- # first, parse results from the inverse part of the run
- if proceed != all:
- proceed = proceed_dialog('Create average 1x1 flux data files? [y|yes, n|no, a|all|yes-to-all ]')
- if proceed != no:
- while DaCycle['time.start'] < DaCycle['time.finish']:
-
- savedas=SaveWeeklyAvg1x1Data(DaCycle, StateVector)
- savedas=SaveWeeklyAvgStateData(DaCycle, StateVector)
- savedas=SaveWeeklyAvgTCData(DaCycle, StateVector)
- savedas=SaveWeeklyAvgExtTCData(DaCycle)
-
- DaCycle.AdvanceCycleTimes()
-
- sys.exit(2)
-
- a=SaveTimeAvgData(rundat,savedas,avg='monthly')
- a=SaveTimeAvgData(rundat,savedas,avg='yearly')
- a=SaveTimeAvgData(rundat,savedas,avg='longterm')
-
- if proceed != all:
- proceed = proceed_dialog('Create average tc flux data files? [y|yes, n|no, a|all|yes-to-all] ')
- if proceed != no:
- a=SaveTimeAvgData(rundat,savedas,avg='monthly')
- a=SaveTimeAvgData(rundat,savedas,avg='seasonal')
- a=SaveTimeAvgData(rundat,savedas,avg='yearly')
- a=SaveTimeAvgData(rundat,savedas,avg='longterm')
- if proceed != all:
- proceed = proceed_dialog('Create average extended tc flux data files? [y|yes, n|no, a|all|yes-to-all] ')
- if proceed != no:
- savedas=SaveTCDataExt(rundat)
- a=SaveTimeAvgData(rundat,savedas,avg='monthly')
- a=SaveTimeAvgData(rundat,savedas,avg='seasonal')
- a=SaveTimeAvgData(rundat,savedas,avg='yearly')
- a=SaveTimeAvgData(rundat,savedas,avg='longterm')
sys.exit(0)
+
--
GitLab