Source code for pycif.plugins.datastreams.fluxes.GCP_1x1_N2O.get_domain

import datetime
import glob
import os
import numpy as np
import xarray as xr
from logging import debug
from .....utils.classes.setup import Setup
from .....utils.classes.domains import Domain
import copy
from .....utils.check.errclass import CifError
from .....utils.hdf5 import _hdf5_lock




def get_domain(ref_dir, ref_file, input_interval, target_dir, tracer=None):
    """Build the horizontal domain from a GCP N2O reference file.

    Locates a reference NetCDF file (either an exact name match for
    ``ref_file`` in ``ref_dir``, or the first file whose name parses
    against ``ref_file`` as a date format), reads its ``LON``/``LAT``
    variables (forcing latitude into increasing order), and computes cell
    corners assuming a regular grid spacing (a known limitation: this is
    not valid for irregular grids). Builds a single-level, longitude-cyclic
    domain.

    Args:
        ref_dir (str): directory where the original files are found.
        ref_file (str): (template) name of the original files.
        input_interval (list): unused, kept for interface compatibility.
        target_dir (str): unused, kept for interface compatibility.
        tracer: unused, kept for interface compatibility.

    Returns:
        Domain: a single-level, longitude-cyclic domain.

    Raises:
        CifError: if no reference file could be found in ``ref_dir``.
    """

    # Looking for a reference file to read lon/lat in
    list_file = glob.glob(f"{ref_dir}/*nc")

    domain_file = None
    # Either a file is specified in the Yaml
    if ref_file in list_file:
        domain_file = f"{ref_dir}/{ref_file}"

    # Or loop over available file regarding file pattern
    else:
        for flx_file in list_file:
            try:
                date = datetime.datetime.strptime(
                    os.path.basename(flx_file), ref_file
                )
                domain_file = flx_file
                break
            except ValueError:
                continue

    if domain_file is None:
        raise CifError(
            "GCP domain could not be initialized as no file was found"
        )
    debug(f'Domain file for GCP fluxes: {domain_file}')

    with _hdf5_lock:
        nc = xr.open_dataset(domain_file, decode_times=False)

        # Read lon/lat in domain_file
        lon = nc["LON"]
        lat = nc["LAT"]
        # must be increasing
        if lat[1] < lat[0]:
            lat2 = np.flip(lat)
            lat = copy.copy(lat2)
        nlon = lon.size
        nlat = lat.size

        lon_min = lon.min()
        lon_max = lon.max()
        lat_min = lat.min()
        lat_max = lat.max()

        # Compute corners #WARNING: only valid for regular grid in lat/lon -> to generalize
        dx = (lon[1] - lon[0]) / 2
        dy = (lat[1] - lat[0]) / 2
        lonc = np.linspace(lon_min - dx/2., lon_max + dx/2., nlon + 1)
        latc = np.linspace(lat_min - dy/2., lat_max + dy/2., nlat + 1)

        zlon, zlat = np.meshgrid(lon, lat)

    zlonc, zlatc = np.meshgrid(lonc, latc)

    # Read vertical information in domain_file
    # if relevant for emissions
    # XXXXXXXXX comment voir si une dimension verticale est presente?
    # else: dummy vertical
    punit = "Pa"
    nlevs = 1
    sigma_a_mid = np.array([0])
    sigma_b_mid = np.array([1])

    # Put it to a domain Plugin
    domain = Domain(nlon=nlon, nlat=nlat,
                    zlon=zlon, zlat=zlat,
                    zlonc=zlonc, zlatc=zlatc,
                    nlev=nlevs, pressure_unit="Pa",
                    sigma_b_mid=sigma_b_mid, sigma_a_mid=sigma_a_mid,
                    lon_cyclic=True)

    return domain