import glob
import os
from logging import debug, info
import xarray as xr
import numpy as np
from .....utils.check.errclass import CifKeyError
# D&B output variable names (see src/ncwrite_output.f90) for each pyCIF
# output component declared in ini_mapper.py / __init__.py::output_components
_DB_VARNAME = {
"NEE": "nee",
}
[docs]
def outputs2native(
self, data2dump, input_type, di, df, runsubdir, mode="fwd",
onlyinit=False, check_transforms=False, **kwargs
):
"""Read D&B NetCDF diagnostics back into pyCIF objects.
D&B writes its diagnostics to
``diagout/dalec-bethy_{hourly,daily}-output_<yyyymmdd_start>-<yyyymmdd_end>.nc``
(see ``src/ncwrite_output.f90``), on the ``nsp`` sample-point dimension
(grid-cell x active-PFT combinations, see ``src/dimensions.f90``). The
stream (``hourly``/``daily``) to read is controlled by
``self.output_resolution``.
The output file also carries, on that same ``nsp`` dimension, the
``gidx`` (1-based grid-cell index, D&B's ``map_sample2grid``) and
``pft_fraction`` (D&B's ``fracv``) static fields. These are used here
to recombine the PFT-level sample points of each grid cell into a
single per-cell value (a ``pft_fraction``-weighted average), so the
result lines up with the pyCIF ``domain`` built by
``pycif.plugins.domains.dalecbethy.read_domain`` (``domain.active``,
on the same grid-cell indexing as ``gidx``).
Args:
self: the dalecbethy model Plugin.
data2dump (dict): output data structure to fill for every
component/tracer declared in the mapper.
input_type (str): the type of model outputs to be processed
(redundant with the components of ``data2dump``).
di, df (datetime.datetime): start/end date of the present
sub-simulation.
runsubdir (str): path to the present sub-simulation work directory.
mode (str): running mode; one of "fwd", "tl" and "adj".
Returns:
dict: a dictionary with structure the components/tracers to be
extracted, with each component's per-grid-cell ``xarray.DataArray``
(PFT sample points recombined onto ``self.domain``'s grid cells).
"""
if mode not in ("fwd",):
debug(f"dalecbethy outputs2native: nothing to do for mode={mode!r}")
return {}
stream = "hourly" if self.output_resolution == "hourly" else "daily"
pattern = os.path.join(
runsubdir, "diagout", f"dalec-bethy_{stream}-output_*.nc")
matches = sorted(glob.glob(pattern))
if not matches:
debug(f"No D&B {stream} output file found matching {pattern}")
return {}
ncfile = matches[0]
info(f"Reading D&B output from {ncfile}")
domain = self.domain
active = domain.active
dataout = {}
with xr.open_dataset(ncfile) as ds:
for key in ("gidx", "pft_fraction"):
if key not in ds:
raise CifKeyError(
f"mandatory variable {key!r} not found in {ncfile}; "
f"cannot recombine D&B's sample points back onto the "
f"pyCIF domain's grid cells."
)
# 0-based grid-cell index of each sample point (D&B writes it
# 1-based, see src/static_fields.f90::map_sample2grid), and the
# PFT area fraction to weight it by when recombining sample points
# sharing the same grid cell (src/static_fields.f90::fracv).
gidx = ds["gidx"].values.astype(int) - 1
weight = ds["pft_fraction"].values
denom = np.zeros(domain.ng)
np.add.at(denom, gidx, weight)
for trid in data2dump:
if trid[0] != input_type:
continue
varname = _DB_VARNAME.get(trid[1])
if varname is None or varname not in ds:
debug(
f"Variable for component {trid} (D&B name "
f"{varname!r}) not found in {ncfile}; skipping."
)
continue
dates = ds.time.values
nee_sp = ds[varname].values # (time, nsp)
numer = np.zeros((dates.size, domain.ng))
np.add.at(numer, (slice(None), gidx), weight[np.newaxis, :] * nee_sp)
with np.errstate(invalid="ignore"):
nee_cell = numer / denom[np.newaxis, :]
xmod = xr.DataArray(
nee_cell[:, active][:, np.newaxis, np.newaxis, :],
coords={"time": dates},
dims=("time", "lev", "lat", "lon"),
)
dataout[trid] = {"spec": xmod}
return dataout
