Commit 0c6ceae1 authored by brunner's avatar brunner
Browse files

Working, 22 stations, different mdms

parent 6b5b4b0f
......@@ -336,7 +336,7 @@ def timehistograms_new(fig, infile, option='final'):
# Get a scaling factor for the x-axis range. Now we will include 5 standard deviations
sc = res.std()
print 'sc',sc
print('sc',sc)
# If there is too little data for a reasonable PDF, skip to the next value in the loop
if res.shape[0] < 10: continue
......
......@@ -97,10 +97,10 @@ def summarize_obs(analysisdir, printfmt='html'):
infiles = [os.path.join(mrdir, f) for f in mrfiles if f.endswith('.nc')]
if printfmt == 'tex':
print '\\begin{tabular*}{\\textheight}{l l l l r r r r}'
print 'Code & Name & Lat, Lon, Elev & Lab & N (flagged) & $\\sqrt{R}$ &Inn \\XS &Bias\\\\'
print '\hline\\\\ \n\multicolumn{8}{ c }{Semi-Continuous Surface Samples}\\\\[3pt] '
fmt = '%8s & ' + ' %55s & ' + '%20s &' + '%6s &' + ' %4d (%d) & ' + ' %5.2f & ' + ' %5.2f & ' + '%+5.2f \\\\'
print('\\begin{tabular*}{\\textheight}{l l l l r r r r}')
print('Code & Name & Lat, Lon, Elev & Lab & N (flagged) & $\\sqrt{R}$ &Inn \\XS &Bias\\\\')
print('\hline\\\\ \n\multicolumn{8}{ c }{Semi-Continuous Surface Samples}\\\\[3pt] ')
fmt =('%8s & ' + ' %55s & ' + '%20s &' + '%6s &' + ' %4d (%d) & ' + ' %5.2f & ' + ' %5.2f & ' + '%+5.2f \\\\')
elif printfmt == 'html':
tablehead = \
"<TR>\n <TH> Site code </TH> \
......@@ -136,7 +136,7 @@ def summarize_obs(analysisdir, printfmt='html'):
<TD>%s</TD>\n \
</TR>\n"""
elif printfmt == 'scr':
print 'Code Site NObs flagged R Inn X2'
print('Code Site NObs flagged R Inn X2')
fmt = '%8s ' + ' %55s %s %s' + ' %4d ' + ' %4d ' + ' %5.2f ' + ' %5.2f'
table = []
......
......@@ -156,6 +156,7 @@ class Optimizer(object):
members = statevector.ensemble_members[n]
for m, mem in enumerate(members):
members[m].param_values[:] = self.X_prime[n * self.nparams:(n + 1) * self.nparams, m] + self.x[n * self.nparams:(n + 1) * self.nparams]
members[m].param_values[:][members[m].param_values[:]<0.] = 0.
logging.debug('Returning optimized data to the StateVector, setting "StateVector.isOptimized = True" ')
......@@ -325,7 +326,6 @@ class Optimizer(object):
res = self.obs[n] - self.Hx[n]
if self.may_reject[n]:
# print('may reject')
threshold = self.rejection_threshold * np.sqrt(self.R[n])
if np.abs(res) > threshold:
logging.debug('Rejecting observation (%s,%i) because residual (%f) exceeds threshold (%f)' % (self.sitecode[n], self.obs_ids[n], res, threshold))
......@@ -336,6 +336,7 @@ class Optimizer(object):
PHt = 1. / (self.nmembers - 1) * np.dot(self.X_prime, self.HX_prime[n, :])
self.HPHR[n] = 1. / (self.nmembers - 1) * (self.HX_prime[n, :] * self.HX_prime[n, :]).sum() + self.R[n]
print('chi2', res*res/self.HPHR[n])
self.KG[:] = PHt / self.HPHR[n]
......@@ -366,8 +367,19 @@ class Optimizer(object):
fac = 1.0 / (self.nmembers - 1) * (self.HX_prime[n, :] * self.HX_prime[m, :]).sum() / self.HPHR[n]
self.Hx[m] = self.Hx[m] + fac * res
self.HX_prime[m, :] = self.HX_prime[m, :] - alpha * fac * self.HX_prime[n, :]
# by pavle
# chi2 = []
# chi2_HPHR = np.zeros((self.nobs,), float)
# for n in range(self.nobs):
# res = self.obs[n] - self.Hx[n]
#chi2_HPHR[n] = 1. / (self.nmembers-1) * (self.HX_prime[n, :] * self.HX_prime[n, :]).sum() + self.R
# print('chi2',res*res/self.HPHR[n])
# #print('chi2',res*res/chi2_HPHR[n])
# chi2.append(res*res/self.HPHR[n])
# chi2 = np.asarray(chi2)
# print('chi2 mean',np.mean(chi2))
# end by pavle
def bulk_minimum_least_squares(self):
""" Make minimum least squares solution by solving matrix equations"""
......
......@@ -76,6 +76,8 @@ class CO2StateVector(StateVector):
(0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 0.16, 0.64, 0.16), \
(0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.16, 0.16, 0.16, 0.64) ])
# L = 300 km
L_matrix = np.array([\
(1.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000), \
(0.000, 1.000, 0.715, 0.650, 0.733, 0.614, 0.909, 0.688, 0.544), \
......@@ -158,6 +160,7 @@ class CO2StateVector(StateVector):
for m in range(self.nmembers):
newmember = EnsembleMember(m)
newmember.param_values = EnsembleMembers[m, :].flatten() + meanstate # add the mean to the deviations to hold the full parameter values
# newmember.param_values[newmember.param_values<0.] = 0
self.ensemble_members[n].append(newmember)
f.close()
......
......@@ -156,6 +156,7 @@ class CO2StateVector(StateVector):
for m in range(self.nmembers):
newmember = EnsembleMember(m)
newmember.param_values = EnsembleMembers[m, :].flatten() + meanstate # add the mean to the deviations to hold the full parameter values
newmember.param_values[newmember.param_values<0.] = 0.
self.ensemble_members[n].append(newmember)
f.close()
......
......@@ -20,7 +20,8 @@ from da.tools.general import to_datetime
identifier = 'ObservationOperator'
version = '10'
cdo = Cdo()
#cdo = Cdo()
cdo = Cdo(logging=True, logFile='cdo_commands.log')
################### Begin Class ObservationOperator ###################
class ObservationOperator(object):
......
......@@ -258,7 +258,6 @@ def prepare_state(dacycle, statevector):
saved_sv = os.path.join(dacycle['dir.restart'], 'savestate_%s.nc' % dacycle['da.restart.tstamp'].strftime('%Y%m%d'))
statevector.read_from_file(saved_sv) # by default will read "opt"(imized) variables, and then propagate
statevector.remove_negs(saved_sv)
# Now propagate the ensemble by one cycle to prepare for the current cycle
statevector.propagate(dacycle)
......
......@@ -425,10 +425,10 @@ class StateVector(object):
for m in range(self.nmembers):
newmember = EnsembleMember(m)
nonegs = ensmembers[n, m, :].flatten() + meanstate[n]
nonegs[nonegs<0.] = 0.
newmember.param_values = nonegs
# newmember.param_values = ensmembers[n, m, :].flatten() + meanstate[n] # add the mean to the deviations to hold the full parameter values
# nonegs = ensmembers[n, m, :].flatten() + meanstate[n]
# nonegs[nonegs<0.] = 0.
# newmember.param_values = nonegs
newmember.param_values = ensmembers[n, m, :].flatten() + meanstate[n] # add the mean to the deviations to hold the full parameter values
self.ensemble_members[n].append(newmember)
logging.info('Successfully read the State Vector from file (%s) ' % filename)
......
......@@ -418,10 +418,10 @@ class StateVector(object):
for m in range(self.nmembers):
newmember = EnsembleMember(m)
# nonegs = ensmembers[n, m, :].flatten() + meanstate[n]
# nonegs[nonegs<0.] = 0.
# newmember.param_values = nonegs
newmember.param_values = ensmembers[n, m, :].flatten() + meanstate[n] # add the mean to the deviations to hold the full parameter values
nonegs = ensmembers[n, m, :].flatten() + meanstate[n]
nonegs[nonegs<0.] = 0.
newmember.param_values = nonegs
# newmember.param_values = ensmembers[n, m, :].flatten() + meanstate[n] # add the mean to the deviations to hold the full parameter values
self.ensemble_members[n].append(newmember)
logging.info('Successfully read the State Vector from file (%s) ' % filename)
......
......@@ -13,7 +13,7 @@
!!! Info for the CarbonTracker data assimilation system
datadir : /store/empa/em05/parsenov/obs
datadir : /store/empa/em05/parsenov/ch_obs
! For ObsPack
obspack.input.id : obs
......
Markdown is supported
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