Commit 131c34c2 authored by Florentie, Liesbeth's avatar Florentie, Liesbeth
Browse files

removed reading of all flask obs in advance step + added a zonal mdm scaling...

removed reading of all flask obs in advance step + added a zonal mdm scaling based on number of observations for column obs
parent 735ba963
......@@ -135,12 +135,15 @@ class TotalColumnObservations(Observations):
self.mdm_weight = float(dacycle.dasystem['obs.column.mdm.weight'])
# Scaling function for mdm to tweak relative weight obs w.r.t. surface obs: none, gaussian (gauss) or inverse gaussian (inv_gauss)
self.mdm_scaling_fnc = dacycle.dasystem.get('obs.column.mdm.scaling.function')
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 '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)))
else:
self.mdm_scaling_fnc = 'none'
# 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?
......@@ -271,6 +274,24 @@ class TotalColumnObservations(Observations):
obs_data = rc.read(self.obs_file)
self.rejection_threshold = int(obs_data['obs.rejection.threshold'])
# if mdm scaling based on Nobs per zone is selected, calculate scaling factors
if self.mdm_scaling_fnc == 'nobs_zones':
nobs_south, nobs_tropic, nobs_north = 0.0, 0.0, 0.0
for obs in self.datalist:
if -90 <= obs.lat < -30:
nobs_south += 1.0
elif -30 <= obs.lat < 30:
nobs_tropic += 1.0
else:
nobs_north += 1.0
del obs
nobs_scaling_south = np.min([ np.sqrt(nobs_south / np.min([nobs_south,nobs_tropic,nobs_north])) , 5 ])
nobs_scaling_tropic = np.min([ np.sqrt(nobs_tropic / np.min([nobs_south,nobs_tropic,nobs_north])) , 5 ])
nobs_scaling_north = np.min([ np.sqrt(nobs_north / np.min([nobs_south,nobs_tropic,nobs_north])) , 5 ])
logging.info('mdm scaling factors based on Nobs per zone calculated: southern = %s, tropics = %s, northern = %s'
%(nobs_scaling_south,nobs_scaling_tropic,nobs_scaling_north))
# At this point mdm is set to the measurement uncertainty only, added in the add_observations function.
# Here this value is used to set the combined mdm by adding an estimate for the model uncertainty as a sum of squares.
if len(self.datalist) == 0: return
......@@ -288,10 +309,17 @@ class TotalColumnObservations(Observations):
# latitude-dependent scaling of mdm to account for varying number of surface observations
# when combining column obs with flask obs
if self.mdm_scaling_fnc is 'gauss':
if self.mdm_scaling_fnc == 'gauss':
obs.mdm = obs.mdm * self.gauss_correction(obs.lat)
elif self.mdm_scaling_fnc is 'inv_gauss':
elif self.mdm_scaling_fnc == 'inv_gauss':
obs.mdm = obs.mdm * self.inv_gauss_correction(obs.lat)
elif self.mdm_scaling_fnc == 'nobs_zones':
if -90 <= obs.lat < -30:
obs.mdm = obs.mdm * nobs_scaling_south
elif -30 <= obs.lat < 30:
obs.mdm = obs.mdm * nobs_scaling_tropic
else:
obs.mdm = obs.mdm * nobs_scaling_north
del obs
meanmdm = np.average(np.array( [obs.mdm for obs in self.datalist] ))
......
......@@ -115,7 +115,7 @@ class ObsPackObservations(Observations):
sitename, sitecategory = key, value
sitename = sitename.strip()
sitecategory = sitecategory.split()[0].strip().lower()
if (sitecategory != 'do-not-use') or advance:
if (sitecategory != 'do-not-use'):
valid_sites.append(sitename)
del key, value
......@@ -376,7 +376,7 @@ class ObsPackObservations(Observations):
sitename, sitecategory = key, value
sitename = sitename.strip()
sitecategory = sitecategory.split()[0].strip().lower()
if (sitecategory == 'do-not-use') and (not advance): continue
if (sitecategory == 'do-not-use'): continue # and (not advance): continue
site_info[sitename] = site_categories[sitecategory]
if 'site.move' in key:
identifier, latmove, lonmove = value.split(';')
......@@ -413,8 +413,9 @@ class ObsPackObservations(Observations):
nr_obs_per_day = 1
else:
nr_obs_per_day = len([c.code for c in self.datalist if c.code == obs.code and c.xdate.day == obs.xdate.day and c.flag == 0])
logging.debug("Observation found (%s, %d), mdm category is: %0.2f, scaled with number of observations per day (%i), final mdm applied is: %0.2f." % (identifier, obs.id, site_info[identifier]['error'],nr_obs_per_day,site_info[identifier]['error']*sqrt(nr_obs_per_day)))
obs.mdm = self.mdm_weight * site_info[identifier]['error'] * sqrt(nr_obs_per_day) * self.global_R_scaling * self.obs_scaling
logging.debug("Observation found (%s, %d), mdm category is: %0.2f, scaled with number of observations per day (%i), final mdm applied is: %0.2f." \
% (identifier, obs.id, site_info[identifier]['error'],nr_obs_per_day,obs.mdm)) #site_info[identifier]['error']*sqrt(nr_obs_per_day)))
obs.may_localize = site_info[identifier]['may_localize']
obs.may_reject = site_info[identifier]['may_reject']
obs.flag = 0
......@@ -432,7 +433,7 @@ class ObsPackObservations(Observations):
obs.height = obs.height + incalt
logging.warning("Observation location for (%s, %d), is moved by %3.2f meters in altitude" % (identifier, obs.id, incalt))
if advance == False and obs.flag > 80:
if obs.flag > 80: # and advance == False:
logging.debug('Dropped observation from datalist, as it is not to be assimilated')
else:
obs_to_keep.append(obs)
......
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