Skip to content
Snippets Groups Projects
Commit 974d6f44 authored by Arne Babenhauserheide's avatar Arne Babenhauserheide
Browse files

refactored analysis.expand_fluxes for better readability

parent 8b066da8
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
refactored/da/analysis/expand_fluxes.py
+ 133
− 181
View file @ 974d6f44
......@@ -26,6 +26,50 @@ def proceed_dialog(txt, yes=['y','yes'], all=['a','all', 'yes-to-all']):
return 2
return 0
def _getpriordata(DaCycle, StateVector, startdate, enddate, nlag, cyclelength):
"""Read gridmean and gridvariance from the oldest available nlag.
:returns: gridmean, gridvariance, used nlag."""
import logging
tried = []
for n in range(nlag,0,-1):
# TODO: This breaks, when enddate - startdate is not a multiple of the cyclelength.
priordate = enddate - timedelta(cyclelength*n)
savedir = DaCycle['dir.output'].replace(startdate.strftime('%Y%m%d'),priordate.strftime('%Y%m%d') )
filename = os.path.join(savedir,'savestate.nc')
tried.append(filename)
if os.path.exists(filename):
dummy = StateVector.ReadFromFile(filename)
gridmean,gridvariance = StateVector.StateToGrid(lag=n)
msg = 'Read prior dataset from file %s, sds %d: '%(filename,n) ; logging.debug(msg)
return gridmean, gridvariance, n
raise ValueError("No savestate.nc found in the possible paths:" + str(tried))
def _saveInto(ncf, standardvar, data, dims,
count=None, name=None, long_name=None,
units=None, comment=None, nsets=None,
standard_name=None):
"""save the given data into the netcdf file ncf as a list.
:param standardvar: The name of the standard definition to use.
:param count: offset for writing the data (==previous data which should not be replaced)."""
savedict=ncf.StandardVar(varname=standardvar)
savedict["values"] = data
savedict["dims"] = dims
for key, val in [("count", count),
("name", name),
("long_name", long_name),
("units", units),
("standard_name", standard_name),
("comment", comment)]:
if val: savedict[key] = val
if nsets:
ncf.AddData(savedict, nsets=nsets)
return ncf.AddData(savedict)
def SaveWeeklyAvg1x1Data(DaCycle, StateVector):
"""
Function creates a NetCDF file with output on 1x1 degree grid. It uses the flux data written by the
......@@ -106,16 +150,7 @@ def SaveWeeklyAvg1x1Data(DaCycle, StateVector):
if prior:
qual_short='prior'
for n in range(nlag,0,-1):
priordate = enddate - timedelta(dt.days*n)
savedir = DaCycle['dir.output'].replace(startdate.strftime('%Y%m%d'),priordate.strftime('%Y%m%d') )
filename = os.path.join(savedir,'savestate.nc')
if os.path.exists(filename):
dummy = StateVector.ReadFromFile(filename)
gridmean,gridvariance = StateVector.StateToGrid(lag=n)
msg = 'Read prior dataset from file %s, sds %d: '%(filename,n) ; logging.debug(msg)
break
gridmean, gridvariance, choicelag = _getpriordata(DaCycle, StateVector, startdate, enddate, nlag, cyclelength=dt.days)
else:
qual_short='opt'
savedir = DaCycle['dir.output']
......@@ -136,49 +171,24 @@ def SaveWeeklyAvg1x1Data(DaCycle, StateVector):
#
# For each dataset, get the standard definitions from the module mysettings, add values, dimensions, and unlimited count, then write
#
savedict=ncf.StandardVar(varname='bio_flux_'+qual_short)
savedict['values']=biomapped.tolist()
savedict['dims']=dimdate+dimgrid
savedict['count']=next
ncf.AddData(savedict)
#
savedict=ncf.StandardVar(varname='ocn_flux_'+qual_short)
savedict['values']=oceanmapped.tolist()
savedict['dims']=dimdate+dimgrid
savedict['count']=next
ncf.AddData(savedict)
savedict=ncf.StandardVar(varname='bio_flux_%s_cov'%qual_short)
savedict['values']=biovarmapped.tolist()
savedict['dims']=dimdate+dimgrid
savedict['count']=next
ncf.AddData(savedict)
#
savedict=ncf.StandardVar(varname='ocn_flux_%s_cov'%qual_short)
savedict['values']=oceanvarmapped.tolist()
savedict['dims']=dimdate+dimgrid
savedict['count']=next
ncf.AddData(savedict)
for name, data in [('bio_flux_'+qual_short, biomapped),
('ocn_flux_'+qual_short, oceanmapped),
('bio_flux_%s_cov'%qual_short, biovarmapped),
('ocn_flux_%s_cov'%qual_short, oceanvarmapped)]:
_saveInto(ncf, name, data.tolist(),
dims=dimdate+dimgrid, count=next)
# End prior/posterior block
savedict=ncf.StandardVar(varname='fire_flux_imp')
savedict['values']=fire.tolist()
savedict['dims']=dimdate+dimgrid
savedict['count']=next
ncf.AddData(savedict)
#
savedict=ncf.StandardVar(varname='fossil_flux_imp')
savedict['values']=fossil.tolist()
savedict['dims']=dimdate+dimgrid
savedict['count']=next
ncf.AddData(savedict)
#
savedict=ncf.StandardVar(varname='date')
savedict['values']=date2num(startdate)-dectime0+dt.days/2.0
savedict['dims']=dimdate
savedict['count']=next
ncf.AddData(savedict)
for name, data in [('fire_flux_imp', fire),
('fossil_flux_imp', fossil)]:
_saveInto(ncf, name, data.tolist(),
dims=dimdate+dimgrid, count=next)
_saveInto(ncf, 'date', date2num(startdate)-dectime0+dt.days/2.0,
dims=dimdate, count=next)
sys.stdout.write('.')
sys.stdout.flush()
......@@ -286,16 +296,7 @@ def SaveWeeklyAvgStateData(DaCycle, StateVector):
if prior:
qual_short='prior'
for n in range(nlag,0,-1):
priordate = enddate - timedelta(dt.days*n)
savedir = DaCycle['dir.output'].replace(startdate.strftime('%Y%m%d'),priordate.strftime('%Y%m%d') )
filename = os.path.join(savedir,'savestate.nc')
if os.path.exists(filename):
dummy = StateVector.ReadFromFile(filename)
choicelag = n
msg = 'Read prior dataset from file %s, sds %d: '%(filename,n) ; logging.debug(msg)
break
dummy, dummy, choicelag = _getpriordata(DaCycle, StateVector, startdate, enddate, nlag, cyclelength=dt.days)
else:
qual_short='opt'
savedir = DaCycle['dir.output']
......@@ -309,83 +310,50 @@ def SaveWeeklyAvgStateData(DaCycle, StateVector):
#
data = StateVector.EnsembleMembers[choicelag-1][0].ParameterValues*vectorbio # units of mole region-1 s-1
savedict = ncf.StandardVar(varname='bio_flux_%s'%qual_short)
savedict['values'] = data
savedict['dims'] = dimdate+dimregs
savedict['count'] = next
ncf.AddData(savedict)
_saveInto(ncf, 'bio_flux_%s'%qual_short, data,
dims=dimdate+dimregs, count=next)
members = StateVector.EnsembleMembers[choicelag-1]
deviations = np.array([mem.ParameterValues*data for mem in members])
savedict=ncf.StandardVar(varname='bio_flux_%s_cov'%qual_short)
savedict['values']=deviations.tolist()
savedict['dims']=dimdate+dimmembers+dimregs
savedict['comment']="This is the matrix square root, use (M x M^T)/(nmembers-1) to make covariance"
savedict['units']="mol region-1 s-1"
savedict['count']=next
ncf.AddData(savedict)
savedict=ncf.StandardVar('unknown')
savedict['name'] = 'bio_flux_%s_std'%qual_short
savedict['long_name'] = 'Biosphere flux standard deviation, %s'%qual_short
savedict['values']=deviations.std(axis=0)
savedict['dims']=dimdate+dimregs
savedict['comment']="This is the standard deviation on each parameter"
savedict['units']="mol region-1 s-1"
savedict['count']=next
ncf.AddData(savedict)
_saveInto(ncf, 'bio_flux_%s_cov'%qual_short, deviations.tolist(),
dims=dimdate+dimmembers+dimregs, count=next,
units="mol region-1 s-1",
comment="This is the matrix square root, use (M x M^T)/(nmembers-1) to make covariance")
_saveInto(ncf, 'unknown', deviations.std(axis=0),
dims=dimdate+dimregs, count=next,
name = 'bio_flux_%s_std'%qual_short,
long_name = 'Biosphere flux standard deviation, %s'%qual_short,
comment = "This is the standard deviation on each parameter",
units = "mol region-1 s-1")
data = StateVector.EnsembleMembers[choicelag-1][0].ParameterValues*vectorocn # units of mole region-1 s-1
savedict=ncf.StandardVar(varname='ocn_flux_%s'%qual_short)
savedict['values']=data
savedict['dims']=dimdate+dimregs
savedict['count']=next
ncf.AddData(savedict)
_saveInto(ncf, 'ocn_flux_%s'%qual_short, data,
dims=dimdate+dimregs, count=next)
deviations = np.array([mem.ParameterValues*data for mem in members])
savedict=ncf.StandardVar(varname='ocn_flux_%s_cov'%qual_short)
savedict['values']=deviations.tolist()
savedict['dims']=dimdate+dimmembers+dimregs
savedict['comment']="This is the matrix square root, use (M x M^T)/(nmembers-1) to make covariance"
savedict['units']="mol region-1 s-1"
savedict['count']=next
ncf.AddData(savedict)
savedict=ncf.StandardVar('unknown')
savedict['name'] = 'ocn_flux_%s_std'%qual_short
savedict['long_name'] = 'Ocean flux standard deviation, %s'%qual_short
savedict['values']=deviations.std(axis=0)
savedict['dims']=dimdate+dimregs
savedict['comment']="This is the standard deviation on each parameter"
savedict['units']="mol region-1 s-1"
savedict['count']=next
ncf.AddData(savedict)
data = vectorfire
savedict=ncf.StandardVar(varname='fire_flux_imp')
savedict['values']=data
savedict['dims']=dimdate+dimregs
savedict['count']=next
ncf.AddData(savedict)
data = vectorfossil
savedict=ncf.StandardVar(varname='fossil_flux_imp')
savedict['values']=data
savedict['dims']=dimdate+dimregs
savedict['count']=next
ncf.AddData(savedict)
savedict=ncf.StandardVar(varname='date')
savedict['values'] = ncfdate
savedict['dims'] = dimdate
savedict['count'] = next
ncf.AddData(savedict)
_saveInto(ncf, 'ocn_flux_%s_cov'%qual_short, deviations.tolist(),
dims=dimdate+dimmembers+dimregs, count=next,
comment="This is the matrix square root, use (M x M^T)/(nmembers-1) to make covariance",
units="mol region-1 s-1")
_saveInto(ncf, 'unknown', deviations.std(axis=0),
dims=dimdate+dimregs, count=next,
name = 'ocn_flux_%s_std'%qual_short,
long_name = 'Ocean flux standard deviation, %s'%qual_short,
comment = "This is the standard deviation on each parameter",
units = "mol region-1 s-1")
for name, data in [('fire_flux_imp', vectorfire),
('fossil_flux_imp', vectorfossil)]:
_saveInto(ncf, name, data,
dims=dimdate+dimregs, count=next)
_saveInto(ncf, 'date', ncfdate,
dims=dimdate, count=next)
sys.stdout.write('.')
sys.stdout.flush()
......@@ -485,11 +453,8 @@ def SaveWeeklyAvgTCData(DaCycle, StateVector):
# First add the date for this cycle to the file, this grows the unlimited dimension
savedict = ncf.StandardVar(varname='date')
savedict['values'] = ncfdate
savedict['dims'] = dimdate
savedict['count'] = index
dummy = ncf.AddData(savedict)
_saveInto(ncf, 'date', ncfdate,
dims=dimdate, count=index)
# Now convert other variables that were inside the flux_1x1 file
......@@ -498,29 +463,25 @@ def SaveWeeklyAvgTCData(DaCycle, StateVector):
data = ncf_in.GetVariable(vname)[index]
savedict = ncf.StandardVar(varname=vname)
if vname=='latitude': continue
elif vname=='longitude': continue
elif vname=='date': continue
elif vname=='idate': continue
elif 'cov' in vname:
cov = data.transpose().dot(data)
#AB? Only works with recent version of numpy, not with 1.4.1 installed on the suns
#cov = data.transpose().dot(data)
cov = np.dot(data.transpose(),data)
tcdata = StateVector.VectorToTC(vectordata=cov,cov=True) # vector to TC
savedict['units'] = '[mol/region/s]**2'
savedict['dims'] = dimdate+dimregs+dimregs
units = '[mol/region/s]**2'
dims = dimdate+dimregs+dimregs
else:
tcdata = StateVector.VectorToTC(vectordata=data) # vector to TC
dims = dimdate+dimregs
units = 'mol/region/s'
savedict['dims'] = dimdate+dimregs
savedict['units'] = 'mol/region/s'
savedict['values'] = tcdata
ncf.AddData(savedict)
_saveInto(ncf, vname, tcdata,
dims=dims, units=units)
ncf_in.close()
ncf.close()
......@@ -582,12 +543,10 @@ def SaveTCDataExt(rundat):
dims=dimdate+dimregs
print vname,data.shape
savedict = ncf.StandardVar(varname=vname)
savedict['values'] = data.tolist()
savedict['dims'] = dims
savedict['units'] = 'mol/region/s'
savedict['count'] = 0
ncf.AddData(savedict,nsets=data.shape[0])
_saveInto(ncf, vname, data.tolist(),
dims=dims, count=0,
units='mol/region/s',
nsets=data.shape[0])
ncf.close()
......@@ -640,11 +599,10 @@ def SaveEcoDataExt(rundat):
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])
_saveInto(ncf, vname, data.tolist(),
dims=dims,
units=atts['units'],
nsets=data.shape[0])
ncf.close()
......@@ -735,18 +693,16 @@ def SaveTimeAvgData(rundat,infile,avg='monthly'):
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()
invar = file.variables[sds]
_saveInto(ncf, sds, data.tolist(),
dims=dims, count=0,
name=sds,
units=invar.units,
long_name=invar.long_name,
comment=invar.comment,
standard_name=invar.standard_name)
sys.stdout.write('.')
sys.stdout.flush()
else:
......@@ -770,22 +726,18 @@ def SaveTimeAvgData(rundat,infile,avg='monthly'):
count=count+1
if not ncf.has_date(date2num(dd)-dectime0):
_saveInto(ncf, 'date', date2num(dd)-dectime0,
dims=dimdate, count=count)
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
units=file.variables[sds].units
if 'cov' in sds:
savedict['dims']=dimdate+dims+dims
dimstmp=dimdate+dims+dims
else:
savedict['dims']=dimdate+dims
savedict['count']=count
ncf.AddData(savedict)
dimstmp=dimdate+dims
_saveInto(ncf, sds, data.tolist(),
dims=dimstmp, count=count,
units=units)
sys.stdout.write('.')
sys.stdout.flush()
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment