From c991da3c52a03c1d3736e1377d54382fc06ff73e Mon Sep 17 00:00:00 2001
From: Wouter Peters <wouter.peters@wur.nl>
Date: Thu, 13 Jun 2013 13:09:04 +0000
Subject: [PATCH] sf6 code to be updated
---
da/sf6/__init__.py | 0
da/sf6/statevector.py | 192 ++++++++++++++++++++++++++++++++++++++++++
2 files changed, 192 insertions(+)
create mode 100644 da/sf6/__init__.py
create mode 100755 da/sf6/statevector.py
diff --git a/da/sf6/__init__.py b/da/sf6/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/da/sf6/statevector.py b/da/sf6/statevector.py
new file mode 100755
index 00000000..ddebf759
--- /dev/null
+++ b/da/sf6/statevector.py
@@ -0,0 +1,192 @@
+#!/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
+from da.baseclasses.statevector import StateVector
+import numpy as np
+import da.tools.io4 as io
+
+identifier = 'SF6 Statevector '
+version = '0.0'
+
+################### Begin Class SF6StateVector ###################
+
+class SF6StateVector(StateVector):
+ """ This is a StateVector object for CarbonTracker. It has a private method to make new ensemble members """
+
+ def get_covariance(self, date):
+ """ 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]
+ """
+
+ import da.tools.io4 as io
+ try:
+ import matplotlib.pyplot as plt
+ except:
+ pass
+
+ # Get the needed matrices from the specified covariance files
+
+ fullcov = np.zeros((self.nparams, self.nparams), float)
+
+ fullcov[:,:] = 0.2**2
+
+ return fullcov
+
+ def Initialize(self):
+ """
+ Initialize the object by specifying the dimensions.
+ There are two major requirements for each statvector that you want to build:
+
+ (1) is that the statevector can map itself onto a regular grid
+ (2) is that the statevector can map itself (mean+covariance) onto TransCom regions
+
+ An example is given below.
+ """
+
+ self.nlag = int(self.DaCycle['time.nlag'])
+ self.nmembers = int(self.DaCycle['da.optimizer.nmembers'])
+ self.nparams = int(self.DaCycle.DaSystem['nparameters'])
+ self.nobs = 0
+
+ self.isOptimized = False
+ self.ObsToAssimmilate = () # empty containter to hold observations to assimilate later on
+
+ # These list objects hold the data for each time step of lag in the system. Note that the ensembles for each time step consist
+ # of lists of EnsembleMember objects, we define member 0 as the mean of the distribution and n=1,...,nmembers as the spread.
+
+ self.EnsembleMembers = range(self.nlag)
+
+ for n in range(self.nlag):
+ self.EnsembleMembers[n] = []
+
+
+ # This specifies the file to read with the gridded mask at 1x1 degrees. Each gridbox holds a number that specifies the parametermember
+ # that maps onto it. From this map, a dictionary is created that allows a reverse look-up so that we can map parameters to a grid.
+
+ mapfile = os.path.join(self.DaCycle.DaSystem['regionsfile'])
+ ncf = io.CT_Read(mapfile, 'read')
+ self.tcmap = ncf.GetVariable('transcom_regions')
+ ncf.close()
+
+ self.gridmap = np.ones((180,360),'float')
+
+ logging.debug("A TransCom map on 1x1 degree was read from file %s" % self.DaCycle.DaSystem['regionsfile'])
+ logging.debug("A parameter map on 1x1 degree was created")
+
+ # Create a dictionary for state <-> gridded map conversions
+
+ nparams = self.gridmap.max()
+ self.griddict = {}
+ for r in range(1, nparams + 1):
+ sel = (self.gridmap.flat == r).nonzero()
+ if len(sel[0]) > 0:
+ self.griddict[r] = sel
+
+ logging.debug("A dictionary to map grids to states and vice versa was created")
+
+ # Create a matrix for state <-> TransCom conversions
+
+ self.tcmatrix = np.zeros((self.nparams, 23), 'float')
+
+ logging.debug("A matrix to map states to TransCom regions and vice versa was created")
+
+ # Create a mask for species/unknowns
+
+ self.make_species_mask()
+
+ def make_species_mask(self):
+
+ """
+
+ This method creates a dictionary with as key the name of a tracer, and as values an array of 0.0/1.0 values
+ specifying which StateVector elements are constrained by this tracer. This mask can be used in
+ the optimization to ensure that certain types of osbervations only update certain unknowns.
+
+ An example would be that the tracer '14CO2' can be allowed to only map onto fossil fuel emissions in the state
+
+ The form of the mask is:
+
+ {'co2': np.ones(self.nparams), 'co2c14', np.zeros(self.nparams) }
+
+ so that 'co2' maps onto all parameters, and 'co2c14' on none at all. These arrays are used in the Class
+ optimizer when state updates are actually performed
+
+ """
+ self.speciesdict = {'sf6': np.ones(self.nparams)}
+ logging.debug("A species mask was created, only the following species are recognized in this system:")
+ for k in self.speciesdict.keys():
+ logging.debug(" -> %s" % k)
+
+ def propagate(self):
+ """
+ :rtype: None
+
+ Propagate the parameter values in the StateVector to the next cycle. This means a shift by one cycle
+ step for all states that will
+ be optimized once more, and the creation of a new ensemble for the time step that just
+ comes in for the first time (step=nlag).
+ In the future, this routine can incorporate a formal propagation of the statevector.
+
+ """
+ from datetime import timedelta
+
+ ensemble_deviations = np.array([p.ParameterValues for p in self.EnsembleMembers[0] ])
+ if ensemble_deviations.std(axis=0) < 0.05 :
+ for m in self.EnsembleMembers[0]:
+ m = 3.0 * m # inflate deviations by a factor of 3.0
+
+ logging.info('The state vector covariance was inflated to 1-sigma of %5.2f ' % (ensemble_deviations.std(axis=0)*3.0))
+
+ logging.info('The state vector remains the same in the SF6 run')
+ logging.info('The state vector has been propagated by one cycle')
+
+
+################### End Class SF6StateVector ###################
+
+
+if __name__ == "__main__":
+ from da.tools.initexit import StartLogger
+ from da.tools.pipeline import start_job
+
+ sys.path.append(os.getcwd())
+
+ opts = ['-v']
+ args = {'rc':'da.rc', 'logfile':'da_initexit.log', 'jobrcfilename':'test.rc'}
+
+ StartLogger()
+
+ DaCycle = start_job(opts, args)
+
+ DaCycle.Initialize()
+
+ StateVector = CtStateVector()
+
+ StateVector.Initialize()
+
+ for n in range(dims[0]):
+ cov = StateVector.get_covariance()
+ dummy = StateVector.make_new_ensemble(n + 1, cov)
+
+ StateVector.propagate()
+
+ savedir = DaCycle['dir.output']
+ filename = os.path.join(savedir, 'savestate.nc')
+
+ dummy = StateVector.WriteToFile()
+
+ StateVector.ReadFromFile(filename)
+
--
GitLab