9 regions, correlated with L

da/cosmo/covariances.py

@@ -27,6 +27,7 @@ sys.path.append(os.getcwd())
 
 import logging
 import numpy as np
+#from da.cosmo.statevector_read_from_output import StateVector, EnsembleMember
 from da.cosmo.statevector import StateVector, EnsembleMember
 
 import da.tools.io4 as io
@@ -46,17 +47,59 @@ class CO2StateVector(StateVector):
 
                                     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
+                                    2. Boardleaf Decidous, Temperate
+                                    3. Boardleaf Decidous, Boreal
+                                    4. Boardleaf Decidous Shrub, Temperate
+                                    5. Boardleaf Decidous Shrub, Boreal
                                     6. C3 Arctic Grass
                                     7. C3 non-Arctic Grass
                                     8. C4 Grass
                                     9. Crop
-                                   10. None
+                                   #10. None (removed)
 
 """
+
+        fullcov = np.zeros(shape=(90,90))
+
+        partcov = np.array([ \
+        (0.64, 0.36, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
+        (0.36, 0.64, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
+        (0.16, 0.16, 0.64, 0.36, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
+        (0.16, 0.16, 0.36, 0.64, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
+        (0.16, 0.16, 0.16, 0.16, 0.64, 0.36, 0.04, 0.04, 0.04, 0.01), \
+        (0.16, 0.16, 0.16, 0.16, 0.36, 0.64, 0.04, 0.04, 0.04, 0.01), \
+        (0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.64, 0.16, 0.16, 0.16), \
+        (0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 0.64, 0.16, 0.16), \
+        (0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 0.16, 0.64, 0.16), \
+        (0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.16, 0.16, 0.16, 0.64) ])
+
+        L_matrix = np.array([\
+        (1.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000), \
+        (0.000, 1.000, 0.715, 0.650, 0.733, 0.614, 0.909, 0.688, 0.544), \
+        (0.000, 0.715, 1.000, 0.862, 0.751, 0.638, 0.695, 0.501, 0.575), \
+        (0.000, 0.650, 0.862, 1.000, 0.776, 0.696, 0.651, 0.472, 0.639), \
+        (0.000, 0.733, 0.751, 0.776, 1.000, 0.827, 0.774, 0.587, 0.732), \
+        (0.000, 0.614, 0.638, 0.696, 0.827, 1.000, 0.660, 0.537, 0.885), \
+        (0.000, 0.909, 0.695, 0.651, 0.774, 0.660, 1.000, 0.721, 0.586), \
+        (0.000, 0.688, 0.501, 0.472, 0.587, 0.537, 0.721, 1.000, 0.489), \
+        (0.000, 0.544, 0.575, 0.639, 0.732, 0.885, 0.586, 0.489, 1.000) ])
+
+        for i in range(9):
+            for j in range(9):
+                fullcov[i*10:(i+1)*10,j*10:(j+1)*10] = partcov * L_matrix[i,j]
+
+#        correl = np.array([ \
+ #       (0.8, 0.6, 0.4, 0.4, 0.4, 0.4, 0.2, 0.2, 0.2, 0.1), \
+  #      (0.6, 0.8, 0.4, 0.4, 0.4, 0.4, 0.2, 0.2, 0.2, 0.1), \
+   #     (0.4, 0.4, 0.8, 0.6, 0.4, 0.4, 0.2, 0.2, 0.2, 0.1), \
+    #    (0.4, 0.4, 0.6, 0.8, 0.4, 0.4, 0.2, 0.2, 0.2, 0.1), \
+     #   (0.4, 0.4, 0.4, 0.4, 0.8, 0.6, 0.2, 0.2, 0.2, 0.1), \
+      #  (0.4, 0.4, 0.4, 0.4, 0.6, 0.8, 0.2, 0.2, 0.2, 0.1), \
+       # (0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.8, 0.4, 0.4, 0.4), \
+        #(0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.8, 0.4, 0.4), \
+#        (0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.4, 0.8, 0.4), \
+ #       (0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.4, 0.4, 0.4, 0.8) ])
+
 #        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), \
@@ -71,17 +114,6 @@ class CO2StateVector(StateVector):
 #  #      (0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.e-10) ])
   #      (0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.00) ])
 
-        fullcov = np.array([ \
-        (1.00, 0.36, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
-        (0.36, 1.00, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
-        (0.16, 0.16, 1.00, 0.36, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
-        (0.16, 0.16, 0.36, 1.00, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
-        (0.16, 0.16, 0.16, 0.16, 1.00, 0.36, 0.04, 0.04, 0.04, 0.01), \
-        (0.16, 0.16, 0.16, 0.16, 0.36, 1.00, 0.04, 0.04, 0.04, 0.01), \
-        (0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 1.00, 0.16, 0.16, 0.16), \
-        (0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 1.00, 0.16, 0.16), \
-        (0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 0.16, 1.00, 0.16), \
-        (0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.16, 0.16, 0.16, 1.00) ])
 
         return fullcov
 

da/cosmo/observationoperator.py

@@ -56,7 +56,7 @@ class ObservationOperator(object):
         self.simulated_file = os.path.join(self.outputdir, 'samples_simulated.%s.nc' % self.dacycle['time.sample.stamp'])
         self.forecast_nmembers = int(self.dacycle['da.optimizer.nmembers'])
 
-    def run(self,lag,dacycle,statevector):
+    def run(self,lag,dacycle,statevector,advance=False):
         members = statevector.ensemble_members[lag]
         self.forecast_nmembers = int(self.dacycle['da.optimizer.nmembers'])
         #self.nparams = int(dacycle.dasystem['nparameters'])
@@ -170,36 +170,37 @@ class ObservationOperator(object):
         logging.info('COSMO done!')
         os.chdir(dacycle['dir.da_run'])
 
-        args = [
-            (dacycle, starth+168*lag, endh+168*lag-1, n)
-            for n in range(0,self.forecast_nmembers)
-        ]
-
-        with Pool(self.forecast_nmembers) as pool:
-            pool.starmap(self.extract_model_data, args)
-
-        for i in range(0,self.forecast_nmembers):
-            idx = str(i).zfill(3)
- #           cosmo_file = os.path.join('/store/empa/em05/parsenov/cosmo_data/OK_DONT_TOUCH/model_'+idx+'_%s.nc' % dacycle['time.sample.stamp'])    # last run with non-frac
-            cosmo_file = os.path.join('/store/empa/em05/parsenov/cosmo_data/model_'+idx+'_%s.nc' % dacycle['time.sample.stamp'])
-            ifile = Dataset(cosmo_file, mode='r')
-            model_data[i,:] = (np.squeeze(ifile.variables['CO2'][:])*29./44.01)*1E6   # in ppm
-            ifile.close()
-
-        for j,data in enumerate(zip(ids,obs,mdm)):
-            f.variables['obs_num'][j] = data[0]
-            f.variables['flask'][j,:] = model_data[:,j]
-        f.close()
+        if not advance:
+            args = [
+                (dacycle, starth+168*lag, endh+168*lag-1, n)
+                for n in range(0,self.forecast_nmembers)
+            ]
+
+            with Pool(self.forecast_nmembers) as pool:
+                pool.starmap(self.extract_model_data, args)
+
+            for i in range(0,self.forecast_nmembers):
+                idx = str(i).zfill(3)
+     #           cosmo_file = os.path.join('/store/empa/em05/parsenov/cosmo_data/OK_DONT_TOUCH/model_'+idx+'_%s.nc' % dacycle['time.sample.stamp'])    # last run with non-frac
+                cosmo_file = os.path.join('/store/empa/em05/parsenov/cosmo_data/model_'+idx+'_%s.nc' % dacycle['time.sample.stamp'])
+                ifile = Dataset(cosmo_file, mode='r')
+                model_data[i,:] = (np.squeeze(ifile.variables['CO2'][:])*29./44.01)*1E6   # in ppm
+                ifile.close()
+
+            for j,data in enumerate(zip(ids,obs,mdm)):
+                f.variables['obs_num'][j] = data[0]
+                f.variables['flask'][j,:] = model_data[:,j]
+            f.close()
 #### WARNING ACHTUNG PAZNJA POZOR VNEMANIE data[2] is model data mismatch (=1000) by default in tools/io4.py!!! pavle
 
 
         logging.info('ObservationOperator finished successfully, output file written (%s)' % self.simulated_file)
 
-    def run_forecast_model(self, lag, dacycle, statevector):
+    def run_forecast_model(self, lag, dacycle, statevector, advance):
         self.prepare_run()
-        self.run(lag, dacycle, statevector)
+        self.run(lag, dacycle, statevector, advance)
 
-    def extract_model_data(self,dacycle,hstart,hstop,ensnum):
+    def extract_model_data(self, dacycle, hstart, hstop, ensnum):
 
         self.dacycle = dacycle
         time_stamp = dacycle['time.sample.stamp']

da/cosmo/pipeline.py

@@ -336,7 +336,7 @@ def sample_step(dacycle, samples, statevector, obsoperator, lag, advance=False):
 
     # Run the observation operator
 
-    obsoperator.run_forecast_model(lag, dacycle, statevector)
+    obsoperator.run_forecast_model(lag, dacycle, statevector, advance)
 
 
     # Read forecast model samples that were written to NetCDF files by each member. Add them to the exisiting

da/cosmo/statevector.py

@@ -296,13 +296,14 @@ class StateVector(object):
             rands_bg = np.random.uniform(low=-0.05, high=0.05, size=1)
 
             newmember = EnsembleMember(member)
-            rands = rands.reshape([9,20])
-            randsC = np.zeros(shape=(9,20))
+#            rands = rands.reshape([9,20])
+#            randsC = np.zeros(shape=(9,20))
 
-            for r in range(0,9):
-                randsC[r,:] = np.hstack((np.dot(C, rands[r,0:10]),np.dot(C, rands[r,10:20])))
+#            for r in range(0,9):
+#                randsC[r,:] = np.hstack((np.dot(C, rands[r,0:10]),np.dot(C, rands[r,10:20])))
+
+            randsC = np.hstack((np.dot(C, rands[0:90]),np.dot(C, rands[90:180])))
             newmember.param_values = (np.hstack((randsC.flatten(),rands_bg)) + newmean).ravel()    #THIS ONE
-            #newmember.param_values = (np.hstack((np.dot(C, rands[0:11]),np.dot(C, rands[11:22]), rands_bg)) + newmean).ravel()    #THIS ONE
             newmember.param_values[newmember.param_values<0.] = 0.
             newmember.param_values[newmember.param_values>2.] = 2.
             self.ensemble_members[lag].append(newmember)