Skip to content
Snippets Groups Projects
Commit a9208947 authored by Peters, Wouter's avatar Peters, Wouter
Browse files

time averaging now included for flux files

parent 22fa0fdc
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
da/analysis/expand_fluxes.py
+ 62
− 161
View file @ a9208947
......@@ -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)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment