.gitignore

+*.pyc
+*/*.pyc*.pyc
+*/*.pyc*.pyc
+*/*.pyc*.pyc
+*/*.pyc*.pyc
+*/*.pyc*.pyc
+*/*.pyc*.pyc
+*/*.pyc*.pyc
+*/*.pyc

da/analysis/expand_molefractions.py

@@ -19,7 +19,6 @@ import shutil
 import logging
 import netCDF4
 import numpy as np
-from string import join
 from datetime import datetime, timedelta
 sys.path.append('../../')
 from da.tools.general import date2num, num2date

da/analysis/tools_transcom.py

@@ -18,7 +18,6 @@ sys.path.append('../../')
 rootdir = os.getcwd().split('da/')[0]
 analysisdir = os.path.join(rootdir, 'da/analysis')
 
-from string import join, split
 from numpy import array, identity, zeros, arange, dot
 import da.tools.io4 as io
 

da/baseclasses/statevector.py

@@ -246,7 +246,7 @@ class StateVector(object):
 
         """    
 
-        if covariancematrix == None: 
+        if covariancematrix is None: 
             covariancematrix = np.identity(self.nparams)
 
         # Make a cholesky decomposition of the covariance matrix
@@ -411,7 +411,7 @@ class StateVector(object):
 
         logging.info('Successfully read the State Vector from file (%s) ' % filename)
 
-    def write_members_to_file(self, lag, outdir,endswith='.nc'):
+    def write_members_to_file(self, lag, outdir,endswith='.nc', obsoperator=None):
         """ 
            :param: lag: Which lag step of the filter to write, must lie in range [1,...,nlag]
            :param: outdir: Directory where to write files

da/carbondioxide/obs.py

@@ -23,9 +23,7 @@ File created on 28 Jul 2010.
 import os
 import sys
 import logging
-#from da.baseclasses.statevector import filename
 import datetime as dtm
-from string import strip
 from numpy import array, logical_and
 
 sys.path.append(os.getcwd())

da/carbondioxide/obspack.py

@@ -24,7 +24,6 @@ import os
 import sys
 import logging
 import datetime as dtm
-#from string import strip
 from numpy import array, logical_and
 
 sys.path.append(os.getcwd())

da/carbondioxide/obspack_geocarbon.py

@@ -25,7 +25,6 @@ import sys
 import logging
         
 import datetime as dtm
-from string import strip
 from numpy import array, logical_and
 sys.path.append(os.getcwd())
 sys.path.append('../../')

da/carbondioxide/obspack_globalviewplus.py

@@ -25,7 +25,6 @@ import sys
 import logging
         
 import datetime as dtm
-#from string import strip
 from numpy import array, logical_and, sqrt
 sys.path.append(os.getcwd())
 sys.path.append('../../')

da/carbondioxide/optimizer.py

-"""CarbonTracker Data Assimilation Shell (CTDAS) Copyright (C) 2017 Wouter Peters. 
+"""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. 
+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. 
+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 
+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
 # optimizer.py
 
 """
-Author : peters 
+Author : peters
 
 Revision History:
 File created on 28 Jul 2010.
@@ -24,6 +24,8 @@ File created on 28 Jul 2010.
 import os
 import sys
 import logging
+import numpy as np
+import da.tools.io4 as io
 sys.path.append(os.getcwd())
 from da.baseclasses.optimizer import Optimizer
 
@@ -41,6 +43,127 @@ class CO2Optimizer(Optimizer):
         All other methods are inherited from the base class Optimizer.
     """
 
