import numpy as np
import pandas as pd
import copy
from ......utils.datastores.empty import init_empty
from ......utils.check.errclass import CifError
[docs]
def make_obs(self, ddi, dict_datastores, runsubdir, mode, list_tracer, input_type, do_simu=True):
"""Write the CHIMERE observation text file ``obs.txt`` for one sub-period.
Accumulates observation metadata (grid indices, time steps, species,
dtstep) into the CHIMERE-formatted ``obs.txt`` file. For adjoint runs
the departure vector (``adj_out``) is appended as an extra column.
Handles:
* **Auxiliary outputs** (pressure, temperature, …): only stores count
metadata in ``self.auxiliary_indexes``; returns early without writing.
* **do_simu=False**: updates obs-count bookkeeping without writing a
file (used when the forward output was already cached).
* Multi-timestep observations (``dtstep > 1``) are expanded row-by-row
so each line covers exactly one CHIMERE sub-time-step.
* Coordinate and level indices are converted from Python (0-based) to
Fortran (1-based) before writing.
Args:
self: CHIMERE model plugin instance (carries ``nhour``,
``subtstep``, ``iniobs``, ``reset_obs``, ``chemistry.outspecies``,
and ``output_components``).
ddi (datetime): sub-simulation period start.
dict_datastores (dict): species-keyed CIF data-store mapping.
runsubdir (str): path to the period run directory.
mode (str): ``'fwd'``, ``'tl'``, or ``'adj'``.
list_tracer (list): species names to write in this call.
input_type (str): output component name (``'concs'``, ``'pressure'``,
etc.).
do_simu (bool): if ``False``, skip writing and only update bookkeeping.
"""
# If empty datastore, do nothing
if np.all([type(dict_datastores[d]) == dict for d in dict_datastores]):
return
if np.all([dict_datastores[d].size == 0 for d in dict_datastores]):
return
# Re-initialize auxiliary_indexes
# if first species processed for the current period
if not hasattr(self, "auxiliary_indexes"):
self.auxiliary_indexes = {
ddi: {
outcomp: {spec: {}
for spec in self.chemistry.outspecies.attributes}
for outcomp in self.output_components
}
}
self.chunk_indexes = {
ddi: {
outcomp: {spec: {}
for spec in self.chemistry.outspecies.attributes}
for outcomp in self.output_components
}
}
if self.reset_obs[ddi]:
self.auxiliary_indexes[ddi] = {
outcomp: {spec: {}
for spec in self.chemistry.outspecies.attributes}
for outcomp in self.output_components
}
self.chunk_indexes[ddi] = {
outcomp: {spec: {}
for spec in self.chemistry.outspecies.attributes}
for outcomp in self.output_components
}
self.reset_obs[ddi] = False
# If input type is not "concs", i.e., auxiliary data such as pressure,
# Keep in memory the length of the datastore, but not more
for trcr in list_tracer:
self.auxiliary_indexes[ddi][input_type][trcr] = \
[len(dict_datastores[trcr])]
# If values to extract come from several transforms, keep in memory
# chunks
if "itransform" in dict_datastores[trcr]["metadata"]:
self.chunk_indexes[ddi][input_type][trcr] = \
dict_datastores[trcr]["metadata"]["itransform"].values
if input_type != "concs":
return
# If do not need to do CHIMERE simulation, just update obs datastore
if not do_simu:
if not self.iniobs[ddi]:
self.nbobs_prior[ddi] = 0
self.nbdatatot_prior[ddi] = 0
# Keep in memory that observations were already dumped for that period
self.iniobs[ddi] = True
for trcr in list_tracer:
self.auxiliary_indexes[ddi][input_type][trcr].extend(
[self.nbobs_prior[ddi], self.nbdatatot_prior[ddi]]
)
self.nbobs_prior[ddi] += len(dict_datastores[trcr])
self.nbdatatot_prior[ddi] += len(dict_datastores[trcr])
return
# Load previous obs if exists
obs_file = f"{runsubdir}/obs.txt"
nbobs_prior = 0
nbdatatot_prior = 0
data_prior = pd.DataFrame(np.empty((0, 5)))
if self.iniobs[ddi]:
nbobs_prior, nbdatatot_prior = pd.read_csv(
obs_file, sep="\s+", header=None, nrows=1
).values[0]
data_prior = pd.concat(
[data_prior, pd.read_csv(
obs_file, sep="\s+", skiprows=1, header=None)]
)
# Loop over tracers
for trcr in list_tracer:
datastore = dict_datastores[trcr]
mask = (
datastore["metadata"]["parameter"]
.str.replace("__sample#", "")
.str.upper()
.isin([s.upper() for s in self.chemistry.outspecies.attributes])
)
data2write = datastore.loc[mask]
# Replace parameter if batch computing
data2write.loc[:, ("metadata", "parameter")] = data2write["metadata"][
"parameter"
].str.replace("__sample#", "")
# For adjoint, check that there is no NaN values
if mode == "adj":
if not np.all(pd.notnull(data2write["maindata"].loc[:, "adj_out"])):
raise CifError(
"WARNING: pycif will drive CHIMERE adjoint "
"with NaNs values! Check prior informations"
)
self.auxiliary_indexes[ddi][input_type][trcr].extend(
[nbobs_prior, nbdatatot_prior]
)
# Reformat data to wrtie
nbobs_prior += len(data2write)
nbdatatot_prior += data2write["metadata"]["dtstep"].sum()
columns = ["tstep", "i", "j", "level", "parameter", "dtstep"]
dataout = copy.deepcopy(data2write["metadata"].loc[:, columns])
# Adding adj_out if adjoint
if mode == "adj":
dataout = pd.concat(
[dataout, data2write["maindata"]["adj_out"]], axis=1)
# Expanding observations spanning several time stamps
if np.any(dataout["dtstep"] != 1):
repeat_index = dataout.reset_index().index.repeat(
dataout["dtstep"].astype(int))
new_tstep = np.zeros(len(repeat_index))
new_tstep[
(np.cumsum(dataout["dtstep"].values) -
dataout["dtstep"].values[0]).astype(int)
] = (np.cumsum(dataout["dtstep"]) - dataout["dtstep"].values[0])
new_tstep = np.arange(len(new_tstep)) - \
np.maximum.accumulate(new_tstep)
dataout = dataout.iloc[repeat_index]
dataout.loc[:, "tstep"] += new_tstep
dataout["dtstep"] = 1
# Convert level and coordinates to fortran
dataout["level"] = np.asarray(dataout["level"]) + 1
dataout["i"] = np.asarray(dataout["i"]) + 1
dataout["j"] = np.asarray(dataout["j"]) + 1
# TODO: deal with altitude; so far, assumes that if no level specified,
# extract first level
mask = ~pd.notnull(dataout["level"]) | (dataout["level"] <= 0)
dataout.loc[mask, "level"] = 1
# Add column with index of current hour
dataout.loc[:, "hour"] = (
np.array(self.nhour[ddi])[dataout["tstep"].values.astype(int)] - 1
)
# Time steps in CHIMERE are time sub time steps in a given hours
# They have been precomputed in self.subtstep
dataout["tstep"] = np.array(self.subtstep[ddi])[
dataout["tstep"].values.astype(int)]
# Parameter should fit OUTPUT_SPECIES in CHIMERE (case sensitive)
output_species = pd.Series(
index=[s.upper() for s in self.chemistry.outspecies.attributes],
data=range(len(self.chemistry.outspecies.attributes)),
)
index_output_species = output_species.loc[dataout["parameter"].str.upper(
)]
dataout.loc[:, "parameter"] = np.array(self.chemistry.outspecies.attributes)[
index_output_species
]
# Concatenate prior values with new values
ind2dump = [7, 0, 1, 2, 3, 4, 6] if mode == "adj" else [
6, 0, 1, 2, 3, 4]
data_prior = pd.concat(
[data_prior, pd.DataFrame(dataout.iloc[:, ind2dump].values)]
).astype({0: "int", 1: "int", 2: "int", 3: "int", 4: "int"})
with open(obs_file, "w") as f:
f.write(str(int(nbobs_prior)) + " " + str(int(nbdatatot_prior)) + "\n")
data_prior.to_csv(f, sep=" ", index=False, header=False)
# Keep in memory that observations were already dumped for that period
self.iniobs[ddi] = True
