Commit 915eb8de authored by Peters, Wouter's avatar Peters, Wouter
Browse files

Merge branch 'py3-statevector' into 'master'

Py3 statevector

See merge request ctdas/CTDAS!15
parents 42c6f431 b6d0761b
...@@ -34,6 +34,7 @@ your own baseclass StateVector we refer to :ref:`tut_chapter5`. ...@@ -34,6 +34,7 @@ your own baseclass StateVector we refer to :ref:`tut_chapter5`.
""" """
import os import os
import sys
import logging import logging
import numpy as np import numpy as np
from datetime import timedelta from datetime import timedelta
...@@ -162,46 +163,14 @@ class StateVector(object): ...@@ -162,46 +163,14 @@ class StateVector(object):
self.ensemble_members[n] = [] self.ensemble_members[n] = []
# This specifies the file to read with the gridded mask at 1x1 degrees. Each gridbox holds a number that specifies the parametermember self.gridmap = np.random.randint(low=1,high=self.nparams+1,size=(180,360,))
# 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(dacycle.dasystem['regionsfile'])
ncf = io.ct_read(mapfile, 'read')
self.gridmap = ncf.get_variable('regions')
self.tcmap = ncf.get_variable('transcom_regions')
ncf.close()
logging.debug("A TransCom map on 1x1 degree was read from file %s" % dacycle.dasystem['regionsfile'])
logging.debug("A parameter map on 1x1 degree was read from file %s" % dacycle.dasystem['regionsfile'])
# Create a dictionary for state <-> gridded map conversions
nparams = self.gridmap.max()
self.griddict = {} self.griddict = {}
for r in range(1, int(nparams) + 1): for r in range(1, self.nparams+1):
sel = np.nonzero(self.gridmap.flat == r) sel = np.nonzero(self.gridmap.flat == r)
if len(sel[0]) > 0: if len(sel[0]) > 0:
self.griddict[r] = sel self.griddict[r] = sel
logging.debug("A dictionary to map grids to states and vice versa was created") 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')
for r in range(1, self.nparams + 1):
sel = np.nonzero(self.gridmap.flat == r)
if len(sel[0]) < 1:
continue
else:
n_tc = set(self.tcmap.flatten().take(sel[0]))
if len(n_tc) > 1:
logging.error("Parameter %d seems to map to multiple TransCom regions (%s), I do not know how to handle this" % (r, n_tc))
raise ValueError
self.tcmatrix[r - 1, n_tc.pop() - 1] = 1.0
logging.debug("A matrix to map states to TransCom regions and vice versa was created")
# Create a mask for species/unknowns # Create a mask for species/unknowns
self.make_species_mask() self.make_species_mask()
...@@ -230,7 +199,7 @@ class StateVector(object): ...@@ -230,7 +199,7 @@ class StateVector(object):
logging.debug(" -> %s" % k) logging.debug(" -> %s" % k)
def make_new_ensemble(self, lag, covariancematrix=None): def make_new_ensemble(self, lag, covariancematrix=[None]):
""" """
:param lag: an integer indicating the time step in the lag order :param lag: an integer indicating the time step in the lag order
:param covariancematrix: a matrix to draw random values from :param covariancematrix: a matrix to draw random values from
......
...@@ -11,15 +11,7 @@ ...@@ -11,15 +11,7 @@
! 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/>. ! program. If not, see <http://www.gnu.org/licenses/>.
datadir : /Storage/CO2/carbontracker/input/ctdas_2016/ nparameters : 100
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 random.seed : 4385
regionsfile : ${datadir}/covariances/gridded_NH/griddedNHparameters.nc
obs.sites.rc : Random obs.sites.rc : Random
...@@ -332,25 +332,20 @@ class CycleControl(dict): ...@@ -332,25 +332,20 @@ class CycleControl(dict):
strippedname = os.path.split(self['jobrcfilename'])[-1] strippedname = os.path.split(self['jobrcfilename'])[-1]
self['jobrcfilename'] = os.path.join(self['dir.exec'], strippedname) self['jobrcfilename'] = os.path.join(self['dir.exec'], strippedname)
# shutil.copy(os.path.join(self.dasystem['regionsfile']),os.path.join(self['dir.exec'],'da','analysis','copied_regions.nc'))
logging.info('Copied regions file to the analysis directory: %s'%os.path.join(self.dasystem['regionsfile']))
if 'extendedregionsfile' in self.dasystem: if 'extendedregionsfile' in self.dasystem:
# shutil.copy(os.path.join(self.dasystem['extendedregionsfile']),os.path.join(self['dir.exec'],'da','analysis','copied_regions_extended.nc')) # shutil.copy(os.path.join(self.dasystem['extendedregionsfile']),os.path.join(self['dir.exec'],'da','analysis','copied_regions_extended.nc'))
logging.info('Copied extended regions file to the analysis directory: %s'%os.path.join(self.dasystem['extendedregionsfile'])) logging.info('Copied extended regions file to the analysis directory: %s'%os.path.join(self.dasystem['extendedregionsfile']))
else: #for filename in glob.glob(os.path.join(self['dir.exec'],'da','analysis','*.pickle')):
# shutil.copy(os.path.join(self['dir.exec'],'da','analysis','olson_extended.nc'),os.path.join(self['dir.exec'],'da','analysis','copied_regions_extended.nc')) # logging.info('Deleting pickle file %s to make sure the correct regions are used'%os.path.split(filename)[1])
logging.info('Copied extended regions within the analysis directory: %s'%os.path.join(self['dir.exec'],'da','analysis','olson_extended.nc')) # os.remove(filename)
for filename in glob.glob(os.path.join(self['dir.exec'],'da','analysis','*.pickle')): #for filename in glob.glob(os.path.join(self['dir.exec'],'*.pickle')):
logging.info('Deleting pickle file %s to make sure the correct regions are used'%os.path.split(filename)[1]) # logging.info('Deleting pickle file %s to make sure the correct regions are used'%os.path.split(filename)[1])
os.remove(filename) # os.remove(filename)
for filename in glob.glob(os.path.join(self['dir.exec'],'*.pickle')):
logging.info('Deleting pickle file %s to make sure the correct regions are used'%os.path.split(filename)[1])
os.remove(filename)
if 'random.seed.init' in self.dasystem: if 'random.seed.init' in self.dasystem:
self.read_random_seed(True) self.read_random_seed(True)
self.parse_times() self.parse_times()
#self.write_rc(self['jobrcfilename'])
def setup_file_structure(self): def setup_file_structure(self):
""" """
......
...@@ -29,7 +29,7 @@ from da.tools.initexit import start_logger, validate_opts_args, parse_options, C ...@@ -29,7 +29,7 @@ from da.tools.initexit import start_logger, validate_opts_args, parse_options, C
from da.tools.pipeline import ensemble_smoother_pipeline, header, footer from da.tools.pipeline import ensemble_smoother_pipeline, header, footer
from da.platform.capegrim import CapeGrimPlatform from da.platform.capegrim import CapeGrimPlatform
from da.baseclasses.dasystem import DaSystem from da.baseclasses.dasystem import DaSystem
from da.co2gridded.statevector import CO2GriddedStateVector from da.baseclasses.statevector import StateVector
from da.baseclasses.obs import Observations from da.baseclasses.obs import Observations
from da.baseclasses.optimizer import Optimizer from da.baseclasses.optimizer import Optimizer
from da.baseclasses.observationoperator import RandomizerObservationOperator from da.baseclasses.observationoperator import RandomizerObservationOperator
...@@ -54,7 +54,7 @@ platform = CapeGrimPlatform() ...@@ -54,7 +54,7 @@ platform = CapeGrimPlatform()
dasystem = DaSystem(dacycle['da.system.rc']) dasystem = DaSystem(dacycle['da.system.rc'])
obsoperator = RandomizerObservationOperator(dacycle['da.obsoperator.rc']) obsoperator = RandomizerObservationOperator(dacycle['da.obsoperator.rc'])
samples = Observations() samples = Observations()
statevector = CO2GriddedStateVector() statevector = StateVector()
optimizer = Optimizer() optimizer = Optimizer()
########################################################################################## ##########################################################################################
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment