made it possible to optimize per site. Made parameters named

da/ffdas/obs.py

@@ -121,7 +121,9 @@ class RINGOObservations(Observations):
 
         ncf = io.ct_read(filename, method='read')
         ids = ncf.get_variable('obs_num')
-        simulated = ncf.get_variable('model')
+        simulated = np.array(ncf.get_variable('model'))
+        self.ids = ids
+
         ncf.close()
         logging.info("Successfully read data from model sample file (%s)" % filename)
 

da/ffdas/observationoperator.py

@@ -41,6 +41,9 @@ except: # The memoization package is not always included. Import the memoization
     cached = cached()
 import xlrd
 
+
+from da.ffdas import category_info
+categories = category_info.categories
 # global constants, which will be used in the following classes
 identifier = 'WRF-STILT'
 version = '1.0'
@@ -91,7 +94,7 @@ class STILTObservationOperator(ObservationOperator):
         logging.info('Observation Operator initialized: %s (%s)' % (self.ID, self.version))
 
         self.load_rc(filename)   # load the specified rc-file
-        
+        self.filename = filename 
         if dacycle != None:
             self.dacycle = dacycle
         else:
@@ -116,8 +119,9 @@ class STILTObservationOperator(ObservationOperator):
         self.bgswitch = int(dacycle.dasystem['obs.bgswitch'])
         self.bgfile = dacycle.dasystem['obs.background']
         self.btime = int(dacycle.dasystem['run.backtime'])
-        self.paramfile = dacycle.dasystem['emis.paramfiles'] + '_NED.csv' # One for checking
-        
+        self.categories = categories
+        self.paramdict = rc.read('/home/awoude/ffdas/RINGO/Data/paramdict.rc')
+
     def load_rc(self, name):
         """This method loads a STILT rc-file with settings for this simulation
         Args:
@@ -198,24 +202,9 @@ class STILTObservationOperator(ObservationOperator):
                         self.lon.append(lon)
                         self.hgt.append(height)
                     
-        #find out which model is to be used for which sector
-        infile = os.path.join(self.obsdir, self.paramfile)
-        f = open(infile, 'r', encoding='utf-8-sig')
-        lines = f.readlines()
-        f.close()
-        self.ops_sector = []
-        self.ops_id = []
         self.temporal_var = []
-        for line in lines:
-            if line[0] == '#': continue
-            else:
-                model = int(line.split(',')[1])
-                temporal_var = line.split(',')[5]
-                if id != 0:
-                    self.ops_sector.append(ct)
-                    self.ops_id.append(id)
-                self.temporal_var.append(temporal_var)
-
+        for k, v in self.categories.items():
+            self.temporal_var.append(v['temporal'])
         #set time control variables for this cycle
         if do_pseudo==0:
             self.timestartkey = dtm.datetime(self.dacycle['time.sample.start'].year,
@@ -235,19 +224,12 @@ class STILTObservationOperator(ObservationOperator):
                                               self.dacycle['time.finish'].month,
                                               self.dacycle['time.finish'].day,
                                               self.dacycle['time.finish'].hour, 0)
-#        dt=dtm.timedelta(0,3600)
-#        dumdate=self.timestartkey
-#        self.datelist=[]
-#        while dumdate<(self.timefinishkey): # self.timefinishkey : dt
-#            self.datelist.append(dumdate)
-#            dumdate=dumdate+dt
-
         self.get_c14_time_profile()
 
     def get_samples(self, samples):
         self.samples = samples.datalist
         datelist = [sample.xdate for sample in self.samples]
-        self.datelist = sorted(set(datelist))
+        self.datelist = sorted(datelist)
         logging.info('Added samples to the observation operator')
 
     def get_c14_time_profile(self):
@@ -345,8 +327,7 @@ class STILTObservationOperator(ObservationOperator):
         
         return background_conc
     
-    @cached
-    def get_background_orig(self, i_species, datepoint):
+    def get_background_orig(self, i_species):
         """Function that finds the background concentration, non-time dependent and hard-coded.
         Input:
             i_species: int: the index of the species for which the background should be found
@@ -458,6 +439,7 @@ class STILTObservationOperator(ObservationOperator):
         if total: return gpp_increase + ter_increase
         else: return gpp_increase, ter_increase
 