+    def create_matrices(self):
+        """ Create Matrix space needed in optimization routine """
+
+        # mean state  [X]
+        self.x = np.zeros((self.nlag * self.nparams,), float)
+        # deviations from mean state  [X']
+        self.X_prime = np.zeros((self.nlag * self.nparams, self.nmembers,), float)
+        # mean state, transported to observation space [ H(X) ]
+        self.Hx = np.zeros((self.nobs,), float)
+        # deviations from mean state, transported to observation space [ H(X') ]
+        self.HX_prime = np.zeros((self.nobs, self.nmembers), float)
+        # observations
+        self.obs = np.zeros((self.nobs,), float)
+        # observation ids
+        self.obs_ids = np.zeros((self.nobs,), float)
+        # covariance of observations
+        # Total covariance of fluxes and obs in units of obs [H P H^t + R]
+        if self.algorithm == 'Serial':
+            self.R = np.zeros((self.nobs,), float)
+            self.HPHR = np.zeros((self.nobs,), float)
+        else:
+            self.R = np.zeros((self.nobs, self.nobs,), float)
+            self.HPHR = np.zeros((self.nobs, self.nobs,), float)
+        # localization of obs
+        self.may_localize = np.zeros(self.nobs, bool)
+        # rejection of obs
+        self.may_reject = np.zeros(self.nobs, bool)
+        # flags of obs
+        self.flags = np.zeros(self.nobs, int)
+        # species type
+        self.species = np.zeros(self.nobs, str)
+        # species type
+        self.sitecode = np.zeros(self.nobs, str)
+        # rejection_threshold
+        self.rejection_threshold = np.zeros(self.nobs, float)
+        # lat and lon
+        self.latitude = np.zeros(self.nobs, float)
+        self.longitude = np.zeros(self.nobs, float)
+        self.date = np.zeros((self.nobs, 6), float)
+
+        # species mask
+        self.speciesmask = {}
+
+        # Kalman Gain matrix
+        #self.KG = np.zeros((self.nlag * self.nparams, self.nobs,), float)
+        self.KG = np.zeros((self.nlag * self.nparams,), float)
+
+
+    def state_to_matrix(self, statevector):
+        allsites     = []  # collect all obs for n=1,..,nlag
+        alllats      = []  # collect all latitudes for n=1,..,nlag
+        alllons      = []  # collect all longitudes for n=1,..,nlag
+        alldates     = []  # collect all date vectors for n=1,..,nlag
+        allobs       = []  # collect all obs for n=1,..,nlag
+        allmdm       = []  # collect all mdm for n=1,..,nlag
+        allids       = []  # collect all model samples for n=1,..,nlag
+        allreject    = []  # collect all model samples for n=1,..,nlag
+        alllocalize  = []  # collect all model samples for n=1,..,nlag
+        allflags     = []  # collect all model samples for n=1,..,nlag
+        allspecies   = []  # collect all model samples for n=1,..,nlag
+        allsimulated = []  # collect all members model samples for n=1,..,nlag
+        allrej_thres = []  # collect all rejection_thresholds, will be the same for all samples of same source
+
+        for n in range(self.nlag):
+
+            samples = statevector.obs_to_assimilate[n]
+            members = statevector.ensemble_members[n]
+            self.x[n * self.nparams:(n + 1) * self.nparams]          = members[0].param_values
+            self.X_prime[n * self.nparams:(n + 1) * self.nparams, :] = np.transpose(np.array([m.param_values for m in members]))
+
+            # Add observation data for all sample objects
+            if samples != None:
+                if type(samples) != list: samples = [samples]
+                for m in range(len(samples)):
+                    sample = samples[m]
+                    logging.debug('Lag %i, sample %i: rejection_threshold = %i, nobs = %i' %(n, m, sample.rejection_threshold, sample.getlength()))
+
+                    allrej_thres.extend([sample.rejection_threshold] * sample.getlength())
+                    allreject.extend(sample.getvalues('may_reject'))
+                    alllocalize.extend(sample.getvalues('may_localize'))
+                    allflags.extend(sample.getvalues('flag'))
+                    allspecies.extend(sample.getvalues('species'))
+                    allobs.extend(sample.getvalues('obs'))
+                    allsites.extend(sample.getvalues('code'))
+                    alllats.extend(sample.getvalues('lat'))
+                    alllons.extend(sample.getvalues('lon'))
+                    alldates.extend([[d.year, d.month, d.day, d.hour, d.minute, d.second] for d in sample.getvalues('xdate')])
+                    allmdm.extend(sample.getvalues('mdm'))
+                    allids.extend(sample.getvalues('id'))
+
+                    simulatedensemble = sample.getvalues('simulated')
+                    for s in range(simulatedensemble.shape[0]):
+                        allsimulated.append(simulatedensemble[s])
+
+        self.rejection_threshold[:] = np.array(allrej_thres)
+
+        self.obs[:]          = np.array(allobs)
+        self.obs_ids[:]      = np.array(allids)
+        self.HX_prime[:, :]  = np.array(allsimulated)
+        self.Hx[:]           = self.HX_prime[:, 0]
+
+        self.may_reject[:]   = np.array(allreject)
+        self.may_localize[:] = np.array(alllocalize)
+        self.flags[:]        = np.array(allflags)
+        self.species[:]      = np.array(allspecies)
+        self.sitecode        = allsites
+        self.latitude[:]     = np.array(alllats)
+        self.longitude[:]    = np.array(alllons)
+        self.date[:,:]       = np.array(alldates)
+
+        self.X_prime  = self.X_prime - self.x[:, np.newaxis] # make into a deviation matrix
+        self.HX_prime = self.HX_prime - self.Hx[:, np.newaxis] # make a deviation matrix
+
+        if self.algorithm == 'Serial':
+            for i, mdm in enumerate(allmdm):
+                self.R[i] = mdm ** 2
+        else:
+            for i, mdm in enumerate(allmdm):
+                self.R[i, i] = mdm ** 2
+
+
     def set_localization(self, loctype='None'):
         """ determine which localization to use """
 
