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




def get_domain(ref_dir, ref_file, input_interval, target_dir, tracer=None):

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

    domain_file = None
    # Either a file is specified in the Yaml
    if ref_file in list_file:
        domain_file = "{}/{}".format(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 Exception(
            "GCP domain could not be initialized as no file was found"
        )
    debug('Domain file for GCP fluxes: {}'.format(domain_file))

    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