Source code for pycif.plugins.datastreams.fields.netcdf_cams.get_domain

import numpy as np
import glob
import datetime
import os
import pandas as pd
import xarray as xr
from .....utils.classes.setup import Setup
from .....utils.classes.domains import Domain
import copy
from logging import debug
import itertools




def get_domain(ref_dir, ref_file, input_interval, target_dir, tracer=None):
    domain_file = list(itertools.chain.from_iterable(
        tracer.input_files.values()))[0]

    if not os.path.isfile(domain_file):
        raise Exception("Could not initialize the domain "
                        "as no reference file is available. "
                        "Expecting the following file: {}".format(ref_file))

    debug('Domain file for CAMS BCs: {}'.format(domain_file))

    nc = xr.open_dataset(domain_file, decode_times=False)

    # Read lon/lat in domain_file
    lon = nc["longitude"]
    lat = nc["latitude"]

    # must be increasing
    if lat[1] < lat[0]:
        lat = np.flip(lat)

    nlon = lon.size
    nlat = lat.size
    zlon, zlat = np.meshgrid(lon, lat)

    # Compute corners
    dlat = np.unique(np.diff(lat))[0]
    dlon = np.unique(np.diff(lon))[0]

    latc = np.append(lat - dlat / 2, lat[-1] + dlat / 2)
    lonc = np.append(lon - dlon / 2, lon[-1] + dlon / 2)
    zlonc, zlatc = np.meshgrid(lonc, latc)
    zlatc = np.minimum(np.maximum(zlatc, -90), 90)

    # Read vertical information in domain_file
    if getattr(tracer, 'aibi_name', False):
        # if tracer.aibi_name:
        if hasattr(tracer, "aibi_file"):
            aibi_file = tracer.aibi_file
            tab = pd.read_csv(aibi_file, sep=',', header=0, usecols=[0, 1, 2])
            sigma_a = tab.values[::-1, 1]
            sigma_b = tab.values[::-1, 2]
        else:
            sigma_a = nc["ap"].values
            sigma_b = nc["bp"].values
        nlevs = sigma_b.size - 1
        domain = Domain(
            nlon=nlon, nlat=nlat, zlon=zlon, zlat=zlat,
            zlonc=zlonc, zlatc=zlatc,
            pressure_unit="Pa", lon_cyclic=True, nlev=nlevs,
            sigma_a=sigma_a, sigma_b=sigma_b
        )

    else:
        sigma_a_mid = nc["hyam"].values
        sigma_b_mid = nc["hybm"].values
        nlevs = sigma_a_mid.size
        domain = Domain(
            nlon=nlon, nlat=nlat, zlon=zlon, zlat=zlat,
            zlonc=zlonc, zlatc=zlatc,
            pressure_unit="Pa", lon_cyclic=True, nlev=nlevs,
            sigma_a_mid=sigma_a_mid, sigma_b_mid=sigma_b_mid
        )

    # Initializes domain
    return domain