@@ -55,12 +178,12 @@ class CO2Optimizer(Optimizer):
             elif self.nmembers == 150:
                 self.tvalue = 1.97591
             elif self.nmembers == 200:
-                self.tvalue = 1.9719    
-            else: self.tvalue = 0    
+                self.tvalue = 1.9719
+            else: self.tvalue = 0
         else:
             self.localization = False
             self.localizetype = 'None'
-    
+
         logging.info("Current localization option is set to %s" % self.localizetype)
         if self.localization == True:
             if self.tvalue == 0:
@@ -72,9 +195,9 @@ class CO2Optimizer(Optimizer):
         """ localize the Kalman Gain matrix """
         import numpy as np
 
-        if not self.localization: 
+        if not self.localization:
             logging.debug('Not localized observation %i' % self.obs_ids[n])
-            return 
+            return
         if self.localizetype == 'CT2007':
             count_localized = 0
             for r in range(self.nlag * self.nparams):
@@ -92,9 +215,184 @@ class CO2Optimizer(Optimizer):
             self.algorithm = 'Serial'
         else:
             self.algorithm = 'Bulk'
-    
+
         logging.info("Current minimum least squares algorithm is set to %s" % self.algorithm)
 
+
+    def write_diagnostics(self, filename, type):
+        """
+            Open a NetCDF file and write diagnostic output from optimization process:
+
+                - calculated residuals
+                - model-data mismatches
+                - HPH^T
+                - prior ensemble of samples
+                - posterior ensemble of samples
+                - prior ensemble of fluxes
+                - posterior ensemble of fluxes
+
+            The type designation refers to the writing of prior or posterior data and is used in naming the variables"
+        """
+
+        # Open or create file
+        if type == 'prior':
+            f = io.CT_CDF(filename, method='create')
+            logging.debug('Creating new diagnostics file for optimizer (%s)' % filename)
+        elif type == 'optimized':
+            f = io.CT_CDF(filename, method='write')
+            logging.debug('Opening existing diagnostics file for optimizer (%s)' % filename)
+
+        # Add dimensions
+        dimparams  = f.add_params_dim(self.nparams)
+        dimmembers = f.add_members_dim(self.nmembers)
+        dimlag     = f.add_lag_dim(self.nlag, unlimited=False)
+        dimobs     = f.add_obs_dim(self.nobs)
+        dimstate   = f.add_dim('nstate', self.nparams * self.nlag)
+        dim200char = f.add_dim('string_of200chars', 200)
+        dimdate    = f.add_dim('calendar_components', 6)
+
+        # Add data, first the ones that are written both before and after the optimization
+        savedict              = io.std_savedict.copy()
+        savedict['name']      = "statevectormean_%s" % type
+        savedict['long_name'] = "full_statevector_mean_%s" % type
+        savedict['units']     = "unitless"
+        savedict['dims']      = dimstate
+        savedict['values']    = self.x.tolist()
+        savedict['comment']   = 'Full %s state vector mean ' % type
+        f.add_data(savedict)
+
+        savedict              = io.std_savedict.copy()
+        savedict['name']      = "statevectordeviations_%s" % type
+        savedict['long_name'] = "full_statevector_deviations_%s" % type
+        savedict['units']     = "unitless"
+        savedict['dims']      = dimstate + dimmembers
+        savedict['values']    = self.X_prime.tolist()
+        savedict['comment']   = 'Full state vector %s deviations as resulting from the optimizer' % type
+        f.add_data(savedict)
+
+        savedict              = io.std_savedict.copy()
+        savedict['name']      = "modelsamplesmean_%s" % type
+        savedict['long_name'] = "modelsamplesforecastmean_%s" % type
+        savedict['units']     = "mol mol-1"
+        savedict['dims']      = dimobs
+        savedict['values']    = self.Hx.tolist()
+        savedict['comment']   = '%s mean mole fractions based on %s state vector' % (type, type)
+        f.add_data(savedict)
+
+        savedict              = io.std_savedict.copy()
+        savedict['name']      = "modelsamplesdeviations_%s" % type
+        savedict['long_name'] = "modelsamplesforecastdeviations_%s" % type
+        savedict['units']     = "mol mol-1"
+        savedict['dims']      = dimobs + dimmembers
+        savedict['values']    = self.HX_prime.tolist()
+        savedict['comment']   = '%s mole fraction deviations based on %s state vector' % (type, type)
+        f.add_data(savedict)
+
+        # Continue with prior only data
+
+        if type == 'prior':
+
+            savedict              = io.std_savedict.copy()
+            savedict['name']      = "sitecode"
+            savedict['long_name'] = "site code propagated from observation file"
+            savedict['dtype']     = "char"
+            savedict['dims']      = dimobs + dim200char
+            savedict['values']    = self.sitecode
+            savedict['missing_value'] = '!'
+            f.add_data(savedict)
+
+            savedict              = io.std_savedict.copy()
+            savedict['name']      = 'date_components'
+            savedict['dims']      = dimobs + dimdate
+            savedict['units']     = "integer components of UTC date/time"
+            savedict['missing_value'] = -999.9
+            savedict['order']     = "year, month, day, hour, minute, second"
+            savedict['values']    = self.date
+            f.add_data(savedict)
+
+            savedict              = io.std_savedict.copy()
+            savedict['name']      = 'latitude'
+            savedict['dims']      = dimobs
+            savedict['units']     = "degrees_north"
+            savedict['missing_value'] = -999.9
+            savedict['values']    = self.latitude
+            f.add_data(savedict)
+
+            savedict              = io.std_savedict.copy()
+            savedict['name']      = 'longitude'
+            savedict['units']     = "degrees_east"
+            savedict['missing_value'] = -999.9
+            savedict['dims']      = dimobs
+            savedict['values']    = self.longitude
+            f.add_data(savedict)
+
+            savedict              = io.std_savedict.copy()
+            savedict['name']      = "observed"
+            savedict['long_name'] = "observedvalues"
+            savedict['units']     = "mol mol-1"
+            savedict['dims']      = dimobs
+            savedict['values']    = self.obs.tolist()
+            savedict['comment']   = 'Observations used in optimization'
+            f.add_data(savedict)
+
+            savedict              = io.std_savedict.copy()
+            savedict['name']      = "obspack_num"
+            savedict['dtype']     = "int64"
+            savedict['long_name'] = "Unique_observation_number"
+            savedict['units']     = ""
+            savedict['dims']      = dimobs
+            savedict['values']    = self.obs_ids.tolist()
+            savedict['comment']   = 'Unique observation number across the entire ObsPack or satellite retrieval distribution'
+            f.add_data(savedict)
+
+            savedict              = io.std_savedict.copy()
+            savedict['name']      = "modeldatamismatchvariance"
+            savedict['long_name'] = "modeldatamismatch variance"
+            savedict['units']     = "[mol mol-1]^2"
+            if self.algorithm == 'Serial':
+                savedict['dims'] = dimobs
+            else: savedict['dims'] = dimobs + dimobs
+            savedict['values']    = self.R.tolist()
+            savedict['comment']   = 'Variance of mole fractions resulting from model-data mismatch'
+            f.add_data(savedict)
+
+        # Continue with posterior only data
+
+        elif type == 'optimized':
+
+            savedict              = io.std_savedict.copy()
+            savedict['name']      = "totalmolefractionvariance"
+            savedict['long_name'] = "totalmolefractionvariance"
+            savedict['units'] = "[mol mol-1]^2"
+            if self.algorithm == 'Serial':
+                savedict['dims'] = dimobs
+            else: savedict['dims'] = dimobs + dimobs
+            savedict['values']    = self.HPHR.tolist()
+            savedict['comment']   = 'Variance of mole fractions resulting from prior state and model-data mismatch'
+            f.add_data(savedict)
+
+            savedict              = io.std_savedict.copy()
+            savedict['name']      = "flag"
+            savedict['long_name'] = "flag_for_obs_model"
+            savedict['units']     = "None"
+            savedict['dims']      = dimobs
+            savedict['values']    = self.flags.tolist()
+            savedict['comment']   = 'Flag (0/1/2/99) for observation value, 0 means okay, 1 means QC error, 2 means rejected, 99 means not sampled'
+            f.add_data(savedict)
+
+            #savedict = io.std_savedict.copy()
+            #savedict['name'] = "kalmangainmatrix"
+            #savedict['long_name'] = "kalmangainmatrix"
+            #savedict['units'] = "unitless molefraction-1"
+            #savedict['dims'] = dimstate + dimobs
+            #savedict['values'] = self.KG.tolist()
+            #savedict['comment'] = 'Kalman gain matrix of all obs and state vector elements'
+            #dummy                   = f.add_data(savedict)
+
+        f.close()
+        logging.debug('Diagnostics file closed')
+
+
 ################### End Class CO2Optimizer ###################
 
 if __name__ == "__main__":

