# Copyright Iris contributors # # This file is part of Iris and is released under the BSD license. # See LICENSE in the root of the repository for full licensing details. r"""Allow the construction of :class:`~iris.experimental.ugrid.mesh.Mesh`. Extensions to Iris' NetCDF loading to allow the construction of :class:`~iris.experimental.ugrid.mesh.Mesh` from UGRID data in the file. Eventual destination: :mod:`iris.fileformats.netcdf`. """ from contextlib import contextmanager from itertools import groupby from pathlib import Path import threading import warnings from ...config import get_logger from ...coords import AuxCoord from ...fileformats._nc_load_rules.helpers import get_attr_units, get_names from ...fileformats.netcdf import loader as nc_loader from ...io import decode_uri, expand_filespecs from ...util import guess_coord_axis from ...warnings import IrisCfWarning, IrisDefaultingWarning, IrisIgnoringWarning from .cf import ( CFUGridAuxiliaryCoordinateVariable, CFUGridConnectivityVariable, CFUGridMeshVariable, CFUGridReader, ) from .mesh import Connectivity, Mesh # Configure the logger. logger = get_logger(__name__, propagate=True, handler=False) class _WarnComboCfDefaulting(IrisCfWarning, IrisDefaultingWarning): """One-off combination of warning classes - enhances user filtering.""" pass class _WarnComboCfDefaultingIgnoring(_WarnComboCfDefaulting, IrisIgnoringWarning): """One-off combination of warning classes - enhances user filtering.""" pass class ParseUGridOnLoad(threading.local): def __init__(self): """Thead-safe state to enable UGRID-aware NetCDF loading. 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 @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): """Mesh from cf, 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 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 def load_meshes(uris, var_name=None): r"""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` 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` 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 == nc_loader.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: with CFUGridReader(source) as cf_reader: meshes_dict = _meshes_from_cf(cf_reader) 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`. 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 = nc_loader._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`. 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 = nc_loader._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`. 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'." warnings.warn( cf_role_message, category=_WarnComboCfDefaulting, ) 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." ) warnings.warn(msg, category=_WarnComboCfDefaulting) 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." ) warnings.warn( msg, category=_WarnComboCfDefaultingIgnoring, ) 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: nc_loader._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`. 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, } location = getattr(cf_var, "location", "") if location is None or location not in element_dimensions: # We should probably issue warnings and recover, but that is too much # work. Raising a more intelligible error is easy to do though. msg = ( f"mesh data variable {cf_var.name!r} has an invalid " f"location={location!r}." ) raise ValueError(msg) mesh_dim_name = element_dimensions.get(location) if mesh_dim_name is None: msg = f"mesh {mesh.name!r} has no {location} dimension." raise ValueError(msg) if mesh_dim_name in cf_var.dimensions: mesh_dim = cf_var.dimensions.index(mesh_dim_name) else: msg = ( f"mesh data variable {cf_var.name!r} does not have the " f"{location} mesh dimension {mesh_dim_name!r}, in its dimensions." ) raise ValueError(msg) mesh_coords = mesh.to_MeshCoords(location=cf_var.location) return mesh_coords, mesh_dim