Source code for pycif.plugins.datastreams.fluxes.GridFED.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.check.errclass import CifError
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 GridFED files and load them into a pyCIF variable.

    For each requested date/file pair, locates the matching time index,
    sums the values of `tracer.sectors` for the `tracer.group_name` species
    group, converts the result from kgCO2/month to kgCO2/h, and finally
    divides by the domain cell areas (computed on the fly if not already
    available) to obtain kgCO2/m2/h.

    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, kept for interface consistency with
            other flux plugins.
        tracer: The flux tracer plugin, providing ``group_name``,
            ``sectors`` and ``domain``.
        model: Unused, kept for interface consistency.
        ddi: Unused, kept for interface consistency.

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

    Raises:
        CifError: If a required sector in `tracer.sectors` is missing from
            a file.
    """

    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:
            times = xr.open_dataset(ff)["time"].values
            ind_time = np.where(times == np.datetime64(dd[0]))[0][0]

            ds = xr.open_dataset(ff, group=tracer.group_name)

            data_tmp = 0
            for sector in tracer.sectors:
                if sector not in ds.variables:
                    raise CifError(f"Required sector {sector} while only the following sectors are available in file {ff}: {list(ds.variables.keys())}")

                data_tmp += ds[sector][ind_time].values

        # Convert outputs from kgCO2/month to kgCO2/h
        data_tmp /= pd.DatetimeIndex([dd[0]]).days_in_month[0] * 24
        
        data.append(data_tmp)
        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"),
    )
    
    # Convert outputs from kgCO2/h to kgCO2/m2/h
    if not hasattr(tracer.domain, "areas"):
        tracer.domain.calc_areas(**kwargs)
    
    xmod /= tracer.domain.areas
    
    return xmod