Source code for pycif.plugins.datastreams.fluxes.orchidee.read

import numpy as np
import xarray as xr
from logging import debug
import xarray as xr
from netCDF4 import Dataset
import pandas as pd

from .....utils.hdf5 import _hdf5_lock




def read(
        self,
        name,
        varnames,
        dates,
        files,
        interpol_flx=False,
        tracer=None,
        model=None,
        ddi=None,
        **kwargs
):
    """Get fluxes from raw ORCHIDEE files and load them into a pyCIF variable.

    For each requested date/file pair, determines the file's time values
    (reading the ``time`` variable, or synthesizing it from
    `tracer.timeresol` if absent), shifting to period start if the values
    mark period midpoints. If `tracer.interpol_resolution` is not set,
    selects the exact matching time index; otherwise, linearly interpolates
    the variable between the two time steps surrounding the midpoint of the
    requested interval.

    Args:
        name (str): name of the component
        varnames (list[str]): original names of variables to read; use `name`
            if `varnames` is empty
        dates (list): list of the date intervals to extract
        files (list): list of the files matching dates
        interpol_flx (bool): Unused directly, kept for interface
            consistency with other flux plugins.
        tracer: The flux tracer plugin, providing ``timeresol`` and,
            optionally, ``interpol_resolution``.
        model: Unused directly, kept for interface consistency.
        ddi: Unused directly, kept for interface consistency.

    Return:
        xr.DataArray: the actual data with dimension:
            time, levels, latitudes, longitudes

    """

    var2extract = varnames if varnames != "" else name

    # Loop over dates/files and import data
    data = []
    out_dates = []
    for dd, ff in zip(dates, files):
        debug(
            f"Reading the file {ff} for the date interval {dd}"
        )

        # Read the file to fetch dates
        with _hdf5_lock:
            with Dataset(ff, "r") as f:
                isvar = "time" in f.variables

            times = xr.open_dataset(ff)["time"].values

        # Define dates if not a variable
        if not isvar:
            times = pd.date_range(ddi, periods=len(times),
                                  freq=tracer.timeresol).values[:, np.newaxis]

        # Shift dates if in variables as the middle of periods is specified
        freq = np.unique(np.diff(times))
        if isvar and not hasattr(tracer, "interpol_resolution"):
            times -= freq / 2
        
        # Fetch correct index
        if not hasattr(tracer, "interpol_resolution"):
            dref = np.datetime64(dd[0])
            ind_time = np.where(times == dref)[0][0]
            with _hdf5_lock:
                data.append(xr.open_dataset(ff)[var2extract][ind_time].values)
        else:
            dtref = np.datetime64(dd[1]) - np.datetime64(dd[0])
            dref = np.datetime64(dd[0]) + dtref / 2
            ind_time = np.argmax(times > dref)
            d0 = times[ind_time - 1]
            d1 = times[ind_time]
            with _hdf5_lock:
                data0 = xr.open_dataset(ff)[var2extract][ind_time - 1].values
                data1 = xr.open_dataset(ff)[var2extract][ind_time].values
            data.append(
                (dref - d0) / (d1 - d0) * data1
                + (d1 - dref) / (d1 - d0) * data0)

        out_dates.append(dd[0])

    # if only one level for emissions, create the axis
    dataout = np.array(data)[:, np.newaxis]
    dataout[np.isnan(dataout)] = 0
    xmod = xr.DataArray(
        dataout,
        coords={"time": out_dates},
        dims=("time", "lev", "lat", "lon"),
    )
    return xmod