Commit bb5bb773 authored by Ingrid Luijkx's avatar Ingrid Luijkx
Browse files

Added mdm scaling with number of observations per day (except for aircraft)

parent a30cdb97
......@@ -14,7 +14,7 @@ import logging
import datetime as dtm
from string import strip
from numpy import array, logical_and
from numpy import array, logical_and, sqrt
sys.path.append(os.getcwd())
sys.path.append('../../')
......@@ -85,9 +85,15 @@ class ObsPackObservations(Observations):
subselect = logical_and(dates >= self.startdate , dates <= self.enddate).nonzero()[0]
dates = dates.take(subselect, axis=0)
obspacknum = ncf.get_variable('obspack_num').take(subselect) # or should we propagate obs_num which is not unique across datasets??
obspackid = ncf.get_variable('obspack_id').take(subselect, axis=0)
if 'merge_num' in ncf.variables:
obspacknum = ncf.get_variable('merge_num').take(subselect)
else:
obspacknum = ncf.get_variable('obspack_num').take(subselect)
if 'ccggAllData' in ncfile:
obspackid = ncf.get_variable('id').take(subselect, axis=0)
else:
obspackid = ncf.get_variable('obspack_id').take(subselect, axis=0)
obspackid = [s.tostring().lower() for s in obspackid]
obspackid = map(strip, obspackid)
datasetname = ncfile # use full name of dataset to propagate for clarity
......@@ -325,7 +331,7 @@ class ObsPackObservations(Observations):
identifier, incalt = value.split(';')
site_incalt[identifier.strip()] = (int(incalt))
for obs in self.datalist: # loop over all available data points
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
if obs.flag == 1: # flag is taken from the gv+ datasets: 1=background/representative, 0=local.
......@@ -334,6 +340,8 @@ class ObsPackObservations(Observations):
obs.flag = 99 # 99 means: do-not-use
else: obs.flag = 99
for obs in self.datalist: # second loop over all available data points to set mdm
identifier = obs.code
species, site, method, lab, datasetnr = identifier.split('_')
......@@ -345,8 +353,12 @@ class ObsPackObservations(Observations):
obs.may_reject = site_info[identifier]['may_reject']
obs.flag = 99
else:
logging.debug("Observation found (%s, %d), mdm used is: %0.2f." % (identifier, obs.id, site_info[identifier]['error']))
obs.mdm = site_info[identifier]['error'] * self.global_R_scaling
if site_info[identifier]['category'] == 'aircraft':
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 = site_info[identifier]['error'] * sqrt(nr_obs_per_day) * self.global_R_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