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

This routine now also opens the optimizer.nc file and retrieves the forecast...

This routine now also opens the optimizer.nc file and retrieves the forecast values. Then these are written to
the observation files too for later use.

Note that with the current processing chain, it is expected that the format of observations is ObsPack!

parent 3adbe550
Branches
No related tags found
No related merge requests found
......@@ -55,7 +55,7 @@ def WriteMixingRatios(DaCycle):
DaCycle['time.sample.stamp'] = "%s_%s"%(startdate.strftime("%Y%m%d%H"),enddate.strftime("%Y%m%d%H"),)
# Step (1): Get the sample output data file for this cycle
# Step (1): Get the posterior sample output data file for this cycle
infile = os.path.join(DaCycle['dir.output'],'sampleinfo_%s.nc'% DaCycle['time.sample.stamp'])
......@@ -63,7 +63,6 @@ def WriteMixingRatios(DaCycle):
obs_num = ncf_in.GetVariable('obs_num')
obs_val = ncf_in.GetVariable('observed')
mdm = ncf_in.GetVariable('modeldatamismatch')
simulated = ncf_in.GetVariable('modelsamples')
infilename= ncf_in.GetVariable('inputfilename')
......@@ -71,6 +70,33 @@ def WriteMixingRatios(DaCycle):
ncf_in.close()
# Step (2): Get the prior sample output data file for this cycle
infile = os.path.join(DaCycle['dir.output'],'optimizer.%s.nc'% startdate.strftime('%Y%m%d') )
ncf_fc_in = io.CT_Read(infile,'read')
fc_obs_num = ncf_fc_in.GetVariable('obspack_num')
fc_obs_val = ncf_fc_in.GetVariable('observed')
fc_simulated = ncf_fc_in.GetVariable('modelsamplesmean_prior')
fc_simulated_ens = ncf_fc_in.GetVariable('modelsamplesdeviations_prior')
fc_flag = ncf_fc_in.GetVariable('flag')
fc_r = ncf_fc_in.GetVariable('modeldatamismatchvariance').diagonal()
fc_hphtr = ncf_fc_in.GetVariable('totalmolefractionvariance').diagonal()
filesitecode = ncf_fc_in.GetVariable('sitecode')
fc_sitecodes = [join(s.compressed(),'') for s in filesitecode]
ncf_fc_in.close()
# Expand the list of input files with those available from the forecast list
infiles_rootdir = os.path.split(infiles[0])[0]
infiles.extend(os.path.join(infiles_rootdir,f+'.nc') for f in fc_sitecodes)
#Step (2): For each observation timeseries we now have data for, open it and fill with data
for orig_file in set(infiles):
if not os.path.exists(orig_file):
......@@ -124,7 +150,15 @@ def WriteMixingRatios(DaCycle):
dimmembers = ('nmembers',)
nmembers = len(dimmembers)
# Create empty arrays
# Create empty arrays for posterior samples, as well as for forecast sample statistics
savedict = io.std_savedict.copy()
savedict['name'] = "flag_forecast"
savedict['long_name'] = "flag_for_obs_model in forecast"
savedict['units'] = "None"
savedict['dims'] = dimid
savedict['comment'] = 'Flag (0/1/2/99) for observation value, 0 means okay, 1 means QC error, 2 means rejected, 99 means not sampled'
dummy = ncf_out.AddVariable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modeldatamismatch"
......@@ -134,6 +168,13 @@ def WriteMixingRatios(DaCycle):
savedict['comment'] = 'Variance of mole fractions resulting from model-data mismatch'
status = ncf_out.AddVariable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "totalmolefractionvariance_forecast"
savedict['long_name'] = "totalmolefractionvariance of forecast"
savedict['units'] = "[mol mol-1]^2"
savedict['dims'] = dimid
savedict['comment'] = 'Variance of mole fractions resulting from prior state and model-data mismatch'
dummy = ncf_out.AddVariable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modelsamplesmean"
......@@ -143,6 +184,14 @@ def WriteMixingRatios(DaCycle):
savedict['comment'] = 'simulated mixing ratios based on optimized state vector'
status = ncf_out.AddVariable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modelsamplesmean_forecast"
savedict['long_name'] = "mean modelsamples from forecast"
savedict['units'] = "mol mol-1"
savedict['dims'] = dimid
savedict['comment'] = 'simulated mixing ratios based on prior state vector'
status = ncf_out.AddVariable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modelsamplesstandarddeviation"
savedict['long_name'] = "standard deviaton of modelsamples over all ensemble members"
......@@ -151,6 +200,14 @@ def WriteMixingRatios(DaCycle):
savedict['comment'] = 'std dev of simulated mixing ratios based on optimized state vector'
status = ncf_out.AddVariable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modelsamplesstandarddeviation_forecast"
savedict['long_name'] = "standard deviaton of modelsamples from forecast over all ensemble members"
savedict['units'] = "mol mol-1"
savedict['dims'] = dimid
savedict['comment'] = 'std dev of simulated mixing ratios based on prior state vector'
status = ncf_out.AddVariable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modelsamplesensemble"
savedict['long_name'] = "modelsamples over all ensemble members"
......@@ -159,6 +216,14 @@ def WriteMixingRatios(DaCycle):
savedict['comment'] = 'ensemble of simulated mixing ratios based on optimized state vector'
status = ncf_out.AddVariable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modelsamplesensemble_forecast"
savedict['long_name'] = "modelsamples from forecast over all ensemble members"
savedict['units'] = "mol mol-1"
savedict['dims'] = dimid+dimmembers
savedict['comment'] = 'ensemble of simulated mixing ratios based on prior state vector'
status = ncf_out.AddVariable(savedict)
else:
msg = "Modifying existing file (%s) in the analysis directory"%copy_file ; logging.debug(msg)
......@@ -175,16 +240,17 @@ def WriteMixingRatios(DaCycle):
# Get all obs_nums related to this file, determine their indices in the local arrays
selected_obs_nums = [num for infile,num in zip(infiles,obs_num) if infile == orig_file]
selected_fc_obs_nums = [num for sitecode,num in zip(fc_sitecodes,fc_obs_num) if sitecode in orig_file]
# For each index, get the data and add to the file in the proper file index location
# Optimized data 1st: For each index, get the data and add to the file in the proper file index location
for num in selected_obs_nums:
model_index = obs_num.tolist().index(num)
file_index = file_obs_nums.tolist().index(num)
var = ncf_out.variables['modeldatamismatch']
var[file_index] = mdm[model_index]
#var = ncf_out.variables['modeldatamismatch'] # Take from optimizer.yyyymmdd.nc file instead
#var[file_index] = mdm[model_index]
var = ncf_out.variables['modelsamplesmean']
var[file_index] = simulated[model_index,0]
......@@ -195,6 +261,32 @@ def WriteMixingRatios(DaCycle):
var = ncf_out.variables['modelsamplesensemble']
var[file_index] = simulated[model_index,:]
# Now forecast data too: For each index, get the data and add to the file in the proper file index location
for num in selected_fc_obs_nums:
model_index = fc_obs_num.tolist().index(num)
file_index = file_obs_nums.tolist().index(num)
var = ncf_out.variables['modeldatamismatch']
var[file_index] = np.sqrt(fc_r[model_index])
var = ncf_out.variables['modelsamplesmean_forecast']
var[file_index] = fc_simulated[model_index]
var = ncf_out.variables['modelsamplesstandarddeviation_forecast']
var[file_index] = fc_simulated_ens[model_index,1:].std()
var = ncf_out.variables['modelsamplesensemble_forecast']
var[file_index] = fc_simulated_ens[model_index,:]
var = ncf_out.variables['totalmolefractionvariance_forecast']
var[file_index] = fc_hphtr[model_index]
var = ncf_out.variables['flag_forecast']
var[file_index] = fc_flag[model_index]
# close the file
status = ncf_out.close()
......@@ -213,11 +305,11 @@ if __name__ == "__main__":
logging.root.setLevel(logging.DEBUG)
DaCycle = CycleControl(args={'rc':'../../da.rc'})
DaCycle = CycleControl(args={'rc':'../../ctdas-od-gfed2-glb6x4-obspack-full.rc'})
DaCycle.Initialize()
DaCycle.ParseTimes()
DaSystem = CtDaSystem('../rc/carbontracker.rc')
DaSystem = CtDaSystem('../rc/carbontracker_ct09_opf.rc')
DaSystem.Initialize()
DaCycle.DaSystem = DaSystem
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment