Commit 83d93d63 authored by Florentie, Liesbeth's avatar Florentie, Liesbeth
Browse files

fix consistency units of samples

parent 1c13b8c2
......@@ -59,6 +59,12 @@ class ObsPackObservations(Observations):
self.obspack_dir = op_dir
self.obspack_id = op_id
# observations & mdm scaling. Choose such that unit is ppm
try:
self.obs_scaling = float(dacycle.dasystem['obs.scaling.factor'])
except KeyError:
self.obs_scaling = 1.0
self.datalist = []
if not os.path.exists(dacycle.dasystem['obs.sites.rc']):
......@@ -148,6 +154,9 @@ class ObsPackObservations(Observations):
flags = ncf.get_variable('obs_flag').take(subselect, axis=0)
ncf.close()
# convert obs unit if required
obs = obs * self.obs_scaling
# only add observations with flag = 1
counter = 0
for n in range(len(dates)):
......@@ -173,7 +182,8 @@ class ObsPackObservations(Observations):
ncf = io.ct_read(filename, method='read')
ids = ncf.get_variable('obs_num')
simulated = ncf.get_variable('flask')
# save samples in ppm, not mol/mol (TM5 output)
simulated = ncf.get_variable('flask') * 1e6
ncf.close()
logging.info("Successfully read data from model sample file (%s)" % filename)
......@@ -377,7 +387,7 @@ class ObsPackObservations(Observations):
del key, value
for obs in self.datalist: # first loop over all available data points to set flags correctly
obs.mdm = 1000.0 # default is very high model-data-mismatch, until explicitly set by script
obs.mdm = 1000.0 * self.obs_scaling # default is very high model-data-mismatch, until explicitly set by script
if obs.flag == 1: # flag is taken from the gv+ datasets: 1=background/representative, 0=local.
obs.flag = 0
elif obs.flag == 0:
......@@ -404,7 +414,7 @@ class ObsPackObservations(Observations):
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
obs.mdm = self.mdm_weight * site_info[identifier]['error'] * sqrt(nr_obs_per_day) * self.global_R_scaling * self.obs_scaling
obs.may_localize = site_info[identifier]['may_localize']
obs.may_reject = site_info[identifier]['may_reject']
obs.flag = 0
......
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