diff --git a/da/cosmo/covariances.py b/da/cosmo/covariances.py index d6d65fd8b932ef913a6c72e106a2cebb153ab24f..8fefeda7712e55ed38d35730bcd4951145bd2a9a 100755 --- a/da/cosmo/covariances.py +++ b/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 diff --git a/da/cosmo/observationoperator.py b/da/cosmo/observationoperator.py index 112de43c18fee0ba10258605c936b921250bfd1b..16cdad8bd85bad768b041d77c0c382ef8feee8b6 100755 --- a/da/cosmo/observationoperator.py +++ b/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'] diff --git a/da/cosmo/pipeline.py b/da/cosmo/pipeline.py index c37802df9b3f0e5d5b28c48354f1878776c71469..ba86435b4fc5102db912400fd081172552fd6d32 100755 --- a/da/cosmo/pipeline.py +++ b/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 diff --git a/da/cosmo/statevector.py b/da/cosmo/statevector.py index be4b818362a93b0d9fda493158bc8c1f4da91da0..892c76e71104f7453d7c9f95de8452f2ec55f122 100644 --- a/da/cosmo/statevector.py +++ b/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)