da/co2gridded/statevectorNHgridded.py

-"""CarbonTracker Data Assimilation Shell (CTDAS) Copyright (C) 2017 Wouter Peters. 
+"""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. 
+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. 
+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 
+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.
@@ -50,12 +50,12 @@ class CO2GriddedStateVector(StateVector):
 
 
     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:
@@ -63,7 +63,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]
@@ -79,7 +79,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:
@@ -87,14 +87,19 @@ class CO2GriddedStateVector(StateVector):
 
             f = io.ct_read(file, 'read')
 
-            if 'pco2' in file or 'cov_ocean' in file: 
-                cov_ocn = f.get_variable('CORMAT')
-                parnr = range(9805,9835)
-                cov = cov_ocn
+            if 'pco2' in file or 'cov_ocean' in file:
+                parnr = list(range(9805,9835))
+                if 'ocn.optimize' in dacycle.dasystem and not dacycle.dasystem['ocn.optimize']:
+                    cov = np.identity(30)*1E-20
+                    logging.debug('Prior ocean covariance matrix = Identity matrix')
+                else:
+                    cov_ocn = f.get_variable('CORMAT')
+                    cov     = cov_ocn
+                    
             elif 'tropic' in file:
                 cov = f.get_variable('covariance')
                 parnr = f.get_variable('parameternumber')
