Commit 2216a271 authored by Woude, Auke van der's avatar Woude, Auke van der
Browse files

fix coordinates of SiB data

parent 9ba799c0
......@@ -137,11 +137,11 @@ def get_spatial_emissions(dacycle, i_member, species, datepoint):
f.close()
return emis
def average2d(arr, new_shape):
""" Function to average the paint-by-colour NEE to the domain size"""
shape = (new_shape[0], arr.shape[0] // new_shape[0],
new_shape[1], arr.shape[1] // new_shape[1])
return arr.reshape(shape).mean(-1).mean(1)
#def average2d(arr, new_shape):
# """ Function to average the paint-by-colour NEE to the domain size"""
# shape = (new_shape[0], arr.shape[0] // new_shape[0],
# new_shape[1], arr.shape[1] // new_shape[1])
# return arr.reshape(shape).mean(-1).mean(1)
def get_time_index_nc(time=None, startdate=None):
"""Function that gets the time index from the flux files
......@@ -307,15 +307,21 @@ class STILTObservationOperator(ObservationOperator):
lon_ind_end = st.lonindex(self.lon_end)
lat_ind_start = st.latindex(self.lat_start)
lat_ind_end = st.latindex(self.lat_end)
# The SiB data is read in correctly, whilst the remainder
# Of the data is read in upside down. THerefore, we flip
# The SiB data also upside down.
nee_2d = np.flip(st.convert_2d(nee, latsib, lonsib), axis=2)
self.nee = nee_2d[:, :,
nee_2d = nee_2d[:, :,
lat_ind_end: lat_ind_start,
lon_ind_start: lon_ind_end].astype(np.float32)
nee_2d = np.flip(nee_2d, axis=2)
self.nee = nee_2d
lu_pref = np.flip(st.convert_2d(lu_pref, latsib, lonsib), axis=1)
self.lu_pref = lu_pref[:,
lat_ind_end: lat_ind_start,
lon_ind_start: lon_ind_end]
self.lu_pref = np.flip(lu_pref[:,
lat_ind_end: lat_ind_start,
lon_ind_start: lon_ind_end],
axis=1)
def prepare_biosphere(self):
"""Function that prepares the biosphere fluxes on a map
......@@ -358,7 +364,7 @@ class STILTObservationOperator(ObservationOperator):
with rasterio.open(self.pftfile) as dataset:
pftdata = dataset.read(out_shape=new_shape,
resampling=Resampling.mode)
pftdata = pftdata.reshape(pftdata.shape[1:])
pftdata = np.flipud(pftdata.reshape(pftdata.shape[1:]))
logging.getLogger().setLevel(logging.DEBUG)
scaling_factors = np.array(np.array(new_shape) / np.array(self.lu_pref[0].shape), dtype=int)
......@@ -366,11 +372,11 @@ class STILTObservationOperator(ObservationOperator):
logging.debug('Starting paint-by-number')
times = len(self.nee)
timeslist = list(range(times))
lus = np.kron(self.lu_pref, np.ones((scaling_factors), dtype=self.lu_pref.dtype))
lus_prefs = np.kron(self.lu_pref, np.ones((scaling_factors), dtype=self.lu_pref.dtype))
pool = Pool(times)
# Create the function that calculates the conentration
func = partial(self.paint_by_number, self.nee, self.param_values, scaling_factors, lus,
pftdata, self.emismodel.find_in_state, self.domain_shape)
func = partial(self.paint_by_number, self.nee, self.param_values, scaling_factors, lus_prefs,
pftdata, self.emismodel.find_in_state, self.average2d, self.domain_shape)
scaled_nees = pool.map(func, timeslist)
pool.close()
pool.join()
......@@ -378,7 +384,15 @@ class STILTObservationOperator(ObservationOperator):
logging.debug('Paint-by-number done, biosphere prepared')
@staticmethod
def paint_by_number(nee, all_param_values, scaling_factors, lus, pftdata, find_in_state, domain_shape, time_ind):
def average2d(arr, new_shape):
""" Function to average the paint-by-colour NEE to the domain size"""
shape = (new_shape[0], arr.shape[0] // new_shape[0],
new_shape[1], arr.shape[1] // new_shape[1])
return arr.reshape(shape).mean(-1).mean(1)
@staticmethod
def paint_by_number(nee, all_param_values, scaling_factors, lus_prefs, pftdata,
find_in_state, average2d, domain_shape, time_ind):
nmembers = len(all_param_values)
scaled_nees_time = np.zeros((nmembers, *domain_shape), dtype=np.float32)
all_lutypes = np.unique(pftdata.reshape(-1))
......@@ -386,7 +400,7 @@ class STILTObservationOperator(ObservationOperator):
nee = np.kron(nee_time, np.ones((scaling_factors), dtype=np.float32))
for lutype_ind, lutype in enumerate(all_lutypes):
mask = np.where(pftdata==lutype, 1, 0)
nee_lut = np.where(lus==lutype, nee, 0)
nee_lut = np.where(lus_prefs==lutype, nee, 0)
for mem, values in enumerate(all_param_values):
param_values = all_param_values[mem]
index_in_state = find_in_state(param_values, 'BIO', str(lutype), return_index=True)
......
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