Commit 021c4336 authored by Florentie, Liesbeth's avatar Florentie, Liesbeth
Browse files

small changes to facilitate use with focal retrieval

parent 131c34c2
......@@ -119,11 +119,11 @@ class TotalColumnObservations(Observations):
# Get names of additional variables required for sampling that should be read from sounding files
# Currently implemented options are: 'surface_pressure', 'pressure_weight', 'SigmaB'
# MAYBE IMPLEMENT ALSO OPTION TO ENTER SATELLITE, AND BASED ON THAT SET DEFAULT ADDITIONAL VARIABLES?
if 'obs.column.extra.variables' in dacycle.dasystem.keys():
self.additional_variables = map(str.strip, dacycle.dasystem['obs.column.extra.variables'].split(','))
logging.debug('Additional sounding variables that will be read from file: %s' %self.additional_variables)
else:
self.additional_variables = []
# if 'obs.column.extra.variables' in dacycle.dasystem.keys():
# self.additional_variables = map(str.strip, dacycle.dasystem['obs.column.extra.variables'].split(','))
# logging.debug('Additional sounding variables that will be read from file: %s' %self.additional_variables)
# else:
self.additional_variables = []
# Model data mismatch approach
self.mdm_calculation = dacycle.dasystem.get('obs.column.mdm.calculation')
......@@ -137,13 +137,14 @@ class TotalColumnObservations(Observations):
# Scaling function for mdm to tweak relative weight obs w.r.t. surface obs: none, gaussian (gauss) or inverse gaussian (inv_gauss)
if 'obs.column.mdm.scaling.function' in dacycle.dasystem.keys():
self.mdm_scaling_fnc = dacycle.dasystem.get('obs.column.mdm.scaling.function').strip()
self.mdm_scaling_peak = float(dacycle.dasystem.get('obs.column.mdm.scaling.peak'))
self.mdm_scaling_mean = float(dacycle.dasystem.get('obs.column.mdm.scaling.mean'))
self.mdm_scaling_std = float(dacycle.dasystem.get('obs.column.mdm.scaling.std'))
logging.info('MDM uniform weight = %s, MDM scaling function = %s with parameters %s, %s, %s' %(str(self.mdm_weight), \
self.mdm_scaling_fnc,str(self.mdm_scaling_peak),str(self.mdm_scaling_mean),str(self.mdm_scaling_std)))
if self.mdm_scaling_fnc in ['gauss', 'inv_gauss']:
self.mdm_scaling_peak = float(dacycle.dasystem.get('obs.column.mdm.scaling.peak'))
self.mdm_scaling_mean = float(dacycle.dasystem.get('obs.column.mdm.scaling.mean'))
self.mdm_scaling_std = float(dacycle.dasystem.get('obs.column.mdm.scaling.std'))
logging.info('MDM scaling parameters: %s, %s, %s' %(str(self.mdm_scaling_peak),str(self.mdm_scaling_mean),str(self.mdm_scaling_std)))
else:
self.mdm_scaling_fnc = 'none'
logging.info('MDM uniform weight = %s, MDM scaling function = %s' %(str(self.mdm_weight), self.mdm_scaling_fnc))
# Path to file with observation error settings for column observations
# Currently the same settings for all assimilated retrieval products: should this be one file per product?
......@@ -212,36 +213,40 @@ class TotalColumnObservations(Observations):
lons = ncf.get_variable('longitude').take(subselect, axis=0)
obs = ncf.get_variable('obs').take(subselect, axis=0)
unc = ncf.get_variable('uncertainty').take(subselect, axis=0)
prior = ncf.get_variable('prior').take(subselect, axis=0)
prior_profile = ncf.get_variable('prior_profile').take(subselect, axis=0)
av_kernel = ncf.get_variable('averaging_kernel').take(subselect, axis=0)
pressure = ncf.get_variable('pressure_levels').take(subselect, axis=0)
dates = ncf.get_variable('date').take(subselect, axis=0)
dates = array([dtm.datetime(*d) for d in dates])
av_kernel = ncf.get_variable('averaging_kernel').take(subselect, axis=0)
prior_profile = ncf.get_variable('prior_profile').take(subselect, axis=0)
pressure = ncf.get_variable('pressure_levels').take(subselect, axis=0)
if 'prior' in ncf.variables:
prior = ncf.get_variable('prior').take(subselect, axis=0)
self.additional_variables.append('prior')
else: prior = [float('nan')]*len(ids)
if 'pressure_weight' in self.additional_variables and 'pressure_weight' in ncf.variables:
if 'pressure_weight' in ncf.variables:
h = ncf.get_variable('pressure_weight').take(subselect, axis=0)
else:
h = [float('nan')]*len(ids)
self.additional_variables.append('pressure_weight')
else: h = [float('nan')]*len(ids)
if ('surface_pressure' in self.additional_variables or level_def == 'layer_average') and 'surface_pressure' in ncf.variables:
if 'surface_pressure' in ncf.variables:
psurf = ncf.get_variable('surface_pressure').take(subselect, axis=0)
if 'surface_pressure' not in self.additional_variables: self.additional_variables.append('surface_pressure')
else:
psurf = [float('nan')]*len(ids)
self.additional_variables.append('surface_pressure')
else: psurf = [float('nan')]*len(ids)
if 'SigmaB' in self.additional_variables and 'SigmaB' in ncf.variables:
if 'SigmaB' in ncf.variables:
self.sigmab = ncf.get_variable('SigmaB')
self.additional_variables.append('SigmaB')
ncf.close()
# Add samples to datalist
# Note that the mdm is initialized here equal to the measurement uncertainty. This value is used in add_model_data_mismatch to calculate the mdm including model error!
# Note that the mdm is initialized here equal to the measurement uncertainty. This value is used in add_model_data_mismatch to calculate the mdm including model error
for n in range(len(ids)):
# Check for every sounding if time is between start and end time (relevant for first and last days of window)
if self.startdate <= dates[n] <= self.enddate:
self.datalist.append(TotalColumnSample(ids[n], code, dates[n], obs[n], None, lats[n], lons[n], unc[n], prior[n], prior_profile[n,:], \
av_kernel[n,:], pressure=pressure[n,:],level_def=level_def,psurf=psurf[n],h=h[n]))
self.datalist.append(TotalColumnSample(ids[n], code, dates[n], obs[n], None, lats[n], lons[n], unc[n], prior=prior[n], prior_profile=prior_profile[n,:], \
av_kernel=av_kernel[n,:], pressure=pressure[n,:],level_def=level_def,psurf=psurf[n],h=h[n]))
logging.debug("Added %d observations to the Data list" % (len(self.datalist)-len_init))
......@@ -305,7 +310,7 @@ class TotalColumnObservations(Observations):
# elif self.mdm_calculation == 'empirical':
# obs.mdm = ...
else: # assume general model uncertainty of 1 ppm (arbitrary value)
obs.mdm = self.mdm_weight * ( obs.mdm*obs.mdm + 1**2 )**0.5
obs.mdm = self.mdm_weight * ( obs.mdm*obs.mdm + 1.0**2 )**0.5
# latitude-dependent scaling of mdm to account for varying number of surface observations
# when combining column obs with flask obs
......@@ -409,9 +414,13 @@ class TotalColumnObservations(Observations):
logging.debug('Creating new observations file for ObservationOperator (%s) containing %d observations' % (obsinputfile,len(self.datalist)))
dimsoundings = f.add_dim('soundings', len(self.datalist))
dimlevels = f.add_dim('levels', self.getvalues('pressure').shape[1])
dimdate = f.add_dim('epoch_dimension', 7)
dimchar = f.add_dim('char', 20)
dimlevels = f.add_dim('levels', self.getvalues('pressure').shape[1])
if self.getvalues('av_kernel').shape[1] != self.getvalues('pressure').shape[1]:
dimlayers = f.add_dim('layers', self.getvalues('pressure').shape[1] - 1)
layers = True
else: layers = False
savedict = io.std_savedict.copy()
savedict['dtype'] = "int64"
......@@ -441,23 +450,23 @@ class TotalColumnObservations(Observations):
f.add_data(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "prior"
savedict['dims'] = dimsoundings
savedict['values'] = self.getvalues('prior').tolist()
savedict['name'] = "averaging_kernel"
if layers:
savedict['dims'] = dimsoundings + dimlayers
else:
savedict['dims'] = dimsoundings + dimlevels
savedict['values'] = self.getvalues('av_kernel').tolist()
f.add_data(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "prior_profile"
savedict['dims'] = dimsoundings + dimlevels
if layers:
savedict['dims'] = dimsoundings + dimlayers
else:
savedict['dims'] = dimsoundings + dimlevels
savedict['values'] = self.getvalues('prior_profile').tolist()
f.add_data(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "averaging_kernel"
savedict['dims'] = dimsoundings + dimlevels
savedict['values'] = self.getvalues('av_kernel').tolist()
f.add_data(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "pressure_levels"
savedict['dims'] = dimsoundings + dimlevels
......@@ -477,6 +486,13 @@ class TotalColumnObservations(Observations):
savedict['values'] = self.getvalues('obs').tolist()
f.add_data(savedict)
if 'prior' in self.additional_variables:
savedict = io.std_savedict.copy()
savedict['name'] = "prior"
savedict['dims'] = dimsoundings
savedict['values'] = self.getvalues('prior').tolist()
f.add_data(savedict)
if 'surface_pressure' in self.additional_variables:
savedict = io.std_savedict.copy()
savedict['name'] = "psurf"
......@@ -494,7 +510,10 @@ class TotalColumnObservations(Observations):
if 'pressure_weight' in self.additional_variables:
savedict = io.std_savedict.copy()
savedict['name'] = "pressure_weight"
savedict['dims'] = dimsoundings + dimlevels
if layers:
savedict['dims'] = dimsoundings + dimlayers
else:
savedict['dims'] = dimsoundings + dimlevels
savedict['values'] = self.getvalues('h').tolist()
f.add_data(savedict)
......
Supports Markdown
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