You are viewing the latest unreleased documentation v3.3.dev0. You may prefer a stable version.

Source code for iris.experimental.ugrid.load

# Copyright Iris contributors
#
# This file is part of Iris and is released under the LGPL license.
# See COPYING and COPYING.LESSER in the root of the repository for full
# licensing details.

"""
Extensions to Iris' NetCDF loading to allow the construction of
:class:`~iris.experimental.ugrid.mesh.Mesh`\\ es from UGRID data in the file.

Eventual destination: :mod:`iris.fileformats.netcdf` (plan to split that module
into ``load`` and ``save`` in future).

"""
from contextlib import contextmanager
from itertools import groupby
from pathlib import Path
import threading

from ...config import get_logger
from ...coords import AuxCoord
from ...fileformats import netcdf
from ...fileformats._nc_load_rules.helpers import get_attr_units, get_names
from ...io import decode_uri, expand_filespecs
from ...util import guess_coord_axis
from .cf import (
    CFUGridAuxiliaryCoordinateVariable,
    CFUGridConnectivityVariable,
    CFUGridMeshVariable,
    CFUGridReader,
)
from .mesh import Connectivity, Mesh

# Configure the logger.
logger = get_logger(__name__, propagate=True, handler=False)


[docs]class ParseUGridOnLoad(threading.local): def __init__(self): """ A flag for dictating whether to use the experimental UGRID-aware version of Iris NetCDF loading. Object is thread-safe. Use via the run-time switch :const:`~iris.experimental.ugrid.load.PARSE_UGRID_ON_LOAD`. Use :meth:`context` to temporarily activate. .. seealso:: The UGRID Conventions, https://ugrid-conventions.github.io/ugrid-conventions/ """ self._state = False def __bool__(self): return self._state
[docs] @contextmanager def context(self): """ Temporarily activate experimental UGRID-aware NetCDF loading. Use the standard Iris loading API while within the context manager. If the loaded file(s) include any UGRID content, this will be parsed and attached to the resultant cube(s) accordingly. Use via the run-time switch :const:`~iris.experimental.ugrid.load.PARSE_UGRID_ON_LOAD`. For example:: with PARSE_UGRID_ON_LOAD.context(): my_cube_list = iris.load([my_file_path, my_file_path2], constraint=my_constraint, callback=my_callback) """ try: self._state = True yield finally: self._state = False
#: Run-time switch for experimental UGRID-aware NetCDF loading. See :class:`~iris.experimental.ugrid.load.ParseUGridOnLoad`. PARSE_UGRID_ON_LOAD = ParseUGridOnLoad() def _meshes_from_cf(cf_reader): """ Common behaviour for extracting meshes from a CFReader. Simple now, but expected to increase in complexity as Mesh sharing develops. """ # Mesh instances are shared between file phenomena. # TODO: more sophisticated Mesh sharing between files. # TODO: access external Mesh cache? mesh_vars = cf_reader.cf_group.meshes meshes = { name: _build_mesh(cf_reader, var, cf_reader.filename) for name, var in mesh_vars.items() } return meshes
[docs]def load_mesh(uris, var_name=None): """ Load a single :class:`~iris.experimental.ugrid.mesh.Mesh` object from one or more NetCDF files. Raises an error if more/less than one :class:`~iris.experimental.ugrid.mesh.Mesh` is found. Parameters ---------- uris : str or iterable of str One or more filenames/URI's. Filenames can include wildcards. Any URI's must support OpenDAP. var_name : str, optional Only return a :class:`~iris.experimental.ugrid.mesh.Mesh` if its var_name matches this value. Returns ------- :class:`iris.experimental.ugrid.mesh.Mesh` """ meshes_result = load_meshes(uris, var_name) result = set([mesh for file in meshes_result.values() for mesh in file]) mesh_count = len(result) if mesh_count != 1: message = ( f"Expecting 1 mesh, but input file(s) produced: {mesh_count} ." ) raise ValueError(message) return result.pop() # Return the single element
[docs]def load_meshes(uris, var_name=None): """ Load :class:`~iris.experimental.ugrid.mesh.Mesh` objects from one or more NetCDF files. Parameters ---------- uris : str or iterable of str One or more filenames/URI's. Filenames can include wildcards. Any URI's must support OpenDAP. var_name : str, optional Only return :class:`~iris.experimental.ugrid.mesh.Mesh`\\ es that have var_names matching this value. Returns ------- dict A dictionary mapping each mesh-containing file path/URL in the input ``uris`` to a list of the :class:`~iris.experimental.ugrid.mesh.Mesh`\\ es returned from each. """ # TODO: rationalise UGRID/mesh handling once experimental.ugrid is folded # into standard behaviour. # No constraints or callbacks supported - these assume they are operating # on a Cube. from ...fileformats import FORMAT_AGENT if not PARSE_UGRID_ON_LOAD: # Explicit behaviour, consistent with netcdf.load_cubes(), rather than # an invisible assumption. message = ( f"PARSE_UGRID_ON_LOAD is {bool(PARSE_UGRID_ON_LOAD)}. Must be " f"True to enable mesh loading." ) raise ValueError(message) if isinstance(uris, str): uris = [uris] # Group collections of uris by their iris handler # Create list of tuples relating schemes to part names. uri_tuples = sorted(decode_uri(uri) for uri in uris) valid_sources = [] for scheme, groups in groupby(uri_tuples, key=lambda x: x[0]): # Call each scheme handler with the appropriate URIs if scheme == "file": filenames = [x[1] for x in groups] sources = expand_filespecs(filenames) elif scheme in ["http", "https"]: sources = [":".join(x) for x in groups] else: message = f"Iris cannot handle the URI scheme: {scheme}" raise ValueError(message) for source in sources: if scheme == "file": with open(source, "rb") as fh: handling_format_spec = FORMAT_AGENT.get_spec( Path(source).name, fh ) else: handling_format_spec = FORMAT_AGENT.get_spec(source, None) if handling_format_spec.handler == netcdf.load_cubes: valid_sources.append(source) else: message = f"Ignoring non-NetCDF file: {source}" logger.info(msg=message, extra=dict(cls=None)) result = {} for source in valid_sources: meshes_dict = _meshes_from_cf(CFUGridReader(source)) meshes = list(meshes_dict.values()) if var_name is not None: meshes = list(filter(lambda m: m.var_name == var_name, meshes)) if meshes: result[source] = meshes return result
############ # Object construction. # Helper functions, supporting netcdf.load_cubes ONLY, expected to # altered/moved when pyke is removed. def _build_aux_coord(coord_var, file_path): """ Construct a :class:`~iris.coords.AuxCoord` from a given :class:`~iris.experimental.ugrid.cf.CFUGridAuxiliaryCoordinateVariable`, and guess its mesh axis. todo: integrate with standard loading API post-pyke. """ # TODO: integrate with standard saving API when no longer 'experimental'. assert isinstance(coord_var, CFUGridAuxiliaryCoordinateVariable) attributes = {} attr_units = get_attr_units(coord_var, attributes) points_data = netcdf._get_cf_var_data(coord_var, file_path) # Bounds will not be loaded: # Bounds may be present, but the UGRID conventions state this would # always be duplication of the same info provided by the mandatory # connectivities. # Fetch climatological - not allowed for a Mesh, but loading it will # mean an informative error gets raised. climatological = False # TODO: use CF_ATTR_CLIMATOLOGY on re-integration, when no longer # 'experimental'. attr_climatology = getattr(coord_var, "climatology", None) if attr_climatology is not None: climatology_vars = coord_var.cf_group.climatology climatological = attr_climatology in climatology_vars standard_name, long_name, var_name = get_names(coord_var, None, attributes) coord = AuxCoord( points_data, standard_name=standard_name, long_name=long_name, var_name=var_name, units=attr_units, attributes=attributes, # TODO: coord_system climatological=climatological, ) axis = guess_coord_axis(coord) if axis is None: if var_name[-2] == "_": # Fall back on UGRID var_name convention. axis = var_name[-1] else: message = f"Cannot guess axis for UGRID coord: {var_name} ." raise ValueError(message) return coord, axis def _build_connectivity(connectivity_var, file_path, element_dims): """ Construct a :class:`~iris.experimental.ugrid.mesh.Connectivity` from a given :class:`~iris.experimental.ugrid.cf.CFUGridConnectivityVariable`, and identify the name of its first dimension. todo: integrate with standard loading API post-pyke. """ # TODO: integrate with standard saving API when no longer 'experimental'. assert isinstance(connectivity_var, CFUGridConnectivityVariable) attributes = {} attr_units = get_attr_units(connectivity_var, attributes) indices_data = netcdf._get_cf_var_data(connectivity_var, file_path) cf_role = connectivity_var.cf_role start_index = connectivity_var.start_index dim_names = connectivity_var.dimensions # Connectivity arrays must have two dimensions. assert len(dim_names) == 2 if dim_names[1] in element_dims: location_axis = 1 else: location_axis = 0 standard_name, long_name, var_name = get_names( connectivity_var, None, attributes ) connectivity = Connectivity( indices=indices_data, cf_role=cf_role, standard_name=standard_name, long_name=long_name, var_name=var_name, units=attr_units, attributes=attributes, start_index=start_index, location_axis=location_axis, ) return connectivity, dim_names[0] def _build_mesh(cf, mesh_var, file_path): """ Construct a :class:`~iris.experimental.ugrid.mesh.Mesh` from a given :class:`~iris.experimental.ugrid.cf.CFUGridMeshVariable`. todo: integrate with standard loading API post-pyke. """ # TODO: integrate with standard saving API when no longer 'experimental'. assert isinstance(mesh_var, CFUGridMeshVariable) attributes = {} attr_units = get_attr_units(mesh_var, attributes) cf_role_message = None if not hasattr(mesh_var, "cf_role"): cf_role_message = f"{mesh_var.cf_name} has no cf_role attribute." cf_role = "mesh_topology" else: cf_role = getattr(mesh_var, "cf_role") if cf_role != "mesh_topology": cf_role_message = ( f"{mesh_var.cf_name} has an inappropriate cf_role: {cf_role}." ) if cf_role_message: cf_role_message += " Correcting to 'mesh_topology'." # TODO: reconsider logging level when we have consistent practice. logger.warning(cf_role_message, extra=dict(cls=None)) if hasattr(mesh_var, "volume_node_connectivity"): topology_dimension = 3 elif hasattr(mesh_var, "face_node_connectivity"): topology_dimension = 2 elif hasattr(mesh_var, "edge_node_connectivity"): topology_dimension = 1 else: # Nodes only. We aren't sure yet whether this is a valid option. topology_dimension = 0 if not hasattr(mesh_var, "topology_dimension"): msg = ( f"Mesh variable {mesh_var.cf_name} has no 'topology_dimension'" f" : *Assuming* topology_dimension={topology_dimension}" ", consistent with the attached connectivities." ) # TODO: reconsider logging level when we have consistent practice. logger.warning(msg, extra=dict(cls=None)) else: quoted_topology_dimension = mesh_var.topology_dimension if quoted_topology_dimension != topology_dimension: msg = ( f"*Assuming* 'topology_dimension'={topology_dimension}" f", from the attached connectivities of the mesh variable " f"{mesh_var.cf_name}. However, " f"{mesh_var.cf_name}:topology_dimension = " f"{quoted_topology_dimension}" " -- ignoring this as it is inconsistent." ) # TODO: reconsider logging level when we have consistent practice. logger.warning(msg=msg, extra=dict(cls=None)) node_dimension = None edge_dimension = getattr(mesh_var, "edge_dimension", None) face_dimension = getattr(mesh_var, "face_dimension", None) node_coord_args = [] edge_coord_args = [] face_coord_args = [] for coord_var in mesh_var.cf_group.ugrid_coords.values(): coord_and_axis = _build_aux_coord(coord_var, file_path) coord = coord_and_axis[0] if coord.var_name in mesh_var.node_coordinates.split(): node_coord_args.append(coord_and_axis) node_dimension = coord_var.dimensions[0] elif ( coord.var_name in getattr(mesh_var, "edge_coordinates", "").split() ): edge_coord_args.append(coord_and_axis) elif ( coord.var_name in getattr(mesh_var, "face_coordinates", "").split() ): face_coord_args.append(coord_and_axis) # TODO: support volume_coordinates. else: message = ( f"Invalid UGRID coord: {coord.var_name} . Must be either a" f"node_, edge_ or face_coordinate." ) raise ValueError(message) if node_dimension is None: message = ( "'node_dimension' could not be identified from mesh node " "coordinates." ) raise ValueError(message) # Used for detecting transposed connectivities. element_dims = (edge_dimension, face_dimension) connectivity_args = [] for connectivity_var in mesh_var.cf_group.connectivities.values(): connectivity, first_dim_name = _build_connectivity( connectivity_var, file_path, element_dims ) assert connectivity.var_name == getattr(mesh_var, connectivity.cf_role) connectivity_args.append(connectivity) # If the mesh_var has not supplied the dimension name, it is safe to # fall back on the connectivity's first dimension's name. if edge_dimension is None and connectivity.location == "edge": edge_dimension = first_dim_name if face_dimension is None and connectivity.location == "face": face_dimension = first_dim_name standard_name, long_name, var_name = get_names(mesh_var, None, attributes) mesh = Mesh( topology_dimension=topology_dimension, node_coords_and_axes=node_coord_args, connectivities=connectivity_args, edge_coords_and_axes=edge_coord_args, face_coords_and_axes=face_coord_args, standard_name=standard_name, long_name=long_name, var_name=var_name, units=attr_units, attributes=attributes, node_dimension=node_dimension, edge_dimension=edge_dimension, face_dimension=face_dimension, ) mesh_elements = ( list(mesh.all_coords) + list(mesh.all_connectivities) + [mesh] ) mesh_elements = filter(None, mesh_elements) for iris_object in mesh_elements: netcdf._add_unused_attributes( iris_object, cf.cf_group[iris_object.var_name] ) return mesh def _build_mesh_coords(mesh, cf_var): """ Construct a tuple of :class:`~iris.experimental.ugrid.mesh.MeshCoord` using from a given :class:`~iris.experimental.ugrid.mesh.Mesh` and :class:`~iris.fileformats.cf.CFVariable`. todo: integrate with standard loading API post-pyke. """ # TODO: integrate with standard saving API when no longer 'experimental'. # Identify the cube's mesh dimension, for attaching MeshCoords. element_dimensions = { "node": mesh.node_dimension, "edge": mesh.edge_dimension, "face": mesh.face_dimension, } mesh_dim_name = element_dimensions[cf_var.location] # (Only expecting 1 mesh dimension per cf_var). mesh_dim = cf_var.dimensions.index(mesh_dim_name) mesh_coords = mesh.to_MeshCoords(location=cf_var.location) return mesh_coords, mesh_dim