Source code for pycif.plugins.domains.unstructured_NetCDF.read_domain

import os
import numpy as np
import xarray as xr

from ..gridded_NetCDF.utils import find_coord
from ..gridded_NetCDF.read_domain import read_vcoord




def read_grid(domain, **kwargs):
    """Read Dynamico grid from a netCDF file.

    Args:
        domain (Domain): The domain object.
    """
    domain.unstructured_domain = True
    file_path = os.path.join(domain.dir, domain.file)

    with xr.open_dataset(file_path) as ds:
        # Getting cell centers coordinates variable names
        lat_coord_name = find_coord(ds, domain.latitude_varname, 'lat')
        lon_coord_name = find_coord(ds, domain.longitude_varname, 'lon')

        # Getting cell centers coordinates
        lat = ds[lat_coord_name]
        lon = ds[lon_coord_name]

        # Checking cell centers coordinates dimensions
        if lat.ndim != 1:
            raise ValueError(
                f"latitude centers variable '{lat_coord_name}' is not 1D.")
        if lon.ndim != 1:
            raise ValueError(
                f"longitude centers variable '{lon_coord_name}' is not 1D.")
        if lat.dims != lon.dims:
            raise ValueError("latitude and longitude centers variables have "
                             "different dimension names.")

        # Getting cell vertices coordinate variable names
        if hasattr(domain, "latitude_bounds_varname"):
            lat_bnds_coord_name = domain.latitude_bounds_varname
        elif "bounds" in lat.attrs:
            lat_bnds_coord_name = ds[lat_coord_name].attrs["bounds"]
        else:
            raise ValueError("could not find latitude bounds coordinate, "
                             "please use the 'latitude_bounds_varname' argument.")

        if lat_bnds_coord_name not in ds:
            raise KeyError("could not find latitude bounds coordinate "
                           f"'{lat_bnds_coord_name}' in file '{file_path}'.")

        if hasattr(domain, "longitude_bounds_varname"):
            lon_bnds_coord_name = domain.longitude_bounds_varname
        elif "bounds" in lon.attrs:
            lon_bnds_coord_name = ds[lon_coord_name].attrs["bounds"]
        else:
            raise ValueError("could not find longitude bounds coordinate, "
                             "please use the 'longitude_bounds_varname' argument.")

        if lon_bnds_coord_name not in ds:
            raise KeyError("could not find latitude bounds coordinate "
                           f"'{lon_bnds_coord_name}' in file '{file_path}'.")

        # Getting cell vertices coordinates
        lat_vertices = ds[lat_bnds_coord_name]
        lon_vertices = ds[lon_bnds_coord_name]

        # Checking cell vertices coordinates dimensions
        if lat_vertices.ndim != 2:
            raise ValueError(
                f"latitude bounds variable '{lat_coord_name}' is not 2D.")
        if lon_vertices.ndim != 2:
            raise ValueError(
                f"longitude bounds variable '{lon_coord_name}' is not 2D.")
        if lat_vertices.dims != lon_vertices.dims:
            raise ValueError("latitude and longitude bounds variables have "
                             "different dimension names.")

        # Dimensions
        cell_dim_name = lat.dims[0]
        if cell_dim_name not in lat_vertices.dims:
            raise ValueError("cell dimension '{cell_dim_name}' not found in "
                             "latitude/longitude bounds variable dimensions.")
        if cell_dim_name == lat_vertices.dims:
            lat_vertices = lat_vertices.T
            lon_vertices = lon_vertices.T

        # Assigning to domain object
        domain.unstructured_domain = True
        domain.nlon = ds.sizes[cell_dim_name]
        domain.nlat = 1
        domain.zlat = lat.values[np.newaxis, :]
        domain.zlon = lon.values[np.newaxis, :]
        domain.zlatc = lat_vertices.values
        domain.zlonc = lon_vertices.values

        read_vcoord(domain)