-            else: 
+            else:
                 cov = f.get_variable('covariance')
                 parnr = f.get_variable('parameternumber')
                 cov_sf = 90.0 / np.sqrt(cov.diagonal().sum())  # this scaling factor makes the total variance close to the value of a single ecoregion
@@ -110,12 +115,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]
 
@@ -123,7 +128,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:
@@ -137,7 +142,7 @@ class CO2GriddedStateVector(StateVector):
 
         dims = 1  # start from 1.0 to account for the last parameter that scales Ice+Non-optimized, we have no covariance matrix for this though
 
-        for matrix,parnr,file in covariancematrixlist: 
+        for matrix,parnr,file in covariancematrixlist:
             dims += matrix.shape[0]
 
         if dims != self.nparams:
@@ -157,7 +162,7 @@ class CO2GriddedStateVector(StateVector):
             dof += np.sum(s) ** 2 / sum(s ** 2)
             try:
                 C = np.linalg.cholesky(matrix)
-            except np.linalg.linalg.LinAlgError, err:
+            except np.linalg.linalg.LinAlgError as err:
                 logging.error('Cholesky decomposition has failed ')
                 logging.error('For a matrix of dimensions: %d' % matrix.shape[0])
                 logging.debug(err)
@@ -190,7 +195,7 @@ class CO2GriddedStateVector(StateVector):
         # 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
 
