Commit 1dbbaae1 authored by weihe's avatar weihe
Browse files

revised additive scheme

parent 9ee4af7f
"""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 #!/usr/bin/env python
# ct_statevector_tools.py # ct_statevector_tools.py
""" """
Author : peters Author : peters
Revision History: Revision History:
File created on 28 Jul 2010. File created on 28 Jul 2010.
...@@ -40,12 +28,12 @@ class CO2GriddedStateVector(StateVector): ...@@ -40,12 +28,12 @@ class CO2GriddedStateVector(StateVector):
""" This is a StateVector object for CarbonTracker. It has a private method to make new ensemble members """ """ This is a StateVector object for CarbonTracker. It has a private method to make new ensemble members """
def get_covariance(self, date, dacycle): 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. 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] The argument is thus referring to the lagged state vector as [1,2,3,4,5,..., nlag]
""" """
try: try:
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
except: except:
...@@ -53,7 +41,7 @@ class CO2GriddedStateVector(StateVector): ...@@ -53,7 +41,7 @@ class CO2GriddedStateVector(StateVector):
# Get the needed matrices from the specified covariance files # 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.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] 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): ...@@ -69,7 +57,7 @@ class CO2GriddedStateVector(StateVector):
covariancematrixlist = [] covariancematrixlist = []
for file in cov_files: for file in cov_files:
if not os.path.exists(file): 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) logging.error(msg)
raise IOError, msg raise IOError, msg
else: else:
...@@ -77,26 +65,25 @@ class CO2GriddedStateVector(StateVector): ...@@ -77,26 +65,25 @@ class CO2GriddedStateVector(StateVector):
f = io.ct_read(file, 'read') f = io.ct_read(file, 'read')
if 'pco2' in file: if 'pco2' in file:
cov_ocn = f.get_variable('CORMAT') cov_ocn = f.get_variable('CORMAT')
cov = cov_ocn cov = cov_ocn
else: else:
cov = f.get_variable('covariance') 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 = 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_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
cov = cov * cov_sf 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() f.close()
covariancematrixlist.append(cov) covariancematrixlist.append(cov1)
# Boundary conditions covariance # 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) covariancematrixlist.append(cov)
logging.debug("Succesfully closed files after retrieving prior covariance matrices") logging.debug("Succesfully closed files after retrieving prior covariance matrices")
...@@ -105,12 +92,12 @@ class CO2GriddedStateVector(StateVector): ...@@ -105,12 +92,12 @@ class CO2GriddedStateVector(StateVector):
return covariancematrixlist return covariancematrixlist
def make_new_ensemble(self, lag, covariancematrixlist=[None]): def make_new_ensemble(self, lag, covariancematrixlist=[None]):
""" """
:param lag: an integer indicating the time step in the lag order :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 :param covariancematrix: a list of matrices specifying the covariance distribution to draw from
:rtype: None :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. 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] The argument is thus referring to the lagged state vector as [1,2,3,4,5,..., nlag]
...@@ -118,7 +105,7 @@ class CO2GriddedStateVector(StateVector): ...@@ -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 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. the StateVector nparams becomes very large.
""" """
try: try:
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
except: except:
...@@ -134,7 +121,7 @@ class CO2GriddedStateVector(StateVector): ...@@ -134,7 +121,7 @@ class CO2GriddedStateVector(StateVector):
dims = 0 dims = 0
for matrix in covariancematrixlist: for matrix in covariancematrixlist:
dims += matrix.shape[0] dims += matrix.shape[0]
if dims != self.nparams: if dims != self.nparams:
...@@ -168,50 +155,52 @@ class CO2GriddedStateVector(StateVector): ...@@ -168,50 +155,52 @@ class CO2GriddedStateVector(StateVector):
npoints = matrix.shape[0] npoints = matrix.shape[0]
istop = istop + npoints 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): for member in range(1, self.nmembers):
rands = np.random.randn(npoints) rands = np.random.randn(npoints)
deviations = np.dot(C, rands) deviations = np.dot(C, rands)
dev_matrix[istart:istop, member - 1] = deviations dev_matrix[istart:istop, member - 1] = deviations
#dev_matrix[istop, member - 1] = 1.e-10 * np.random.randn() #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): for member in range(1, self.nmembers):
rands = np.random.randn(npoints) rands = np.random.randn(npoints)
deviations = np.dot(C, rands) 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 dev_matrix[istart:istop, member - 1] = deviations
istart = istart + npoints 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.debug('Successfully constructed a deviation matrix from covariance structure')
logging.info('Appr. degrees of freedom in full covariance matrix is %s' % (int(dof))) 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 # 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 = 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.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 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: if lag == self.nlag - 1 and self.nlag >= 3:
new_mean += self.ensemble_members[lag - 1][0].param_values + \ 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 new_mean = new_mean / 3.0
# Create the first ensemble member with a deviation of 0.0 and add to list # Create the first ensemble member with a deviation of 0.0 and add to list
...@@ -226,7 +215,6 @@ class CO2GriddedStateVector(StateVector): ...@@ -226,7 +215,6 @@ class CO2GriddedStateVector(StateVector):
new_member = EnsembleMember(member) new_member = EnsembleMember(member)
new_member.param_values = dev_matrix[:, member - 1] + new_mean new_member.param_values = dev_matrix[:, member - 1] + new_mean
self.ensemble_members[lag].append(new_member) self.ensemble_members[lag].append(new_member)
logging.debug('%d new ensemble members were added to the state vector # %d' % (self.nmembers, (lag + 1))) logging.debug('%d new ensemble members were added to the state vector # %d' % (self.nmembers, (lag + 1)))
...@@ -259,8 +247,10 @@ class CO2GriddedStateVector(StateVector): ...@@ -259,8 +247,10 @@ class CO2GriddedStateVector(StateVector):
for mem in members: for mem in members:
filename = os.path.join(outdir, 'parameters.%03d%s' % (mem.membernumber, endswith)) filename = os.path.join(outdir, 'parameters.%03d%s' % (mem.membernumber, endswith))
ncf = io.CT_CDF(filename, method='create') ncf = io.CT_CDF(filename, method='create')
dimparams = ncf.add_params_dim(self.nparams-1) #dimparams = ncf.add_params_dim(self.nparams-1)
dimparams_bc = ncf.add_dim('nparameters_bc',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() dimgrid = ncf.add_latlon_dim()
data = mem.param_values data = mem.param_values
...@@ -270,7 +260,8 @@ class CO2GriddedStateVector(StateVector): ...@@ -270,7 +260,8 @@ class CO2GriddedStateVector(StateVector):
savedict['long_name'] = "parameter_values_for_member_%d" % mem.membernumber savedict['long_name'] = "parameter_values_for_member_%d" % mem.membernumber
savedict['units'] = "unitless" savedict['units'] = "unitless"
savedict['dims'] = dimparams 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 savedict['comment'] = 'These are parameter values to use for member %d' % mem.membernumber
ncf.add_data(savedict) ncf.add_data(savedict)
...@@ -290,7 +281,8 @@ class CO2GriddedStateVector(StateVector): ...@@ -290,7 +281,8 @@ class CO2GriddedStateVector(StateVector):
savedict['long_name'] = "parameter_values_for_member_%d" % mem.membernumber savedict['long_name'] = "parameter_values_for_member_%d" % mem.membernumber
savedict['units'] = "unitless" savedict['units'] = "unitless"
savedict['dims'] = dimparams_bc 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 savedict['comment'] = 'These are parameter values to use for member %d' % mem.membernumber
ncf.add_data(savedict) ncf.add_data(savedict)
......
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