+    @cached
     def get_temporal_profiles(self, datepoint, station_name,  i_member):
         """ Function that reads in the temporal profiles for the current timestep
         Input:
@@ -471,8 +453,9 @@ class STILTObservationOperator(ObservationOperator):
         temporal_prior = os.path.join(self.model_settings['datadir'],'prior_temporal_{0}_{1:03}.nc'.format(station_name, i_member))
         temporal_prior = io.ct_read(temporal_prior, 'read')
         temp_profiles = []
-        for cat in range(self.nrcat):
-            temporal_variation = temporal_prior.get_variable(self.temporal_var[cat])[time_indices]
+        for category, values in self.categories.items():
+            temporal_var_name = values['temporal']
+            temporal_variation = temporal_prior.get_variable(temporal_var_name)[time_indices]
             temp_profiles.append(temporal_variation)
         #temporal_prior.close()
         temp_profiles = np.array(temp_profiles)
@@ -480,7 +463,7 @@ class STILTObservationOperator(ObservationOperator):
         return temp_profiles.T # Transpose to be similar to spatial data in dimensions
 
     @cached
-    def get_spatial_emissions(self, i_member, i_species):
+    def get_spatial_emissions(self, i_member, station,  i_species):
         """ Function that gets the spatial emissions
         Input:
             i_member: int: the index of the member for which the simulation is run
@@ -488,7 +471,8 @@ class STILTObservationOperator(ObservationOperator):
         Returns:
             np.ndarray (3d): the spatial emissions per category, lat and lon"""
         #read spatially distributed emissions calculated with the emission model
-        emisfile = os.path.join(self.model_settings['datadir'],'prior_spatial_%03d.nc'%i_member) #format: species[SNAP,lon,lat]
+        emisfile = os.path.join(self.model_settings['datadir'],'prior_spatial_{0}_{1:03d}.nc'.format(station, i_member)) #format: species[SNAP,lon,lat]
+        self.emisfile = emisfile
         f = io.ct_read(emisfile,method='read')
         emis=f.get_variable(self.spname[i_species])
         f.close()
@@ -507,7 +491,7 @@ class STILTObservationOperator(ObservationOperator):
             float: the concentration increase at the respective location due to the emissions from the respective species"""
         # get the date:
         temp_profiles = self.get_temporal_profiles(datepoint, site,  i_member)
-        spatial_emissions = self.get_spatial_emissions(i_member, i_species)
+        spatial_emissions = self.get_spatial_emissions(i_member, site, i_species)
         #find correct background concentrationa
         footprint = self.get_foot(site, datepoint)
 
@@ -527,12 +511,14 @@ class STILTObservationOperator(ObservationOperator):
             self.prepare_run(do_pseudo,adv)
             self.run(dacycle,do_pseudo)
             self.save_obs()
+        self.run_breakdown()
 
     def run(self,dacycle,do_pseudo):
         """Function that loops over all members, locations, species and datepoints and calculates the concentration increase at that point (in time)
         Adds the data as attribute to this object."""
 
-        totser=np.array(int(self.dacycle['da.optimizer.nmembers'])*[(self.nrloc*self.nrspc*len(self.datelist))*[0.]])
+        totser=np.array(int(self.dacycle['da.optimizer.nmembers'])*[len(self.datelist)*[0.]])
+        self.totser = totser
         for i_member in range(int(self.dacycle['da.optimizer.nmembers'])):
             logging.debug('Calculating concentration for member {}'.format(i_member))
             conc_STILT=[]
@@ -542,19 +528,37 @@ class STILTObservationOperator(ObservationOperator):
                     species = sample.species; i_species = self.spname.index(species.upper())
                     datepoint = sample.xdate
                     site = sample.code.split('/')[-1].split('_')[1]
-                    logging.debug('working on time {}'.format(datepoint))
-                    background_concentration = self.get_background_orig(i_species, datepoint)
+                    logging.debug('working on site {} time {}'.format(site, datepoint))
+                    background_concentration = self.get_background_orig(i_species)
                     conc_increase_stilt = self.run_STILT(datepoint, i_member, site, i_species)
                     if species.upper() == 'CO2':
                         gpp, ter = self.get_biosphere_concentration(site, datepoint)
                         logging.debug('{0:.3} for FF and {1:.3} for bio for {2}'.format(conc_increase_stilt, gpp+ter, self.spname[i_species]))
-                    else: gpp, ter = 0, 0
+                    else:
+                        gpp, ter = 0, 0
+                        logging.debug('{0:.3} for CO '.format(conc_increase_stilt))
                     conc_STILT.append(conc_increase_stilt + gpp + ter + background_concentration)
             totser[i_member] = np.array(conc_STILT)
             logging.info('Concentration calculated for member {}'.format(i_member))
         self.mod = np.array(totser).T
         logging.debug('Finished doing simulations per back trajectory time')
 