@@ -198,7 +203,7 @@ class CO2GriddedStateVector(StateVector):
 
         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

da/platform/cartesius.py

@@ -132,9 +132,9 @@ class CartesiusPlatform(Platform):
             logging.info("A new task will be started (%s)" % cmd)
             output = subprocess.Popen(cmd, stdout=subprocess.PIPE).communicate()[0]
             logging.info(output)
-            print 'output', output
+            print('output', output)
             jobid = output.split()[-1]             
-            print 'jobid', jobid
+            print('jobid', jobid)
         else:
             cmd = ["sbatch", jobfile]
             logging.info("A new job will be submitted (%s)" % cmd)

da/platform/cartesius_short.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
+# cartesius.py
+
+"""
+Author : peters 
+
+Revision History:
+File created on 06 Sep 2010.
+
+"""
+
+import logging
+import subprocess
+
+from da.baseclasses.platform import Platform
+from da.platform.cartesius import CartesiusPlatform
+
+std_joboptions = {'jobname':'test', 'jobaccount':'co2', 'jobtype':'serial', 'jobshell':'/bin/sh', 'depends':'', 'jobtime':'01:00:00', 'jobinput':'/dev/null', 'jobnodes':'1', 'jobtasks':'', 'modulenetcdf':'netcdf/4.1.2', 'networkMPI':'','jobqueue': 'normal'}
+
+
+class CartesiusPlatformShort(Platform):
+
+    def get_job_template(self, joboptions={}, block=False):
+        """ 
+        Returns the job template for a given computing system, and fill it with options from the dictionary provided as argument.
+        The job template should return the preamble of a job that can be submitted to a queue on your platform, 
+        examples of popular queuing systems are:
+            - SGE
+            - MOAB
+            - XGrid
+            -
+
+        A list of job options can be passed through a dictionary, which are then filled in on the proper line,
+        an example is for instance passing the dictionary {'account':'co2'} which will be placed 
+        after the ``-A`` flag in a ``qsub`` environment.
+
+        An extra option ``block`` has been added that allows the job template to be configured to block the current
+        job until the submitted job in this template has been completed fully.
+        """
+        
+        #template = """## \n"""+ \
+        #           """## This is a set of dummy names, to be replaced by values from the dictionary \n"""+ \
+        #           """## Please make your own platform specific template with your own keys and place it in a subfolder of the da package.\n """+ \
+        #           """## \n"""+ \
+        #           """ \n"""+ \
+        #           """#$ jobname \n"""+ \
+        #           """#$ jobaccount \n"""+ \
+        #           """#$ jobnodes \n"""+ \
+        #           """#$ jobtime \n"""+ \
+        #           """#$ jobshell \n"""+ \
+        #           """\n"""+ \
+        #           """source /usr/bin/sh\n"""+ \
+        #           """module load python\n"""+ \
+        #           """\n"""
+
+       
+        template = """#!/bin/bash \n""" + \
+                   """## \n""" + \
+                   """## This is a set of dummy names, to be replaced by values from the dictionary \n""" + \
+                   """## Please make your own platform specific template with your own keys and place it in a subfolder of the da package.\n """ + \
+                   """## \n""" + \
+                   """#SBATCH -J jobname \n""" + \
+                   """#SBATCH -p short \n""" + \
+                   """#SBATCH -n jobnodes \n""" + \
+                   """#SBATCH -t 01:00:00 \n""" + \
+                   """#SBATCH -o joblog \n""" + \
+                   """module load pre2019\n""" + \
+		   """module load python\n""" + \
+                   """module load nco\n""" + \
+		   """\n"""
+
+        if 'depends' in joboptions:
+            template += """#$ -hold_jid depends \n"""
+
+        # First replace from passed dictionary
+        for k, v in joboptions.items():
+            while k in template:
+                template = template.replace(k, v)
+
+        # Fill remaining values with std_options
+        for k, v in std_joboptions.items():
+            while k in template:
+                template = template.replace(k, v)
+
+        return template
+ 
+
+
+if __name__ == "__main__":
+    pass

