diff --git a/da/stilt/statevector.py b/da/stilt/statevector.py
index d24902e7aaae9874037e48532778a9d197674de5..65198b9cd3061c349f512afaf1df518ef96d57a8 100755
--- a/da/stilt/statevector.py
+++ b/da/stilt/statevector.py
@@ -1,20 +1,8 @@
-"""CarbonTracker Data Assimilation Shell (CTDAS) Copyright (C) 2017 Wouter Peters.
-Users are recommended to contact the developers (wouter.peters@wur.nl) to receive
-updates of the code. See also: http://www.carbontracker.eu.
-
-This program is free software: you can redistribute it and/or modify it under the
-terms of the GNU General Public License as published by the Free Software Foundation,
-version 3. This program is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along with this
-program. If not, see <http://www.gnu.org/licenses/>."""
#!/usr/bin/env python
# ct_statevector_tools.py
"""
-Author : peters
+Author : peters
Revision History:
File created on 28 Jul 2010.
@@ -40,12 +28,12 @@ class CO2GriddedStateVector(StateVector):
""" This is a StateVector object for CarbonTracker. It has a private method to make new ensemble members """
def get_covariance(self, date, dacycle):
- """ Make a new ensemble from specified matrices, the attribute lag refers to the position in the state vector.
+ """ Make a new ensemble from specified matrices, the attribute lag refers to the position in the state vector.
Note that lag=1 means an index of 0 in python, hence the notation lag-1 in the indexing below.
The argument is thus referring to the lagged state vector as [1,2,3,4,5,..., nlag]
- """
+ """
+
-
try:
import matplotlib.pyplot as plt
except:
@@ -53,7 +41,7 @@ class CO2GriddedStateVector(StateVector):
# Get the needed matrices from the specified covariance files
- #file_ocn_cov = dacycle.dasystem['ocn.covariance']
+ #file_ocn_cov = dacycle.dasystem['ocn.covariance']
cov_files = os.listdir(dacycle.dasystem['bio.cov.dir'])
cov_files = [os.path.join(dacycle.dasystem['bio.cov.dir'], f) for f in cov_files if dacycle.dasystem['bio.cov.prefix'] in f]
@@ -69,7 +57,7 @@ class CO2GriddedStateVector(StateVector):
covariancematrixlist = []
for file in cov_files:
if not os.path.exists(file):
- msg = "Cannot find the specified file %s" % file
+ msg = "Cannot find the specified file %s" % file
logging.error(msg)
raise IOError, msg
else:
@@ -77,26 +65,25 @@ class CO2GriddedStateVector(StateVector):
f = io.ct_read(file, 'read')
- if 'pco2' in file:
+ if 'pco2' in file:
cov_ocn = f.get_variable('CORMAT')
cov = cov_ocn
- else:
+ else:
cov = f.get_variable('covariance')
#cov_sf = 10.0/np.sqrt(cov.diagonal().sum()) # this scaling factor makes the total variance close to the value of a single ecoregion
- cov_sf = 20.0 / np.sqrt(cov.diagonal().sum()) # this scaling factor makes the total variance close to the value of a single ecoregion
- cov = cov * cov_sf
+ cov_sf = 360. / np.sqrt(cov.diagonal().sum()) # this scaling factor makes the total variance close to the value of a single ecoregion #I use 360 to boost up the P matrix uncertainty
+ cov1 = cov * cov_sf * (1.e-6)**2 # here you assume that your P matrix has units of mol m-2 s-1 squared.
f.close()
- covariancematrixlist.append(cov)
+ covariancematrixlist.append(cov1)
# Boundary conditions covariance
- cov = np.array([[2]]) #np.ones((1,1),)
+ cov = np.array([[2*2]])
+ covariancematrixlist.append(cov)
+ covariancematrixlist.append(cov)
+ covariancematrixlist.append(cov)
covariancematrixlist.append(cov)
-
-
-
-
logging.debug("Succesfully closed files after retrieving prior covariance matrices")
@@ -105,12 +92,12 @@ class CO2GriddedStateVector(StateVector):
return covariancematrixlist
def make_new_ensemble(self, lag, covariancematrixlist=[None]):
- """
+ """
:param lag: an integer indicating the time step in the lag order
:param covariancematrix: a list of matrices specifying the covariance distribution to draw from
:rtype: None
-
- Make a new ensemble, the attribute lag refers to the position in the state vector.
+
+ Make a new ensemble, the attribute lag refers to the position in the state vector.
Note that lag=1 means an index of 0 in python, hence the notation lag-1 in the indexing below.
The argument is thus referring to the lagged state vector as [1,2,3,4,5,..., nlag]
@@ -118,7 +105,7 @@ class CO2GriddedStateVector(StateVector):
used to draw ensemblemembers from. Each draw is done on a matrix from the list, to make the computational burden smaller when
the StateVector nparams becomes very large.
- """
+ """
try:
import matplotlib.pyplot as plt
except:
@@ -134,7 +121,7 @@ class CO2GriddedStateVector(StateVector):
dims = 0
- for matrix in covariancematrixlist:
+ for matrix in covariancematrixlist:
dims += matrix.shape[0]
if dims != self.nparams:
@@ -168,50 +155,52 @@ class CO2GriddedStateVector(StateVector):
npoints = matrix.shape[0]
istop = istop + npoints
- logging.info('i %s,%s'%(i,matrix ))
-
-
- if i < len(covariancematrixlist)-1:
+ #if i < len(covariancematrixlist)-1:
+ if i < len(covariancematrixlist)-4:
for member in range(1, self.nmembers):
rands = np.random.randn(npoints)
deviations = np.dot(C, rands)
dev_matrix[istart:istop, member - 1] = deviations
#dev_matrix[istop, member - 1] = 1.e-10 * np.random.randn()
- logging.info('deviation %s,%s,%s'%(deviations[0],istart,istop))
- #cov2 = np.dot(dev_matrix[istart:istop,:],np.transpose(dev_matrix[istart:istop,:])) / (self.nmembers-1)
- #print matrix.sum(),cov2.sum(),abs(matrix.diagonal()-cov2.diagonal()).max(), matrix.shape,cov2.shape
- #istart = istart + npoints
- #logging.info('i %s,%s'%(i,len(covariancematrixlist)-1), )
- if i == len(covariancematrixlist)-1:
+ #if i == len(covariancematrixlist)-1:
+ if i >= len(covariancematrixlist)-4 and i <= len(covariancematrixlist)-1:
for member in range(1, self.nmembers):
rands = np.random.randn(npoints)
deviations = np.dot(C, rands)
- logging.info('deviation boundary conditions %s,%s,%s'%(deviations,istart,istop))
+ #logging.info('deviation boundary conditions %s,%s,%s'%(deviations,istart,istop))
dev_matrix[istart:istop, member - 1] = deviations
istart = istart + npoints
+ #for i in range(self.nparams):
+ #logging.info('i=%s,deviation=%s'%(i,dev_matrix[i,:]))
+
+
logging.debug('Successfully constructed a deviation matrix from covariance structure')
logging.info('Appr. degrees of freedom in full covariance matrix is %s' % (int(dof)))
# Now fill the ensemble members with the deviations we have just created
- # Create mean values
+ # Create mean values
- new_mean = np.ones(self.nparams, float) # standard value for a new time step is 1.0
- new_mean[-1] = np.zeros((1),float) # standard value for boundary conditions for new time step is 0.0
+ #new_mean = np.ones(self.nparams, float) # standard value for a new time step is 1.0
+ new_mean = np.zeros(self.nparams, float)
+ new_mean[-4] = np.zeros((1),float) # standard value for boundary conditions for new time step is 0.0
+ new_mean[-3] = np.zeros((1),float)
+ new_mean[-2] = np.zeros((1),float)
+ new_mean[-1] = np.zeros((1),float)
# If this is not the start of the filter, average previous two optimized steps into the mix
if lag == self.nlag - 1 and self.nlag >= 3:
new_mean += self.ensemble_members[lag - 1][0].param_values + \
- self.ensemble_members[lag - 2][0].param_values
+ self.ensemble_members[lag - 2][0].param_values
new_mean = new_mean / 3.0
# Create the first ensemble member with a deviation of 0.0 and add to list
@@ -226,7 +215,6 @@ class CO2GriddedStateVector(StateVector):
new_member = EnsembleMember(member)
new_member.param_values = dev_matrix[:, member - 1] + new_mean
self.ensemble_members[lag].append(new_member)
-
logging.debug('%d new ensemble members were added to the state vector # %d' % (self.nmembers, (lag + 1)))
@@ -259,8 +247,10 @@ class CO2GriddedStateVector(StateVector):
for mem in members:
filename = os.path.join(outdir, 'parameters.%03d%s' % (mem.membernumber, endswith))
ncf = io.CT_CDF(filename, method='create')
- dimparams = ncf.add_params_dim(self.nparams-1)
- dimparams_bc = ncf.add_dim('nparameters_bc',1)
+ #dimparams = ncf.add_params_dim(self.nparams-1)
+ dimparams = ncf.add_params_dim(self.nparams-4)
+ #dimparams_bc = ncf.add_dim('nparameters_bc',1)
+ dimparams_bc = ncf.add_dim('nparameters_bc',4)
dimgrid = ncf.add_latlon_dim()
data = mem.param_values
@@ -270,7 +260,8 @@ class CO2GriddedStateVector(StateVector):
savedict['long_name'] = "parameter_values_for_member_%d" % mem.membernumber
savedict['units'] = "unitless"
savedict['dims'] = dimparams
- savedict['values'] = data[:self.nparams-1]
+ #savedict['values'] = data[:self.nparams-1]
+ savedict['values'] = data[:self.nparams-4]
savedict['comment'] = 'These are parameter values to use for member %d' % mem.membernumber
ncf.add_data(savedict)
@@ -290,7 +281,8 @@ class CO2GriddedStateVector(StateVector):
savedict['long_name'] = "parameter_values_for_member_%d" % mem.membernumber
savedict['units'] = "unitless"
savedict['dims'] = dimparams_bc
- savedict['values'] = data[-1]
+ #savedict['values'] = data[-1]
+ savedict['values'] = data[-4:] #data[-4:-1]
savedict['comment'] = 'These are parameter values to use for member %d' % mem.membernumber
ncf.add_data(savedict)