No commit message

No commit message

da/cosmo/statevector.py → da/cosmo/base_statevector.py

@@ -251,8 +251,8 @@ class StateVector(object):
 
         """    
 
-        if covariancematrix == None: 
-            covariancematrix = np.identity(self.nparams)
+#        if covariancematrix == None: 
+ #           covariancematrix = np.identity(self.nparams)
 
         # Make a cholesky decomposition of the covariance matrix
 
@@ -273,7 +273,6 @@ class StateVector(object):
         newmean = np.ones(self.nparams, float) # standard value for a new time step is 1.0
 
         # 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:
             newmean += self.ensemble_members[lag - 1][0].param_values + \
                                            self.ensemble_members[lag - 2][0].param_values 
@@ -291,17 +290,18 @@ class StateVector(object):
             rands = np.random.uniform(low=-1., high=1., size=self.nparams)
  #           rands = np.random.randn(self.nparams)
 #            rands = np.clip(rands,-1,1)
-#            if (member==1):
- #               rands = np.array([1.836451,1.91572,1.719371,1.947495,1.788597,1.109507,0.9401712,2,0.003201536,1.2276,0])-1.
-  #          elif (member==2):
-   #             rands = np.array([2,0.04773442,1.809479,1.136452,2,0.8753765,1.488298,2,1.908293,0.02784704,0])-1.
-    #        elif (member==3):
-     #           rands = np.array([0.9136019,1.386204,2,0.1250395,1.815127,2,0.3002012,1.666945,0,0,0])-1.
+            if (member==1):
+                rands = np.array([1.836451,1.91572,1.719371,1.947495,1.788597,1.109507,0.9401712,2,0.003201536,1.2276,0])-1.
+            elif (member==2):
+                rands = np.array([2,0.04773442,1.809479,1.136452,2,0.8753765,1.488298,2,1.908293,0.02784704,0])-1.
+            elif (member==3):
+                rands = np.array([0.9136019,1.386204,2,0.1250395,1.815127,2,0.3002012,1.666945,0,0,0])-1.
 
             newmember = EnsembleMember(member)
             newmember.param_values = np.dot(C, rands) + newmean
             newmember.param_values[-1] = 0
             self.ensemble_members[lag].append(newmember)
+            newmember.param_values[newmember.param_values<0.] = 0.
 
         logging.debug('%d new ensemble members were added to the state vector # %d' % (self.nmembers, (lag + 1)))
 
@@ -452,7 +452,7 @@ class StateVector(object):
         members = self.ensemble_members[lag]
 
         for mem in members:
-            filename = os.path.join(outdir, 'parameters.%03d%s' % (mem.membernumber, endswith))
+            filename = os.path.join(outdir, 'parameters_'+str(lag)+'.%03d%s' % (mem.membernumber, endswith))
             ncf = io.CT_CDF(filename, method='create')
             dimparams = ncf.add_params_dim(self.nparams)
             dimgrid = ncf.add_latlon_dim()

da/cosmo/covariances.py

+"""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 
+
+Revision History:
+File created on 28 Jul 2010.
+
+"""
+
+import os
+import sys
+sys.path.append(os.getcwd())
+
+import logging
+import numpy as np
+from da.cosmo.base_statevector import StateVector, EnsembleMember
+
+import da.tools.io4 as io
+
+identifier = 'CarbonTracker Statevector '
+version = '0.0'
+
+################### Begin Class CO2StateVector ###################
+
+class CO2StateVector(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. 
+            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]
+
+                                    0. Needleleaf Evergreen, Temperate
+                                    1. Needleleaf Evergreen, Boreal
+                                    2. Boradleaf Decidous, Temperate
+                                    3. Boradleaf Decidous, Boreal
+                                    4. Boradleaf Decidous Shrub, Temperate
+                                    5. Boradleaf Decidous Shrub, Boreal
+                                    6. C3 Arctic Grass
+                                    7. C3 non-Arctic Grass
+                                    8. C4 Grass
+                                    9. Crop
+                                   10. None
+
+"""
+        fullcov = np.array([ \
+        (1.00, 0.36, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
+        (0.36, 1.00, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
+        (0.16, 0.16, 1.00, 0.36, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
+        (0.16, 0.16, 0.36, 1.00, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
+        (0.16, 0.16, 0.16, 0.16, 1.00, 0.36, 0.04, 0.04, 0.04, 0.01, 0.00), \
+        (0.16, 0.16, 0.16, 0.16, 0.36, 1.00, 0.04, 0.04, 0.04, 0.01, 0.00), \
+        (0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 1.00, 0.16, 0.16, 0.16, 0.00), \
+        (0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 1.00, 0.16, 0.16, 0.00), \
+        (0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 0.16, 1.00, 0.16, 0.00), \
+        (0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.16, 0.16, 0.16, 1.00, 0.00), \
+        (0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.e-10) ])
+
+        return fullcov
+
+    def read_from_legacy_file(self, filename, qual='opt'):
+        """ 
+        :param filename: the full filename for the input NetCDF file
+        :param qual: a string indicating whether to read the 'prior' or 'opt'(imized) StateVector from file
+        :rtype: None
+
+        Read the StateVector information from a NetCDF file and put in a StateVector object
+        In principle the input file will have only one four datasets inside 
+        called:
+            * `meanstate_prior`, dimensions [nlag, nparamaters]
+            * `ensemblestate_prior`, dimensions [nlag,nmembers, nparameters]
+            * `meanstate_opt`, dimensions [nlag, nparamaters]
+            * `ensemblestate_opt`, dimensions [nlag,nmembers, nparameters]
+
+        This NetCDF information can be written to file using 
+        :meth:`~da.baseclasses.statevector.StateVector.write_to_file`
+
+        """
+        
+
+        f = io.ct_read(filename, 'read')
+
+        for n in range(self.nlag):
+            if qual == 'opt':
+                meanstate = f.get_variable('statevectormean_opt')
+                EnsembleMembers = f.get_variable('statevectorensemble_opt')
+
+            elif qual == 'prior':
+                meanstate = f.get_variable('statevectormean_prior')
+                EnsembleMembers = f.get_variable('statevectorensemble_prior')
+
+            if not self.ensemble_members[n] == []:
+                self.ensemble_members[n] = []
+                logging.warning('Existing ensemble for lag=%d was removed to make place for newly read data' % (n + 1))
+
+            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
+                self.ensemble_members[n].append(newmember)
+
+        f.close()
+
+        logging.info('Successfully read the State Vector from file (%s) ' % filename)
+    
+    def read_from_file_exceptsam(self, filename, qual='opt'):
+        """ 
+        :param filename: the full filename for the input NetCDF file
+        :param qual: a string indicating whether to read the 'prior' or 'opt'(imized) StateVector from file
+        :rtype: None
+
+        Read the StateVector information from a NetCDF file and put in a StateVector object
+        In principle the input file will have only one four datasets inside 
+        called:
+            * `meanstate_prior`, dimensions [nlag, nparamaters]
+            * `ensemblestate_prior`, dimensions [nlag,nmembers, nparameters]
+            * `meanstate_opt`, dimensions [nlag, nparamaters]
+            * `ensemblestate_opt`, dimensions [nlag,nmembers, nparameters]
+
+        This NetCDF information can be written to file using 
+        :meth:`~da.baseclasses.statevector.StateVector.write_to_file`
+
+        """
+
+        f = io.ct_read(filename, 'read')
+        
+        meanstate = f.get_variable('statevectormean_' + qual)
+#        meanstate[:,39:77] = 1
+        ensmembers = f.get_variable('statevectorensemble_' + qual)
+        f.close()
+
+        for n in range(self.nlag):
+            if not self.ensemble_members[n] == []:
+                self.ensemble_members[n] = []
+                logging.warning('Existing ensemble for lag=%d was removed to make place for newly read data' % (n + 1))
+
+            for m in range(self.nmembers):
+                newmember = EnsembleMember(m)
+                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)
+#        logging.info('State Vector set to 1 for South American regions')
+
+################### End Class CO2StateVector ###################
+
+
+if __name__ == "__main__":
+    pass