Skip to content
Snippets Groups Projects

Py3 statevector

Merged Peters, Wouter requested to merge peter050/CTDAS:py3-statevector into master
4 files
+ 15
59
Compare changes
  • Side-by-side
  • Inline
Files
4
@@ -34,6 +34,7 @@ your own baseclass StateVector we refer to :ref:`tut_chapter5`.
"""
import os
import sys
import logging
import numpy as np
from datetime import timedelta
@@ -162,46 +163,14 @@ class StateVector(object):
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
# 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.gridmap = np.random.randint(low=1,high=self.nparams+1,size=(180,360,))
self.griddict = {}
for r in range(1, int(nparams) + 1):
sel = np.nonzero(self.gridmap.flat == r)
for r in range(1, self.nparams+1):
sel = np.nonzero(self.gridmap.flat == r)
if len(sel[0]) > 0:
self.griddict[r] = sel
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
self.make_species_mask()
@@ -230,7 +199,7 @@ class StateVector(object):
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 covariancematrix: a matrix to draw random values from
Loading