+    def run_breakdown(self):
+        """Function that clears cache for all functions that are cached
+        so that next time, run will make new forecasts
+
+        commented out functions are cached, but not optimised and do 
+        therefore not need an empty cache."""
+        self.run_STILT.cache_clear()
+        self.get_spatial_emissions.cache_clear()
+        self.get_nc_variable.cache_clear()
+        # self.get_biosphere_concentration.cache_clear()
+        # self.get_background.cache_clear()
+        # self.self.get_background_orig.cache_clear()
+        self.get_c14_concentration.cache_clear()
+        # self.get_foot.cache.clear()
+        self.get_temporal_profiles.cache_clear()
+
     def save_data(self):
         """ Write the data that is needed for a restart or recovery of the Observation Operator to the save directory """
 

da/ffdas/statevector.py

@@ -202,27 +202,27 @@ class CO2StateVector(StateVector):
             dt=selectdate.strftime('%Y%m%d')
             file=os.path.join(dacycle['dir.da_run'], 'output', selectdate.strftime('%Y%m%d'),'optimizer.%s.nc' % selectdate.strftime('%Y%m%d'))
             f = io.ct_read(file, 'read')
-            prmval = f.get_variable('statevectormean_optimized')[:]
+            prmval = f.get_variable('statevectormean_optimized')[:][:self.nparams]
             devs = f.get_variable('statevectordeviations_optimized')[:]
             f.close()
-            covariancematrix = (np.dot(devs,devs.T)/(devs.shape[1]-1))
+
             import pickle
+            covariancematrix = (np.dot(devs,devs.T)/(devs.shape[1]-1))
+            covariancematrix = covariancematrix[:len(covariancematrix)//self.nlag, :len(covariancematrix)//self.nlag]
             with open('covmatrix.pkl', 'wb') as tofile:
                 pickle.dump(covariancematrix, tofile)
+            diag_indices = np.diag_indices_from(covariancematrix)
+            covariancematrix[diag_indices] = min(covariancematrix.max()*2, 1.)
+            with open('covmatrix_diag.pkl', 'wb') as tofile:
+                pickle.dump(covariancematrix, tofile)
     
-        # Check dimensions of covariance matrix list, must add up to nparams
         covariancematrix = np.array(covariancematrix)
         self.covariancematrix = covariancematrix
-#        self.devs = devs
-#        diag_covmatrix = np.zeros_like(covariancematrix)
-#        diag_covmatrix[np.diag_indices_from(covariancematrix)] = np.diagonal(covariancematrix)
-#        covariancematrix = diag_covmatrix
-#        logging.warning('Using only the diagonal values of the covariancematrix')
 
         dims = covariancematrix.shape[0]
-        if dims != self.nparams:
-            logging.error("The total dimension of the covariance matrices passed (%d) does not add up to the prescribed nparams (%d), exiting..." % (dims, self.nparams))
-            raise ValueError
+        #if dims != self.nparams:
+        #    logging.error("The total dimension of the covariance matrices passed (%d) does not add up to the prescribed nparams (%d), exiting..." % (dims, self.nparams))
+        #    raise ValueError
 
         # Make a cholesky decomposition of the covariance matrix
 
@@ -239,7 +239,7 @@ class CO2StateVector(StateVector):
         logging.info('Appr. degrees of freedom in covariance matrix is %s' % (int(dof)))
 
         # Create mean values 
-
+        self.prmval = prmval
         newmean = np.ones(self.nparams, float) * prmval # 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

da/ffdas/stilt-ops_urbanall.rc

@@ -5,16 +5,15 @@ datadir         : /home/awoude/ffdas/RINGO/Data
 ! list of all observation sites
 obs.input.id   : obsfiles.csv
 ! number of observation sites included; number of species included and to be used in inversion
-obs.input.nr   : 2
-obs.spec.nr    : 1
-obs.spec.name  : co2
+obs.input.nr   : 6
+obs.spec.nr    : 2
+obs.spec.name  : CO2,CO
 ! number of emission categories defined in the emission model
 obs.cat.nr     : 14
 ! For Rdam obs
 obs.sites.rc        : ${datadir}/sites_weights.rc
-obs.sites.times     : 12; 16
 ! number of parameters
-nparameters     : 50
+nparameters     : 22
 ! set fixed seed for random number generator, or use 0 if you want to use any random seed
 random.seed     : 4385
 !file with prior estimate of scaling factors (statevector) and covariances
@@ -50,4 +49,4 @@ run.backtime            : 24
 ! 1: no propagation, start each cycle with the same prior parameter values and covariance matrix
 ! 2: propagation of optimized parameter values, but not of the covariance matrix
 ! 3: propagation of both optimized parameter values and covariance matrix
-run.propscheme          : 2
+run.propscheme          : 3