import numpy as np
import pandas as pd
from scipy import ndimage
from logging import debug, warning
from .....plugins.transforms.system.fromcontrol.utils.scalemaps \
import map2scale, vmap2vaggreg
from .....utils.dataarrays.reindex import reindex
from .....utils.check.errclass import CifError
[docs]
def get_physical(controlvect, tracer, comp, trcr, **kwargs):
debug(f"Loading physical values for {comp}/{trcr}. This can take a while")
trid = (comp, trcr)
datavect = controlvect.datavect
# Scale xb if prescribed in Yaml
xb_scale = getattr(tracer, "xb_scale", 1.0)
# Filling with defined value if prescribed
xb_value = getattr(tracer, "xb_value", 0.0)
# Raise warning here if merged_dates was not created
if not hasattr(tracer, "merged_dates"):
warning(
f"WARNING! 'merged_dates' was not initialized for {comp} / {trcr}. "
f"This could be due to the tracer not being used in your observation "
f"operator pipe line. Please check your yaml and your pipeline."
)
return np.zeros(tracer.dim), np.zeros(tracer.dim)
merged_dates = tracer.merged_dates
nb_inputs = np.zeros(tracer.dim)
xb = np.zeros(tracer.dim)
std = np.zeros(tracer.dim)
for di in tracer.input_dates:
debug(f"Loading physical values for {di}")
# Skip if empty period
if len(tracer.input_dates[di]) == 0:
continue
# Turn dates to datetime for consistency in read
# TODO: This should be standardized in the future
dates2read = [
[x.to_pydatetime() for x in row]
for row in tracer.input_dates[di].itertuples(index=False, name=None)
]
inputs = tracer.read(
trcr,
tracer.varname,
dates2read,
tracer.input_files[di],
comp_type=comp,
tracer=tracer,
model=controlvect.model,
ddi=di,
**kwargs
)
# Reprojecting inputs and scaling factors to the same merged
# time steps
outdates = pd.to_datetime(
np.unique(np.array(tracer.merged_dates[di]).flatten())
)
if len(outdates) == 1:
outdates = outdates.append(outdates)
# Reindex inputs to merged dates
inputs = reindex(
inputs,
levels={"time": outdates[:-1]},
)
# Corresponding indexes in control vector and inputs
control_indexes = \
pd.Series(range(tracer.ndates),
index=tracer.dates).reindex(
outdates, method="ffill")
# Find indexes in the control vector corresponding to each data
# periods
period_dates = merged_dates[di]
for k, period in period_dates.iterrows():
# Either take the corresponding slice of time,
# or take the exact date
# if the control variable is on a time stamp
try:
dd0, dd1 = period
mask = (outdates[:-1] >= dd0) & (outdates[:-1] < dd1)
# If period is a time stamp, i.e., dd0 = dd1, adapt the mask
if dd0 == dd1:
raise CifError("Test what happens here")
mask = outdates == dd0
except TypeError:
raise CifError("Test what happens here")
dd0 = period
mask = outdates == dd0
control_slice = control_indexes.iloc[k]
# If control_slice is nan,
# it means that the date in the sensitivity
# is not in the control vector, hence skipping
if np.isnan(control_slice):
continue
# Reference matrix with error values
errtype = getattr(tracer, "errtype", "")
err_scale = tracer.err_scale
err_value = getattr(tracer, "err_value", 0)
errdata = err_scale * np.abs(inputs[mask].values) + err_value
if errtype == "max":
errdata = ndimage.maximum_filter(errdata, size=3)
# If errtype is 'avg', prescribes uniform uncertainties
elif errtype == "avg":
errdata[:] = errdata.mean() / 2
# It errtype is "threshold", put all errors below the threshold
# (in absolute value) to 0
elif errtype == "threshold":
errdata = np.where(
errdata <= tracer.err_threshold,
tracer.err_threshold,
errdata
)
# If inputs are below a certain value, put to zero the error
elif errtype == "threshold_input":
errdata = np.where(
np.abs(inputs[mask].values) <= tracer.err_threshold,
0,
errdata
)
# Temporal aggregation
vdata = np.sum(inputs[mask], axis=0).values[np.newaxis]
vcounter = np.sum(np.ones_like(inputs[mask]), axis=0)[np.newaxis]
verrdata = np.sum(errdata, axis=0)[np.newaxis]
# Vertical aggregation depending on vertical stacks
vaggreg = vmap2vaggreg(vdata, tracer, tracer.domain, trid)
vcounteraggreg = vmap2vaggreg(
vcounter, tracer, tracer.domain, trid)
verraggreg = vmap2vaggreg(verrdata, tracer, tracer.domain, trid)
# Horizontal aggregation depending on regions
xstack = map2scale(vaggreg, tracer, tracer.domain,
region_scale_area=True)
xcounter = map2scale(vcounteraggreg, tracer, tracer.domain,
region_scale_area=True)
xerr = map2scale(verraggreg, tracer, tracer.domain,
region_scale_area=True)
# Putting flatten values into xb and std
xb[int(control_slice) * tracer.hresoldim * tracer.vresoldim:
int(control_slice + 1) * tracer.hresoldim * tracer.vresoldim
] += xstack.flatten() * xb_scale + xb_value
std[int(control_slice) * tracer.hresoldim * tracer.vresoldim:
int(control_slice + 1) * tracer.hresoldim * tracer.vresoldim
] += np.abs(xerr).flatten()
# Keep in memory how much inputs values were used for each date
nb_inputs[int(control_slice) * tracer.hresoldim * tracer.vresoldim:
int(control_slice + 1) *
tracer.hresoldim * tracer.vresoldim
] += xcounter.flatten()
# Averaging temporally according to nb inputs
xb /= nb_inputs
std /= nb_inputs
xb = np.where(nb_inputs == 0, 0, xb)
std = np.where(nb_inputs == 0, 0, std)
# TODO: generalize for regions and aggregated temporal resolutions
return xb, std