da/rc/ctdas_column_oco2.rc

+!!! Info for the CarbonTracker data assimilation system
+ 
+datadir             : /projects/0/ctdas/input/ctdas_2012
+
+! XCO2 samples
+obs.column.input.dir          : ${datadir}/xco2/oco2_l2_lite_fp_9r_10secs/
+obs.column.ncfile             : oco2_LtCO2_<YYYYMMDD>_10sec.nc
+obs.column.selection.variables: land_water_indicator, operation_mode
+obs.column.selection.criteria : == 0, <2
+obs.column.extra.variables    : pressure_weight
+obs.column.mdm.weight         : 1
+obs.column.rc	              : ${datadir}/xco2/oco2_weights.rc
+obs.column.mdm.calculation    : parametrization
+
+! For ObsPack
+obspack.input.id   : obspack_co2_111_MERGE_GLOBALVIEWplus_v4.2_NRT_v4.4_2019-05-16
+obspack.input.dir  : ${datadir}/obspacks/${obspack.input.id}
+obs.sites.rc       : /projects/0/ctdas/liesbeth/ctdas_xco2/sites_weights_validation_noaircraft.rc
+obspack.mdm.weight : 1
+obs.scaling.factor : 1e6
+
+ocn.covariance  : ${datadir}/oceans/oif/cov_ocean.2000.01.nc 
+deltaco2.prefix : oif_p3_era40.dpco2
+
+bio.cov.dir     : ${datadir}/covariances/gridded_NH/
+bio.cov.prefix  : cov_ecoregion
+
+regtype         : gridded_oif30
+nparameters     : 9835 
+random.seed     : 4385
+regionsfile     : ${datadir}/covariances/gridded_NH/griddedNHparameters.nc
+!random.seed.init: ${datadir}/randomseedinit.pickle
+
+! Include a naming scheme for the variables
+
+#include NamingScheme.wp_Mar2011.rc 
+
+! Info on the sites file used
+
+ocn.optimize       : True

da/tools/initexit.py

@@ -75,7 +75,6 @@ import copy
 import getopt
 import pickle
 import numpy as np
-#from string import join
 
 import da.tools.rc as rc
 from da.tools.general import create_dirs, to_datetime, advance_time