diff --git a/docs/src/common_links.inc b/docs/src/common_links.inc index eb1ea60b7a..67fc493e3e 100644 --- a/docs/src/common_links.inc +++ b/docs/src/common_links.inc @@ -37,6 +37,7 @@ .. _test-iris-imagehash: https://github.com/SciTools/test-iris-imagehash .. _using git: https://docs.github.com/en/github/using-git .. _requirements/ci/: https://github.com/SciTools/iris/tree/main/requirements/ci +.. _CF-UGRID: https://ugrid-conventions.github.io/ugrid-conventions/ .. comment diff --git a/docs/src/conf.py b/docs/src/conf.py index 5a436f86cb..138e0aa7fc 100644 --- a/docs/src/conf.py +++ b/docs/src/conf.py @@ -200,7 +200,9 @@ def _dotv(version): # -- copybutton extension ----------------------------------------------------- # See https://sphinx-copybutton.readthedocs.io/en/latest/ -copybutton_prompt_text = ">>> " +copybutton_prompt_text = r">>> |\.\.\. " +copybutton_prompt_is_regexp = True +copybutton_line_continuation_character = "\\" # sphinx.ext.todo configuration ----------------------------------------------- # See https://www.sphinx-doc.org/en/master/usage/extensions/todo.html diff --git a/docs/src/further_topics/index.rst b/docs/src/further_topics/index.rst index dc162d6a1e..4f9aadbc79 100644 --- a/docs/src/further_topics/index.rst +++ b/docs/src/further_topics/index.rst @@ -20,6 +20,7 @@ that may be of interest to the more advanced or curious user. * :doc:`metadata` * :doc:`lenient_metadata` * :doc:`lenient_maths` +* :ref:`ugrid` .. _GitHub Documentation Issue: https://github.com/SciTools/iris/issues/new?assignees=&labels=New%3A+Documentation%2C+Type%3A+Documentation&template=documentation.md&title= diff --git a/docs/src/further_topics/ugrid/data_model.rst b/docs/src/further_topics/ugrid/data_model.rst new file mode 100644 index 0000000000..4a2f64f627 --- /dev/null +++ b/docs/src/further_topics/ugrid/data_model.rst @@ -0,0 +1,566 @@ +.. include:: ../../common_links.inc + +.. _ugrid model: + +The Mesh Data Model +******************* + +.. important:: + + This page is intended to summarise the essentials that Iris users need + to know about meshes. For exhaustive details on UGRID itself: + `visit the official UGRID conventions site`__. + +Evolution, not revolution +========================= +Mesh support has been designed wherever possible to fit within the existing +Iris model. Meshes concern only the spatial geography of data, and can +optionally be limited to just the horizontal geography (e.g. X and Y). Other +dimensions such as time or ensemble member (and often vertical levels) +retain their familiar structured format. + +The UGRID conventions themselves are designed as an addition to the existing CF +conventions, which are at the core of Iris' philosophy. + +What's Different? +================= + +The mesh format represents data's geography using an **unstructured +mesh**. This has significant pros and cons when compared to a structured grid. + +.. contents:: + :local: + +The Detail +---------- +.. + The diagram images are SVG's, so editable by any graphical software + (e.g. Inkscape). They were originally made in MS PowerPoint. + + Uses the IBM Colour Blind Palette (see + http://ibm-design-language.eu-de.mybluemix.net/design/language/resources/color-library + ) + +Structured Grids (the old world) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Assigning data to locations using a structured grid is essentially an act of +matching coordinate arrays to each dimension of the data array. The data can +also be represented as an area (instead of a point) by including a bounds array +for each coordinate array. :numref:`data_structured_grid` visualises an +example. + +.. _data_structured_grid: +.. figure:: images/data_structured_grid.svg + :alt: Diagram of how data is represented on a structured grid + :align: right + :width: 1280 + + Data on a structured grid. + + 1D coordinate arrays (pink circles) are combined to construct a structured + grid of points (pink crosses). 2D bounds arrays (blue circles) can also be + used to describe the 1D boundaries (blue lines) at either side of each + rank of points; each point therefore having four bounds (x+y, upper+lower), + together describing a quadrilateral area around that point. Data from the + 2D data array (orange circles) can be assigned to these point locations + (orange diamonds) or area locations (orange quads) by matching the relative + positions in the data array to the relative spatial positions - see the + black outlined shapes as examples of this in action. + +Unstructured Meshes (the new world) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +A mesh is made up of different types of **element**: + +.. list-table:: + :widths: 15, 15, 70 + + * - 0D + - ``node`` + - The 'core' of the mesh. A point position in space, constructed from + 2 or 3 coordinates (2D or 3D space). + * - 1D + - ``edge`` + - Constructed by connecting 2 nodes. + * - 2D + - ``face`` + - Constructed by connecting 3 or more nodes. + * - 3D + - ``volume`` + - Constructed by connecting 4 or more nodes (which must each have 3 + coordinates - 3D space). + +Every node in the mesh is defined by indexing the 1-dimensional X and Y (and +optionally Z) coordinate arrays (the ``node_coordinates``) - e.g. +``(x[3], y[3])`` gives the position of the fourth node. Note that this means +each node has its own coordinates, independent of every other node. + +Any higher dimensional element - an edge/face/volume - is described by a +sequence of the indices of the nodes that make up that element. E.g. a +triangular face made from connecting the first, third and fourth nodes: +``[0, 2, 3]``. These 1D sequences combine into a 2D array enumerating **all** +the elements of that type - edge/face/volume - called a **connectivity**. +E.g. we could make a mesh of 4 nodes, with 2 triangles described using this +``face_node_connectivity``: ``[[0, 2, 3], [3, 2, 1]]`` (note the shared nodes). + +.. note:: More on Connectivities: + + * The element type described by a connectivity is known as its + **location**; ``edge`` in ``edge_node_connectivity``. + * According to the UGRID conventions, the nodes in a face should be + listed in "anti-clockwise order from above". + * Connectivities also exist to connect the higher dimensional elements, + e.g. ``face_edge_connectivity``. These are optional conveniences to + speed up certain operations and will not be discussed here. + +.. important:: + + **Meshes are unstructured**. The mesh elements - represented in the + coordinate and connectivity arrays detailed above - are enumerated + along a single **unstructured dimension**. An element's position along + this dimension has nothing to do with its spatial position. + +A data variable associated with a mesh has a **location** of either ``node``, +``edge``, ``face`` or ``volume``. The data is stored in a 1D array with one +datum per element, matched to its element by matching the datum index with the +coordinate or connectivity index along the **unstructured dimension**. So for +an example data array called ``foo``: +``foo[3]`` would be at position ``(x[3], y[3])`` if it were node-located, or at +``faces[3]`` if it were face-located. :numref:`data_ugrid_mesh` visualises an +example of what is described above. + +.. _data_ugrid_mesh: +.. figure:: images/data_ugrid_mesh.svg + :alt: Diagram of how data is represented on an unstructured mesh + :align: right + :width: 1280 + + Data on an unstructured mesh + + 1D coordinate arrays (pink circles) describe node positions in space (pink + crosses). A 2D connectivity array (blue circles) describes faces by + connecting four nodes - by referencing their indices - into a face outline + (blue outlines on the map). Data from the 1D data array (orange circles) + can be assigned to these node locations (orange diamonds) or face locations + (orange quads) by matching the indices in the data array to the indices in + the coordinate arrays (for nodes) or connectivity array (for faces). See + the black outlined shapes as examples of index matching in action, and the + black stippled shapes to demonstrate that relative array position confers + no relative spatial information. + +---- + +The mesh model also supports edges/faces/volumes having associated 'centre' +coordinates - to allow point data to be assigned to these elements. 'Centre' is +just a convenience term - the points can exist anywhere within their respective +elements. See :numref:`ugrid_element_centres` for a visualised example. + +.. _ugrid_element_centres: +.. figure:: images/ugrid_element_centres.svg + :alt: Diagram demonstrating mesh face-centred data. + :align: right + :width: 1280 + + Data can be assigned to mesh edge/face/volume 'centres' + + 1D *node* coordinate arrays (pink circles) describe node positions in + space (pink crosses). A 2D connectivity array (blue circles) describes + faces by connecting four nodes into a face outline (blue outlines on the + map). Further 1D *face* coordinate arrays (pink circles) describe a + 'centre' point position (pink stars) for each face enumerated in the + connectivity array. + +Mesh Flexibility +++++++++++++++++ +Above we have seen how one could replicate data on a structured grid using +a mesh instead. But the utility of a mesh is the extra flexibility it offers. +Here are the main examples: + +Every node is completely independent - every one can have unique X andY (and Z) coordinate values. See :numref:`ugrid_node_independence`. + +.. _ugrid_node_independence: +.. figure:: images/ugrid_node_independence.svg + :alt: Diagram demonstrating the independence of each mesh node + :align: right + :width: 300 + + Every mesh node is completely independent + + The same array shape and structure used to describe the node positions + (pink crosses) in a regular grid (left-hand maps) is equally able to + describe **any** position for these nodes (e.g. the right-hand maps), + simply by changing the array values. The quadrilateral faces (blue + outlines) can therefore be given any quadrilateral shape by re-positioning + their constituent nodes. + +Faces and volumes can have variable node counts, i.e. different numbers of +sides. This is achieved by masking the unused 'slots' in the connectivity +array. See :numref:`ugrid_variable_faces`. + +.. _ugrid_variable_faces: +.. figure:: images/ugrid_variable_faces.svg + :alt: Diagram demonstrating mesh faces with variable node counts + :align: right + :width: 300 + + Mesh faces can have different node counts (using masking) + + The 2D connectivity array (blue circles) describes faces by connecting + nodes (pink crosses) to make up a face (blue outlines). The faces can use + different numbers of nodes by shaping the connectivity array to accommodate + the face with the most nodes, then masking unused node 'slots' + (black circles) for faces with fewer nodes than the maximum. + +Data can be assigned to lines (edges) just as easily as points (nodes) or +areas (faces). See :numref:`ugrid_edge_data`. + +.. _ugrid_edge_data: +.. figure:: images/ugrid_edge_data.svg + :alt: Diagram demonstrating data assigned to mesh edges + :align: right + :width: 300 + + Data can be assigned to mesh edges + + The 2D connectivity array (blue circles) describes edges by connecting 2 + nodes (pink crosses) to make up an edge (blue lines). Data can be assigned + to the edges (orange lines) by matching the indices of the 1D data array + (not shown) to the indices in the connectivity array. + +.. _ugrid implications: + +What does this mean? +-------------------- +Meshes can represent much more varied spatial arrangements +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +The highly specific way of recording position (geometry) and shape +(topology) allows meshes to represent essentially **any** spatial arrangement +of data. There are therefore many new applications that aren't possible using a +structured grid, including: + +* `The UK Met Office's LFRic cubed-sphere `_ +* `Oceanic model outputs `_ + +.. todo: + a third example! + +Mesh 'payload' is much larger than with structured grids +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Coordinates are recorded per-node, and connectivities are recorded per-element. +This is opposed to a structured grid, where a single coordinate value is shared +by every data point/area along that line. + +For example: representing the surface of a cubed-sphere using a mesh leads to +coordinates and connectivities being **~8 times larger than the data itself**, +as opposed to a small fraction of the data size when dividing a spherical +surface using a structured grid of longitudes and latitudes. + +This further increases the emphasis on lazy loading and processing of data +using packages such as Dask. + +.. note:: + + The large, 1D data arrays associated with meshes are a very different + shape to what Iris users and developers are used to. It is suspected + that optimal performance will need new chunking strategies, but at time + of writing (``Jan 2022``) experience is still limited. + +.. todo: + Revisit when we have more information. + +Spatial operations on mesh data are more complex +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Detail: :doc:`operations` + +Indexing a mesh data array cannot be used for: + +#. Region selection +#. Neighbour identification + +This is because - unlike with a structured data array - relative position in +a mesh's 1-dimensional data arrays has no relation to relative position in +space. We must instead perform specialised operations using the information in +the mesh's connectivities, or by translating the mesh into a format designed +for mesh analysis such as VTK. + +Such calculations can still be optimised to avoid them slowing workflows, but +the important take-away here is that **adaptation is needed when working mesh +data**. + + +How Iris Represents This +======================== + +.. + Include API links to the various classes + + Include Cube/Mesh printout(s) + +.. seealso:: + + Remember this is a prose summary. Precise documentation is at: + :mod:`iris.experimental.ugrid`. + +.. note:: + + At time of writing (``Jan 2022``), neither 3D meshes nor 3D elements + (volumes) are supported. + +The Basics +---------- +The Iris :class:`~iris.cube.Cube` has several new members: + +* | :attr:`~iris.cube.Cube.mesh` + | The :class:`iris.experimental.ugrid.Mesh` that describes the + :class:`~iris.cube.Cube`\'s horizontal geography. +* | :attr:`~iris.cube.Cube.location` + | ``node``/``edge``/``face`` - the mesh element type with which this + :class:`~iris.cube.Cube`\'s :attr:`~iris.cube.Cube.data` is associated. +* | :meth:`~iris.cube.Cube.mesh_dim` + | The :class:`~iris.cube.Cube`\'s **unstructured dimension** - the one that + indexes over the horizontal :attr:`~iris.cube.Cube.data` positions. + +These members will all be ``None`` for a :class:`~iris.cube.Cube` with no +associated :class:`~iris.experimental.ugrid.Mesh`. + +This :class:`~iris.cube.Cube`\'s unstructured dimension has multiple attached +:class:`iris.experimental.ugrid.MeshCoord`\s (one for each axis e.g. +``x``/``y``), which can be used to infer the points and bounds of any index on +the :class:`~iris.cube.Cube`\'s unstructured dimension. + +.. testsetup:: ugrid_summaries + + import numpy as np + + from iris.coords import AuxCoord, DimCoord + from iris.cube import Cube + from iris.experimental.ugrid import Connectivity, Mesh + + node_x = AuxCoord( + points=[0.0, 5.0, 0.0, 5.0, 8.0], + standard_name="longitude", + units="degrees_east", + ) + node_y = AuxCoord( + points=[3.0, 3.0, 0.0, 0.0, 0.0], + standard_name="latitude", + units="degrees_north", + ) + + edge_node_c = Connectivity( + indices=[[0, 1], [0, 2], [1, 3], [1, 4], [2, 3], [3, 4]], + cf_role="edge_node_connectivity", + ) + + face_indices = np.ma.masked_equal([[0, 1, 3, 2], [1, 4, 3, 999]], 999) + face_node_c = Connectivity( + indices=face_indices, cf_role="face_node_connectivity" + ) + + def centre_coords(conn): + indexing = np.ma.filled(conn.indices, 0) + x, y = [ + AuxCoord( + node_coord.points[indexing].mean(axis=conn.connected_axis), + node_coord.standard_name, + units=node_coord.units, + ) + for node_coord in (node_x, node_y) + ] + return [(x, "x"), (y, "y")] + + my_mesh = Mesh( + long_name="my_mesh", + topology_dimension=2, + node_coords_and_axes=[(node_x, "x"), (node_y, "y")], + connectivities=[edge_node_c, face_node_c], + edge_coords_and_axes=centre_coords(edge_node_c), + face_coords_and_axes=centre_coords(face_node_c), + ) + + vertical_levels = DimCoord([0, 1, 2], "height") + + def location_cube(conn): + location = conn.location + mesh_coord_x, mesh_coord_y = my_mesh.to_MeshCoords(location) + data_shape = (conn.shape[conn.location_axis], len(vertical_levels.points)) + data_array = np.arange(np.prod(data_shape)).reshape(data_shape) + + return Cube( + data=data_array, + long_name=f"{location}_data", + units="K", + dim_coords_and_dims=[(vertical_levels, 1)], + aux_coords_and_dims=[(mesh_coord_x, 0), (mesh_coord_y, 0)], + ) + + edge_cube = location_cube(edge_node_c) + face_cube = location_cube(face_node_c) + +.. doctest:: ugrid_summaries + + >>> print(edge_cube) + edge_data / (K) (-- : 6; height: 3) + Dimension coordinates: + height - x + Mesh coordinates: + latitude x - + longitude x - + + >>> print(edge_cube.location) + edge + + >>> print(edge_cube.mesh_dim()) + 0 + + >>> print(edge_cube.mesh.summary(shorten=True)) + + +The Detail +---------- +How UGRID information is stored +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +* | :class:`iris.experimental.ugrid.Mesh` + | Contains all information about the mesh. + | Includes: + + * | :attr:`~iris.experimental.ugrid.Mesh.topology_dimension` + | The maximum dimensionality of shape (1D=edge, 2D=face) supported + by this :class:`~iris.experimental.ugrid.Mesh`. Determines which + :class:`~iris.experimental.ugrid.Connectivity`\s are required/optional + (see below). + + * 1-3 collections of :class:`iris.coords.AuxCoord`\s: + + * | **Required**: :attr:`~iris.experimental.ugrid.Mesh.node_coords` + | The nodes that are the basis for the mesh. + * | Optional: :attr:`~iris.experimental.ugrid.Mesh.edge_coords`, + :attr:`~iris.experimental.ugrid.Mesh.face_coords` + | For indicating the 'centres' of the edges/faces. + | **NOTE:** generating a :class:`~iris.experimental.ugrid.MeshCoord` from + a :class:`~iris.experimental.ugrid.Mesh` currently (``Jan 2022``) + requires centre coordinates for the given ``location``; to be rectified + in future. + + * 1 or more :class:`iris.experimental.ugrid.Connectivity`\s: + + * | **Required for 1D (edge) elements**: + :attr:`~iris.experimental.ugrid.Mesh.edge_node_connectivity` + | Define the edges by connecting nodes. + * | **Required for 2D (face) elements**: + :attr:`~iris.experimental.ugrid.Mesh.face_node_connectivity` + | Define the faces by connecting nodes. + * Optional: any other connectivity type. See + :attr:`iris.experimental.ugrid.mesh.Connectivity.UGRID_CF_ROLES` for the + full list of types. + +.. doctest:: ugrid_summaries + + >>> print(edge_cube.mesh) + Mesh : 'my_mesh' + topology_dimension: 2 + node + node_dimension: 'Mesh2d_node' + node coordinates + + + edge + edge_dimension: 'Mesh2d_edge' + edge_node_connectivity: + edge coordinates + + + face + face_dimension: 'Mesh2d_face' + face_node_connectivity: + face coordinates + + + long_name: 'my_mesh' + +* | :class:`iris.experimental.ugrid.MeshCoord` + | Described in detail in `MeshCoords`_. + | Stores the following information: + + * | :attr:`~iris.experimental.ugrid.MeshCoord.mesh` + | The :class:`~iris.experimental.ugrid.Mesh` associated with this + :class:`~iris.experimental.ugrid.MeshCoord`. This determines the + :attr:`~iris.cube.Cube.mesh` attribute of any :class:`~iris.cube.Cube` + this :class:`~iris.experimental.ugrid.MeshCoord` is attached to (see + `The Basics`_) + + * | :attr:`~iris.experimental.ugrid.MeshCoord.location` + | ``node``/``edge``/``face`` - the element detailed by this + :class:`~iris.experimental.ugrid.MeshCoord`. This determines the + :attr:`~iris.cube.Cube.location` attribute of any + :class:`~iris.cube.Cube` this + :class:`~iris.experimental.ugrid.MeshCoord` is attached to (see + `The Basics`_). + +.. _ugrid MeshCoords: + +MeshCoords +~~~~~~~~~~ +Links a :class:`~iris.cube.Cube` to a :class:`~iris.experimental.ugrid.Mesh` by +attaching to the :class:`~iris.cube.Cube`\'s unstructured dimension, in the +same way that all :class:`~iris.coords.Coord`\s attach to +:class:`~iris.cube.Cube` dimensions. This allows a single +:class:`~iris.cube.Cube` to have a combination of unstructured and structured +dimensions (e.g. horizontal mesh plus vertical levels and a time series), +using the same logic for every dimension. + +:class:`~iris.experimental.ugrid.MeshCoord`\s are instantiated using a given +:class:`~iris.experimental.ugrid.Mesh`, ``location`` +("node"/"edge"/"face") and ``axis``. The process interprets the +:class:`~iris.experimental.ugrid.Mesh`\'s +:attr:`~iris.experimental.ugrid.Mesh.node_coords` and if appropriate the +:attr:`~iris.experimental.ugrid.Mesh.edge_node_connectivity`/ +:attr:`~iris.experimental.ugrid.Mesh.face_node_connectivity` and +:attr:`~iris.experimental.ugrid.Mesh.edge_coords`/ +:attr:`~iris.experimental.ugrid.Mesh.face_coords` +to produce a :class:`~iris.coords.Coord` +:attr:`~iris.coords.Coord.points` and :attr:`~iris.coords.Coord.bounds` +representation of all the :class:`~iris.experimental.ugrid.Mesh`\'s +nodes/edges/faces for the given axis. + +The method :meth:`iris.experimental.ugrid.Mesh.to_MeshCoords` is available to +create a :class:`~iris.experimental.ugrid.MeshCoord` for +every axis represented by that :class:`~iris.experimental.ugrid.Mesh`, +given only the ``location`` argument + +.. doctest:: ugrid_summaries + + >>> for coord in edge_cube.coords(mesh_coords=True): + ... print(coord) + MeshCoord : latitude / (degrees_north) + mesh: + location: 'edge' + points: [3. , 1.5, 1.5, 1.5, 0. , 0. ] + bounds: [ + [3., 3.], + [3., 0.], + [3., 0.], + [3., 0.], + [0., 0.], + [0., 0.]] + shape: (6,) bounds(6, 2) + dtype: float64 + standard_name: 'latitude' + axis: 'y' + MeshCoord : longitude / (degrees_east) + mesh: + location: 'edge' + points: [2.5, 0. , 5. , 6.5, 2.5, 6.5] + bounds: [ + [0., 5.], + [0., 0.], + [5., 5.], + [5., 8.], + [0., 5.], + [5., 8.]] + shape: (6,) bounds(6, 2) + dtype: float64 + standard_name: 'longitude' + axis: 'x' + + +__ CF-UGRID_ \ No newline at end of file diff --git a/docs/src/further_topics/ugrid/images/data_structured_grid.svg b/docs/src/further_topics/ugrid/images/data_structured_grid.svg new file mode 100644 index 0000000000..2f3a1ce342 --- /dev/null +++ b/docs/src/further_topics/ugrid/images/data_structured_grid.svg @@ -0,0 +1 @@ +23, 28-19,-21101525-5-15-20-30xyCoordinate ArraysxyCoordinate Arrays23, 28-19, -21xyBounds Arraysderive point locationsassign data using dimensional indices,position in array == relative spatial positionderive area locations & shapesPoint DataArea DataData Array(bounded coordsalways have points too)my_variable* x+yare not lons+lats, just a demonstration! \ No newline at end of file diff --git a/docs/src/further_topics/ugrid/images/data_ugrid_mesh.svg b/docs/src/further_topics/ugrid/images/data_ugrid_mesh.svg new file mode 100644 index 0000000000..ab7302346b --- /dev/null +++ b/docs/src/further_topics/ugrid/images/data_ugrid_mesh.svg @@ -0,0 +1 @@ +5, 7, 8, 14`xy1212`node_coordinates`every node has its own x + y coordinatesderive node locations1515xy`node_coordinates`[5][7][8][14]construct faces by connecting nodesderive ‘corner’ node locationsassign data using 1D indexing,position in array unrelated to spatial positionmatch indices with facesmatch indices with nodesNode DataFace Data12Data Arraymy_variable12 ×4`face_node_connectivity`face_nodes \ No newline at end of file diff --git a/docs/src/further_topics/ugrid/images/geovistalogo.svg b/docs/src/further_topics/ugrid/images/geovistalogo.svg new file mode 100644 index 0000000000..4c68f0ee3f --- /dev/null +++ b/docs/src/further_topics/ugrid/images/geovistalogo.svg @@ -0,0 +1,573 @@ + + + + + + + + + + + + + + + + + + + + + + + + Cartographic rendering and mesh analytics powered by PyVista. + GeoVista + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + GeoVista + + \ No newline at end of file diff --git a/docs/src/further_topics/ugrid/images/iris-esmf-regrid.svg b/docs/src/further_topics/ugrid/images/iris-esmf-regrid.svg new file mode 100644 index 0000000000..e70a9386a7 --- /dev/null +++ b/docs/src/further_topics/ugrid/images/iris-esmf-regrid.svg @@ -0,0 +1,93 @@ + + + + + + + + + + image/svg+xml + + + + + + + + + Iris + + diff --git a/docs/src/further_topics/ugrid/images/plotting_basic.png b/docs/src/further_topics/ugrid/images/plotting_basic.png new file mode 100644 index 0000000000..ba2b0b3329 Binary files /dev/null and b/docs/src/further_topics/ugrid/images/plotting_basic.png differ diff --git a/docs/src/further_topics/ugrid/images/plotting_global.png b/docs/src/further_topics/ugrid/images/plotting_global.png new file mode 100644 index 0000000000..62fb56d974 Binary files /dev/null and b/docs/src/further_topics/ugrid/images/plotting_global.png differ diff --git a/docs/src/further_topics/ugrid/images/ugrid_edge_data.svg b/docs/src/further_topics/ugrid/images/ugrid_edge_data.svg new file mode 100644 index 0000000000..374ef57388 --- /dev/null +++ b/docs/src/further_topics/ugrid/images/ugrid_edge_data.svg @@ -0,0 +1 @@ +`edge_node_connectivity`12 ×2 \ No newline at end of file diff --git a/docs/src/further_topics/ugrid/images/ugrid_element_centres.svg b/docs/src/further_topics/ugrid/images/ugrid_element_centres.svg new file mode 100644 index 0000000000..13b885d600 --- /dev/null +++ b/docs/src/further_topics/ugrid/images/ugrid_element_centres.svg @@ -0,0 +1 @@ +`face_node_connectivity`xy`node_coordinates`xy`face_coordinates`151512 ×41212`face_coordinates``node_coordinates` \ No newline at end of file diff --git a/docs/src/further_topics/ugrid/images/ugrid_node_independence.svg b/docs/src/further_topics/ugrid/images/ugrid_node_independence.svg new file mode 100644 index 0000000000..ba72c42ffa --- /dev/null +++ b/docs/src/further_topics/ugrid/images/ugrid_node_independence.svg @@ -0,0 +1 @@ +` \ No newline at end of file diff --git a/docs/src/further_topics/ugrid/images/ugrid_variable_faces.svg b/docs/src/further_topics/ugrid/images/ugrid_variable_faces.svg new file mode 100644 index 0000000000..378978abc3 --- /dev/null +++ b/docs/src/further_topics/ugrid/images/ugrid_variable_faces.svg @@ -0,0 +1 @@ +`face_node_connectivity`12 ×6 \ No newline at end of file diff --git a/docs/src/further_topics/ugrid/index.rst b/docs/src/further_topics/ugrid/index.rst new file mode 100644 index 0000000000..81ba24428a --- /dev/null +++ b/docs/src/further_topics/ugrid/index.rst @@ -0,0 +1,54 @@ +.. include:: ../../common_links.inc + +.. _ugrid: + +Mesh Support +************ + +Iris includes specialised handling of mesh-located data (as opposed to +grid-located data). Iris and its :ref:`partner packages ` are +designed to make working with mesh-located data as simple as possible, with new +capabilities being added all the time. More detail is in this section and in +the :mod:`iris.experimental.ugrid` API documentation. + +This mesh support is based on the `CF-UGRID Conventions`__; UGRID-conformant +meshes + data can be loaded from a file into Iris' data model, and meshes + +data represented in Iris' data model can be saved as a UGRID-conformant file. + +---- + +Meshes are different + Mesh-located data is fundamentally different to grid-located data. + Many of Iris' existing operations need adapting before they can work with + mesh-located data, and in some cases entirely new concepts are needed. + **Read the detail in these pages before jumping into your own code.** +Iris' mesh support is experimental + This is a rapidly evolving part of the codebase at time of writing + (``Jan 2022``), as we continually expand the operations that work with mesh + data. **Be prepared for breaking changes even in minor releases.** +:ref:`Get involved! ` + We know meshes are an exciting new area for much of Earth science, so we hope + there are a lot of you with new files/ideas/wishlists, and we'd love to hear + more 🙂. + +---- + +Read on to find out more... + +* :doc:`data_model` - learn why the mesh experience is so different. +* :doc:`partner_packages` - meet some optional dependencies that provide powerful mesh operations. +* :doc:`operations` - experience how your workflows will look when written for mesh data. + +.. + Need an actual TOC to get Sphinx working properly, but have hidden it in + favour of the custom bullets above. + +.. toctree:: + :hidden: + :maxdepth: 1 + + data_model + partner_packages + operations + +__ CF-UGRID_ diff --git a/docs/src/further_topics/ugrid/operations.rst b/docs/src/further_topics/ugrid/operations.rst new file mode 100644 index 0000000000..84e817bf48 --- /dev/null +++ b/docs/src/further_topics/ugrid/operations.rst @@ -0,0 +1,995 @@ +.. _ugrid operations: + +Working with Mesh Data +********************** + +.. note:: Several of the operations below rely on the optional dependencies + mentioned in :doc:`partner_packages`. + +Operations Summary +------------------ +.. list-table:: + :align: left + :widths: 35, 75 + + * - `Making a Mesh`_ + - |tagline: making a mesh| + * - `Making a Cube`_ + - |tagline: making a cube| + * - `Save`_ + - |tagline: save| + * - `Load`_ + - |tagline: load| + * - `Plotting`_ + - |tagline: plotting| + * - `Region Extraction`_ + - |tagline: region extraction| + * - `Regridding`_ + - |tagline: regridding| + * - `Equality`_ + - |tagline: equality| + * - `Combining Cubes`_ + - |tagline: combining cubes| + * - `Arithmetic`_ + - |tagline: arithmetic| + +.. + Below: use demo code over prose wherever workable. Headings aren't an + exhaustive list (can you think of any other popular operations?). + +Making a Mesh +------------- +.. |tagline: making a mesh| replace:: |new| + +.. rubric:: |tagline: making a mesh| + +**Already have a file?** Consider skipping to `Load`_. + +Creating Iris objects from scratch is a highly useful skill for testing code +and improving understanding of how Iris works. This knowledge will likely prove +particularly useful when converting data into the Iris mesh data model from +structured formats and non-UGRID mesh formats. + +The objects created in this example will be used where possible in the +subsequent example operations on this page. + +.. dropdown:: :opticon:`code` + + .. doctest:: ugrid_operations + + >>> import numpy as np + + >>> from iris.coords import AuxCoord + >>> from iris.experimental.ugrid import Connectivity, Mesh + + # Going to create the following mesh + # (node indices are shown to aid understanding): + # + # 0----1 + # | |\ + # | + |+\ + # 2----3--4 + + >>> node_x = AuxCoord( + ... points=[0.0, 5.0, 0.0, 5.0, 8.0], + ... standard_name="longitude", + ... units="degrees_east", + ... long_name="node_x_coordinates", + ... ) + >>> node_y = AuxCoord(points=[3.0, 3.0, 0.0, 0.0, 0.0], standard_name="latitude") + + >>> face_x = AuxCoord([2.0, 6.0], "longitude") + >>> face_y = AuxCoord([1.0, 1.0], "latitude") + + >>> edge_node_c = Connectivity( + ... indices=[[0, 1], [0, 2], [1, 3], [1, 4], [2, 3], [3, 4]], + ... cf_role="edge_node_connectivity", + ... attributes={"demo": "Supports every standard CF property"}, + ... ) + + # Create some dead-centre edge coordinates. + >>> edge_x, edge_y = [ + ... AuxCoord( + ... node_coord.points[edge_node_c.indices_by_location()].mean(axis=1), + ... node_coord.standard_name, + ... ) + ... for node_coord in (node_x, node_y) + ... ] + + >>> face_indices = np.ma.masked_equal([[0, 1, 3, 2], [1, 4, 3, 999]], 999) + >>> face_node_c = Connectivity( + ... indices=face_indices, cf_role="face_node_connectivity" + ... ) + + >>> my_mesh = Mesh( + ... long_name="my_mesh", + ... topology_dimension=2, # Supports 2D (face) elements. + ... node_coords_and_axes=[(node_x, "x"), (node_y, "y")], + ... connectivities=[edge_node_c, face_node_c], + ... edge_coords_and_axes=[(edge_x, "x"), (edge_y, "y")], + ... face_coords_and_axes=[(face_x, "x"), (face_y, "y")], + ... ) + + >>> print(my_mesh) + Mesh : 'my_mesh' + topology_dimension: 2 + node + node_dimension: 'Mesh2d_node' + node coordinates + + + edge + edge_dimension: 'Mesh2d_edge' + edge_node_connectivity: + edge coordinates + + + face + face_dimension: 'Mesh2d_face' + face_node_connectivity: + face coordinates + + + long_name: 'my_mesh' + + +.. _making a cube: + +Making a Cube (with a Mesh) +--------------------------- +.. |tagline: making a cube| replace:: |unchanged| + +.. rubric:: |tagline: making a cube| + +Creating a :class:`~iris.cube.Cube` is unchanged; the +:class:`~iris.experimental.ugrid.Mesh` is linked via a +:class:`~iris.experimental.ugrid.MeshCoord` (see :ref:`ugrid MeshCoords`): + +.. dropdown:: :opticon:`code` + + .. doctest:: ugrid_operations + + >>> import numpy as np + + >>> from iris.coords import DimCoord + >>> from iris.cube import Cube, CubeList + + >>> vertical_levels = DimCoord([0, 1, 2], "height") + + >>> my_cubelist = CubeList() + >>> for conn in (edge_node_c, face_node_c): + ... location = conn.location + ... mesh_coord_x, mesh_coord_y = my_mesh.to_MeshCoords(location) + ... data_shape = (len(conn.indices_by_location()), len(vertical_levels.points)) + ... data_array = np.arange(np.prod(data_shape)).reshape(data_shape) + ... + ... my_cubelist.append( + ... Cube( + ... data=data_array, + ... long_name=f"{location}_data", + ... units="K", + ... dim_coords_and_dims=[(vertical_levels, 1)], + ... aux_coords_and_dims=[(mesh_coord_x, 0), (mesh_coord_y, 0)], + ... ) + ... ) + + >>> print(my_cubelist) + 0: edge_data / (K) (-- : 6; height: 3) + 1: face_data / (K) (-- : 2; height: 3) + + >>> for cube in my_cubelist: + ... print(f"{cube.name()}: {cube.mesh.name()}, {cube.location}") + edge_data: my_mesh, edge + face_data: my_mesh, face + + >>> print(my_cubelist.extract_cube("edge_data")) + edge_data / (K) (-- : 6; height: 3) + Dimension coordinates: + height - x + Mesh coordinates: + latitude x - + longitude x - + + +Save +---- +.. |tagline: save| replace:: |unchanged| + +.. rubric:: |tagline: save| + +.. note:: UGRID saving support is limited to the NetCDF file format. + +The Iris saving process automatically detects if the :class:`~iris.cube.Cube` +has an associated :class:`~iris.experimental.ugrid.Mesh` and automatically +saves the file in a UGRID-conformant format: + +.. dropdown:: :opticon:`code` + + .. doctest:: ugrid_operations + + >>> from subprocess import run + + >>> from iris import save + + >>> cubelist_path = "my_cubelist.nc" + >>> save(my_cubelist, cubelist_path) + + >>> ncdump_result = run(["ncdump", "-h", cubelist_path], capture_output=True) + >>> print(ncdump_result.stdout.decode().replace("\t", " ")) + netcdf my_cubelist { + dimensions: + Mesh2d_node = 5 ; + Mesh2d_edge = 6 ; + Mesh2d_face = 2 ; + height = 3 ; + my_mesh_face_N_nodes = 4 ; + my_mesh_edge_N_nodes = 2 ; + variables: + int my_mesh ; + my_mesh:cf_role = "mesh_topology" ; + my_mesh:topology_dimension = 2 ; + my_mesh:long_name = "my_mesh" ; + my_mesh:node_coordinates = "longitude latitude" ; + my_mesh:edge_coordinates = "longitude_0 latitude_0" ; + my_mesh:face_coordinates = "longitude_1 latitude_1" ; + my_mesh:face_node_connectivity = "mesh2d_face" ; + my_mesh:edge_node_connectivity = "mesh2d_edge" ; + double longitude(Mesh2d_node) ; + longitude:units = "degrees_east" ; + longitude:standard_name = "longitude" ; + longitude:long_name = "node_x_coordinates" ; + double latitude(Mesh2d_node) ; + latitude:standard_name = "latitude" ; + double longitude_0(Mesh2d_edge) ; + longitude_0:standard_name = "longitude" ; + double latitude_0(Mesh2d_edge) ; + latitude_0:standard_name = "latitude" ; + double longitude_1(Mesh2d_face) ; + longitude_1:standard_name = "longitude" ; + double latitude_1(Mesh2d_face) ; + latitude_1:standard_name = "latitude" ; + int64 mesh2d_face(Mesh2d_face, my_mesh_face_N_nodes) ; + mesh2d_face:_FillValue = -1LL ; + mesh2d_face:cf_role = "face_node_connectivity" ; + mesh2d_face:start_index = 0LL ; + int64 mesh2d_edge(Mesh2d_edge, my_mesh_edge_N_nodes) ; + mesh2d_edge:demo = "Supports every standard CF property" ; + mesh2d_edge:cf_role = "edge_node_connectivity" ; + mesh2d_edge:start_index = 0LL ; + int64 edge_data(Mesh2d_edge, height) ; + edge_data:long_name = "edge_data" ; + edge_data:units = "K" ; + edge_data:mesh = "my_mesh" ; + edge_data:location = "edge" ; + int64 height(height) ; + height:standard_name = "height" ; + int64 face_data(Mesh2d_face, height) ; + face_data:long_name = "face_data" ; + face_data:units = "K" ; + face_data:mesh = "my_mesh" ; + face_data:location = "face" ; + + // global attributes: + :Conventions = "CF-1.7" ; + } + + +The :func:`iris.experimental.ugrid.save_mesh` function allows +:class:`~iris.experimental.ugrid.Mesh`\es to be saved to file without +associated :class:`~iris.cube.Cube`\s: + +.. dropdown:: :opticon:`code` + + .. doctest:: ugrid_operations + + >>> from subprocess import run + + >>> from iris.experimental.ugrid import save_mesh + + >>> mesh_path = "my_mesh.nc" + >>> save_mesh(my_mesh, mesh_path) + + >>> ncdump_result = run(["ncdump", "-h", mesh_path], capture_output=True) + >>> print(ncdump_result.stdout.decode().replace("\t", " ")) + netcdf my_mesh { + dimensions: + Mesh2d_node = 5 ; + Mesh2d_edge = 6 ; + Mesh2d_face = 2 ; + my_mesh_face_N_nodes = 4 ; + my_mesh_edge_N_nodes = 2 ; + variables: + int my_mesh ; + my_mesh:cf_role = "mesh_topology" ; + my_mesh:topology_dimension = 2 ; + my_mesh:long_name = "my_mesh" ; + my_mesh:node_coordinates = "longitude latitude" ; + my_mesh:edge_coordinates = "longitude_0 latitude_0" ; + my_mesh:face_coordinates = "longitude_1 latitude_1" ; + my_mesh:face_node_connectivity = "mesh2d_face" ; + my_mesh:edge_node_connectivity = "mesh2d_edge" ; + double longitude(Mesh2d_node) ; + longitude:units = "degrees_east" ; + longitude:standard_name = "longitude" ; + longitude:long_name = "node_x_coordinates" ; + double latitude(Mesh2d_node) ; + latitude:standard_name = "latitude" ; + double longitude_0(Mesh2d_edge) ; + longitude_0:standard_name = "longitude" ; + double latitude_0(Mesh2d_edge) ; + latitude_0:standard_name = "latitude" ; + double longitude_1(Mesh2d_face) ; + longitude_1:standard_name = "longitude" ; + double latitude_1(Mesh2d_face) ; + latitude_1:standard_name = "latitude" ; + int64 mesh2d_face(Mesh2d_face, my_mesh_face_N_nodes) ; + mesh2d_face:_FillValue = -1LL ; + mesh2d_face:cf_role = "face_node_connectivity" ; + mesh2d_face:start_index = 0LL ; + int64 mesh2d_edge(Mesh2d_edge, my_mesh_edge_N_nodes) ; + mesh2d_edge:demo = "Supports every standard CF property" ; + mesh2d_edge:cf_role = "edge_node_connectivity" ; + mesh2d_edge:start_index = 0LL ; + + // global attributes: + :Conventions = "CF-1.7" ; + } + + +Load +---- +.. |tagline: load| replace:: |different| - UGRID parsing is opt-in + +.. rubric:: |tagline: load| + +.. note:: UGRID loading support is limited to the NetCDF file format. + +While Iris' UGRID support remains :mod:`~iris.experimental`, parsing UGRID when +loading a file remains **optional**. To load UGRID data from a file into the +Iris mesh data model, use the +:const:`iris.experimental.ugrid.PARSE_UGRID_ON_LOAD` context manager: + +.. dropdown:: :opticon:`code` + + .. doctest:: ugrid_operations + + >>> from iris import load + >>> from iris.experimental.ugrid import PARSE_UGRID_ON_LOAD + + >>> with PARSE_UGRID_ON_LOAD.context(): + ... loaded_cubelist = load(cubelist_path) + + # Sort CubeList to ensure consistent result. + >>> loaded_cubelist.sort(key=lambda cube: cube.name()) + >>> print(loaded_cubelist) + 0: edge_data / (K) (-- : 6; height: 3) + 1: face_data / (K) (-- : 2; height: 3) + +All the existing loading functionality still operates on UGRID-compliant +data - :class:`~iris.Constraint`\s, callbacks, :func:`~iris.load_cube` +etcetera: + +.. dropdown:: :opticon:`code` + + .. doctest:: ugrid_operations + + >>> from iris import Constraint, load_cube + + >>> with PARSE_UGRID_ON_LOAD.context(): + ... ground_cubelist = load(cubelist_path, Constraint(height=0)) + ... face_cube = load_cube(cubelist_path, "face_data") + + # Sort CubeList to ensure consistent result. + >>> ground_cubelist.sort(key=lambda cube: cube.name()) + >>> print(ground_cubelist) + 0: edge_data / (K) (-- : 6) + 1: face_data / (K) (-- : 2) + + >>> print(face_cube) + face_data / (K) (-- : 2; height: 3) + Dimension coordinates: + height - x + Mesh coordinates: + latitude x - + longitude x - + Attributes: + Conventions CF-1.7 + +.. note:: + + We recommend caution if constraining on coordinates associated with a + :class:`~iris.experimental.ugrid.Mesh`. An individual coordinate value + might not be shared by any other data points, and using a coordinate range + will demand notably higher performance given the size of the dimension + versus structured grids + (:ref:`see the data model detail `). + +The :func:`iris.experimental.ugrid.load_mesh` and +:func:`~iris.experimental.ugrid.load_meshes` functions allow only +:class:`~iris.experimental.ugrid.Mesh`\es to be loaded from a file without +creating any associated :class:`~iris.cube.Cube`\s: + +.. dropdown:: :opticon:`code` + + .. doctest:: ugrid_operations + + >>> from iris.experimental.ugrid import load_mesh + + >>> with PARSE_UGRID_ON_LOAD.context(): + ... loaded_mesh = load_mesh(cubelist_path) + + >>> print(loaded_mesh) + Mesh : 'my_mesh' + topology_dimension: 2 + node + node_dimension: 'Mesh2d_node' + node coordinates + shape(5,)> + shape(5,)> + edge + edge_dimension: 'Mesh2d_edge' + edge_node_connectivity: shape(6, 2)> + edge coordinates + shape(6,)> + shape(6,)> + face + face_dimension: 'Mesh2d_face' + face_node_connectivity: shape(2, 4)> + face coordinates + shape(2,)> + shape(2,)> + long_name: 'my_mesh' + var_name: 'my_mesh' + +Plotting +-------- +.. |tagline: plotting| replace:: |different| - plot with GeoVista + +.. rubric:: |tagline: plotting| + +The Cartopy-Matplotlib combination is not optimised for displaying the high +number of irregular shapes associated with meshes. Thankfully mesh +visualisation is already popular in many other fields (e.g. CGI, gaming, +SEM microscopy), so there is a wealth of tooling available, which +:ref:`ugrid geovista` harnesses for cartographic plotting. + +GeoVista's default behaviour is to convert lat-lon information into full XYZ +coordinates so the data is visualised on the surface of a 3D globe. The plots +are interactive by default, so it's easy to explore the data in detail. + +2D projections have also been demonstrated in proofs of concept, and will +be added to API in the near future. + +This first example uses GeoVista to plot the ``face_cube`` that we created +earlier: + +.. dropdown:: :opticon:`code` + + .. code-block:: python + + >>> from geovista import GeoPlotter, Transform + >>> from geovista.common import to_xyz + + + # We'll re-use this to plot some real global data later. + >>> def cube_faces_to_polydata(cube): + ... lons, lats = cube.mesh.node_coords + ... face_node = cube.mesh.face_node_connectivity + ... indices = face_node.indices_by_location() + ... + ... mesh = Transform.from_unstructured( + ... lons.points, + ... lats.points, + ... indices, + ... data=cube.data, + ... name=f"{cube.name()} / {cube.units}", + ... start_index=face_node.start_index, + ... ) + ... return mesh + + >>> print(face_cube) + face_data / (K) (-- : 2; height: 3) + Dimension coordinates: + height - x + Mesh coordinates: + latitude x - + longitude x - + Attributes: + Conventions CF-1.7 + + # Convert our mesh+data to a PolyData object. + # Just plotting a single height level. + >>> face_polydata = cube_faces_to_polydata(face_cube[:, 0]) + >>> print(face_polydata) + PolyData (0x7ff4861ff4c0) + N Cells: 2 + N Points: 5 + X Bounds: 9.903e-01, 1.000e+00 + Y Bounds: 0.000e+00, 1.392e-01 + Z Bounds: 6.123e-17, 5.234e-02 + N Arrays: 2 + + # Create the GeoVista plotter and add our mesh+data to it. + >>> my_plotter = GeoPlotter() + >>> my_plotter.add_coastlines(color="black") + >>> my_plotter.add_base_layer(color="grey") + >>> my_plotter.add_mesh(face_polydata) + + # Centre the camera on the data. + >>> camera_region = to_xyz( + ... face_cube.coord("longitude").points, + ... face_cube.coord("latitude").points, + ... radius=3, + ... ) + >>> camera_pos = camera_region.mean(axis=0) + >>> my_plotter.camera.position = camera_pos + + >>> my_plotter.show() + + .. image:: images/plotting_basic.png + :alt: A GeoVista plot of the basic example Mesh. + + This artificial data makes West Africa rather chilly! + +Here's another example using a global cubed-sphere data set: + +.. dropdown:: :opticon:`code` + + .. code-block:: python + + >>> from iris import load_cube + >>> from iris.experimental.ugrid import PARSE_UGRID_ON_LOAD + + # Demonstrating with a global data set. + # You could also download this file from github.com/SciTools/iris-test-data. + >>> from iris.tests import get_data_path + >>> file_path = get_data_path( + ... [ + ... "NetCDF", + ... "unstructured_grid", + ... "lfric_surface_mean.nc", + ... ] + ... ) + >>> with PARSE_UGRID_ON_LOAD.context(): + ... global_cube = load_cube(file_path, "tstar_sea") + >>> print(global_cube) + sea_surface_temperature / (K) (-- : 1; -- : 13824) + Mesh coordinates: + latitude - x + longitude - x + Auxiliary coordinates: + time x - + Cell methods: + mean time (300 s) + mean time_counter + Attributes: + Conventions UGRID + description Created by xios + interval_operation 300 s + interval_write 1 d + name lfric_surface + online_operation average + timeStamp 2020-Feb-07 16:23:14 GMT + title Created by xios + uuid 489bcef5-3d1c-4529-be42-4ab5f8c8497b + + >>> global_polydata = cube_faces_to_polydata(global_cube) + >>> print(global_polydata) + PolyData (0x7f761b536160) + N Cells: 13824 + N Points: 13826 + X Bounds: -1.000e+00, 1.000e+00 + Y Bounds: -1.000e+00, 1.000e+00 + Z Bounds: -1.000e+00, 1.000e+00 + N Arrays: 2 + + >>> my_plotter = GeoPlotter() + >>> my_plotter.add_coastlines() + >>> my_plotter.add_mesh(global_polydata, show_edges=True) + + >>> my_plotter.show() + + .. image:: images/plotting_global.png + :alt: A GeoVista plot of a global sea surface temperature Mesh. + +Region Extraction +----------------- +.. |tagline: region extraction| replace:: |different| - use GeoVista for mesh analysis + +.. rubric:: |tagline: region extraction| + +As described in :doc:`data_model`, indexing for a range along a +:class:`~iris.cube.Cube`\'s :meth:`~iris.cube.Cube.mesh_dim` will not provide +a contiguous region, since **position on the unstructured dimension is +unrelated to spatial position**. This means that subsetted +:class:`~iris.experimental.ugrid.MeshCoord`\s cannot be reliably interpreted +as intended, and subsetting a :class:`~iris.experimental.ugrid.MeshCoord` is +therefore set to return an :class:`~iris.coords.AuxCoord` instead - breaking +the link between :class:`~iris.cube.Cube` and +:class:`~iris.experimental.ugrid.Mesh`: + +.. dropdown:: :opticon:`code` + + .. doctest:: ugrid_operations + + >>> edge_cube = my_cubelist.extract_cube("edge_data") + >>> print(edge_cube) + edge_data / (K) (-- : 6; height: 3) + Dimension coordinates: + height - x + Mesh coordinates: + latitude x - + longitude x - + + # Sub-setted MeshCoords have become AuxCoords. + >>> print(edge_cube[:-1]) + edge_data / (K) (-- : 5; height: 3) + Dimension coordinates: + height - x + Auxiliary coordinates: + latitude x - + longitude x - + +Extracting a region therefore requires extra steps - to determine the spatial +position of the data points before they can be analysed as inside/outside the +selected region. The recommended way to do this is using tools provided by +:ref:`ugrid geovista`, which is optimised for performant mesh analysis. + +This approach centres around using :meth:`geovista.geodesic.BBox.enclosed` to +get the subset of the original mesh that is inside the +:class:`~geovista.geodesic.BBox`. This subset :class:`pyvista.PolyData` object +includes the original indices of each datapoint - the ``vtkOriginalCellIds`` +array, which can be used to index the original :class:`~iris.cube.Cube`. Since +we **know** that this subset :class:`~iris.cube.Cube` represents a regional +mesh, we then reconstruct a :class:`~iris.experimental.ugrid.Mesh` from the +:class:`~iris.cube.Cube`\'s :attr:`~iris.cube.Cube.aux_coords` using +:meth:`iris.experimental.ugrid.Mesh.from_coords`: + +.. + Not using doctest here as want to keep GeoVista as optional dependency. + +.. dropdown:: :opticon:`code` + + .. code-block:: python + + >>> from geovista import Transform + >>> from geovista.geodesic import BBox + >>> from iris import load_cube + >>> from iris.experimental.ugrid import Mesh, PARSE_UGRID_ON_LOAD + + # Need a larger dataset to demonstrate this operation. + # You could also download this file from github.com/SciTools/iris-test-data. + >>> from iris.tests import get_data_path + >>> file_path = get_data_path( + ... [ + ... "NetCDF", + ... "unstructured_grid", + ... "lfric_ngvat_2D_72t_face_half_levels_main_conv_rain.nc", + ... ] + ... ) + + >>> with PARSE_UGRID_ON_LOAD.context(): + ... global_cube = load_cube(file_path, "conv_rain") + >>> print(global_cube) + surface_convective_rainfall_rate / (kg m-2 s-1) (-- : 72; -- : 864) + Mesh coordinates: + latitude - x + longitude - x + Auxiliary coordinates: + time x - + Cell methods: + point time + Attributes: + Conventions UGRID + description Created by xios + interval_operation 300 s + interval_write 300 s + name lfric_ngvat_2D_72t_face_half_levels_main_conv_rain + online_operation instant + timeStamp 2020-Oct-18 21:18:35 GMT + title Created by xios + uuid b3dc0fb4-9828-4663-a5ac-2a5763280159 + + # Convert the Mesh to a GeoVista PolyData object. + >>> lons, lats = global_cube.mesh.node_coords + >>> face_node = global_cube.mesh.face_node_connectivity + >>> indices = face_node.indices_by_location() + >>> global_polydata = Transform.from_unstructured( + ... lons.points, lats.points, indices, start_index=face_node.start_index + ... ) + + # Define a region of 4 corners connected by great circles. + # Specialised sub-classes of BBox are also available e.g. panel/wedge. + >>> region = BBox(lons=[0, 70, 70, 0], lats=[-25, -25, 45, 45]) + # 'Apply' the region to the PolyData object. + >>> region_polydata = region.enclosed(global_polydata, preference="center") + # Get the remaining face indices, to use for indexing the Cube. + >>> indices = region_polydata["vtkOriginalCellIds"] + + >>> print(type(indices)) + + # 101 is smaller than the original 864. + >>> print(len(indices)) + 101 + >>> print(indices[:10]) + [ 6 7 8 9 10 11 18 19 20 21] + + # Use the face indices to subset the global cube. + >>> region_cube = global_cube[:, indices] + + # In this case we **know** the indices correspond to a contiguous + # region, so we will convert the sub-setted Cube back into a + # Cube-with-Mesh. + >>> new_mesh = Mesh.from_coords(*region_cube.coords(dimensions=1)) + >>> new_mesh_coords = new_mesh.to_MeshCoords(global_cube.location) + >>> for coord in new_mesh_coords: + ... region_cube.remove_coord(coord.name()) + ... region_cube.add_aux_coord(coord, 1) + + # A Mesh-Cube with a subset (101) of the original 864 faces. + >>> print(region_cube) + surface_convective_rainfall_rate / (kg m-2 s-1) (-- : 72; -- : 101) + Mesh coordinates: + latitude - x + longitude - x + Auxiliary coordinates: + time x - + Cell methods: + point time + Attributes: + Conventions UGRID + description Created by xios + interval_operation 300 s + interval_write 300 s + name lfric_ngvat_2D_72t_face_half_levels_main_conv_rain + online_operation instant + timeStamp 2020-Oct-18 21:18:35 GMT + title Created by xios + uuid b3dc0fb4-9828-4663-a5ac-2a5763280159 + +Regridding +---------- +.. |tagline: regridding| replace:: |different| - use iris-esmf-regrid for mesh regridders + +.. rubric:: |tagline: regridding| + +Regridding to or from a mesh requires different logic than Iris' existing +regridders, which are designed for structured grids. For this we recommend +ESMF's powerful regridding tools, which integrate with Iris' mesh data model +via the :ref:`ugrid iris-esmf-regrid` package. + +.. todo: inter-sphinx links when available. + +Regridding is achieved via the +:class:`esmf_regrid.experimental.unstructured_scheme.MeshToGridESMFRegridder` +and +:class:`~esmf_regrid.experimental.unstructured_scheme.GridToMeshESMFRegridder` +classes. Regridding from a source :class:`~iris.cube.Cube` to a target +:class:`~iris.cube.Cube` involves initialising and then calling one of these +classes. Initialising is done by passing in the source and target +:class:`~iris.cube.Cube` as arguments. The regridder is then called by passing +the source :class:`~iris.cube.Cube` as an argument. We can demonstrate this +with the +:class:`~esmf_regrid.experimental.unstructured_scheme.MeshToGridESMFRegridder`: + +.. + Not using doctest here as want to keep iris-esmf-regrid as optional dependency. + +.. dropdown:: :opticon:`code` + + .. code-block:: python + + >>> from esmf_regrid.experimental.unstructured_scheme import MeshToGridESMFRegridder + >>> from iris import load, load_cube + >>> from iris.experimental.ugrid import PARSE_UGRID_ON_LOAD + + # You could also download these files from github.com/SciTools/iris-test-data. + >>> from iris.tests import get_data_path + >>> mesh_file = get_data_path( + ... ["NetCDF", "unstructured_grid", "lfric_surface_mean.nc"] + ... ) + >>> grid_file = get_data_path( + ... ["NetCDF", "regrid", "regrid_template_global_latlon.nc"] + ... ) + + # Load a list of cubes defined on the same Mesh. + >>> with PARSE_UGRID_ON_LOAD.context(): + ... mesh_cubes = load(mesh_file) + + # Extract a specific cube. + >>> mesh_cube1 = mesh_cubes.extract_cube("sea_surface_temperature") + >>> print(mesh_cube1) + sea_surface_temperature / (K) (-- : 1; -- : 13824) + Mesh coordinates: + latitude - x + longitude - x + Auxiliary coordinates: + time x - + Cell methods: + mean time (300 s) + mean time_counter + Attributes: + Conventions UGRID + description Created by xios + interval_operation 300 s + interval_write 1 d + name lfric_surface + online_operation average + timeStamp 2020-Feb-07 16:23:14 GMT + title Created by xios + uuid 489bcef5-3d1c-4529-be42-4ab5f8c8497b + + # Load the target grid. + >>> sample_grid = load_cube(grid_file) + >>> print(sample_grid) + sample_grid / (unknown) (latitude: 180; longitude: 360) + Dimension coordinates: + latitude x - + longitude - x + Attributes: + Conventions CF-1.7 + + # Initialise the regridder. + >>> rg = MeshToGridESMFRegridder(mesh_cube1, sample_grid) + + # Regrid the mesh cube cube. + >>> result1 = rg(mesh_cube1) + >>> print(result1) + sea_surface_temperature / (K) (-- : 1; latitude: 180; longitude: 360) + Dimension coordinates: + latitude - x - + longitude - - x + Auxiliary coordinates: + time x - - + Cell methods: + mean time (300 s) + mean time_counter + Attributes: + Conventions UGRID + description Created by xios + interval_operation 300 s + interval_write 1 d + name lfric_surface + online_operation average + timeStamp 2020-Feb-07 16:23:14 GMT + title Created by xios + uuid 489bcef5-3d1c-4529-be42-4ab5f8c8497b + +.. note:: + + **All** :class:`~iris.cube.Cube` :attr:`~iris.cube.Cube.attributes` are + retained when regridding, so watch out for any attributes that reference + the format (there are several in these examples) - you may want to manually + remove them to avoid later confusion. + +The initialisation process is computationally expensive so we use caching to +improve performance. Once a regridder has been initialised, it can be used on +any :class:`~iris.cube.Cube` which has been defined on the same +:class:`~iris.experimental.ugrid.Mesh` (or on the same **grid** in the case of +:class:`~esmf_regrid.experimental.unstructured_scheme.GridToMeshESMFRegridder`). +Since calling a regridder is usually a lot faster than initialising, reusing +regridders can save a lot of time. We can demonstrate the reuse of the +previously initialised regridder: + +.. dropdown:: :opticon:`code` + + .. code-block:: python + + # Extract a different cube defined on te same Mesh. + >>> mesh_cube2 = mesh_cubes.extract_cube("precipitation_flux") + >>> print(mesh_cube2) + precipitation_flux / (kg m-2 s-1) (-- : 1; -- : 13824) + Mesh coordinates: + latitude - x + longitude - x + Auxiliary coordinates: + time x - + Cell methods: + mean time (300 s) + mean time_counter + Attributes: + Conventions UGRID + description Created by xios + interval_operation 300 s + interval_write 1 d + name lfric_surface + online_operation average + timeStamp 2020-Feb-07 16:23:14 GMT + title Created by xios + uuid 489bcef5-3d1c-4529-be42-4ab5f8c8497b + + # Regrid the new mesh cube using the same regridder. + >>> result2 = rg(mesh_cube2) + >>> print(result2) + precipitation_flux / (kg m-2 s-1) (-- : 1; latitude: 180; longitude: 360) + Dimension coordinates: + latitude - x - + longitude - - x + Auxiliary coordinates: + time x - - + Cell methods: + mean time (300 s) + mean time_counter + Attributes: + Conventions UGRID + description Created by xios + interval_operation 300 s + interval_write 1 d + name lfric_surface + online_operation average + timeStamp 2020-Feb-07 16:23:14 GMT + title Created by xios + uuid 489bcef5-3d1c-4529-be42-4ab5f8c8497b + +Support also exists for saving and loading previously initialised regridders - +:func:`esmf_regrid.experimental.io.save_regridder` and +:func:`~esmf_regrid.experimental.io.load_regridder` - so that they can be +re-used by future scripts. + +Equality +-------- +.. |tagline: equality| replace:: |unchanged| + +.. rubric:: |tagline: equality| + +:class:`~iris.experimental.ugrid.Mesh` comparison is supported, and comparing +two ':class:`~iris.experimental.ugrid.Mesh`-:class:`~iris.cube.Cube`\s' will +include a comparison of the respective +:class:`~iris.experimental.ugrid.Mesh`\es, with no extra action needed by the +user. + +.. note:: + + Keep an eye on memory demand when comparing large + :class:`~iris.experimental.ugrid.Mesh`\es, but note that + :class:`~iris.experimental.ugrid.Mesh`\ equality is enabled for lazy + processing (:doc:`/userguide/real_and_lazy_data`), so if the + :class:`~iris.experimental.ugrid.Mesh`\es being compared are lazy the + process will use less memory than their total size. + +Combining Cubes +--------------- +.. |tagline: combining cubes| replace:: |pending| + +.. rubric:: |tagline: combining cubes| + +Merging or concatenating :class:`~iris.cube.Cube`\s (described in +:doc:`/userguide/merge_and_concat`) with two different +:class:`~iris.experimental.ugrid.Mesh`\es is not possible - a +:class:`~iris.cube.Cube` must be associated with just a single +:class:`~iris.experimental.ugrid.Mesh`, and merge/concatenate are not yet +capable of combining multiple :class:`~iris.experimental.ugrid.Mesh`\es into +one. + +:class:`~iris.cube.Cube`\s that include +:class:`~iris.experimental.ugrid.MeshCoord`\s can still be merged/concatenated +on dimensions other than the :meth:`~iris.cube.Cube.mesh_dim`, since such +:class:`~iris.cube.Cube`\s will by definition share the same +:class:`~iris.experimental.ugrid.Mesh`. + +.. seealso:: + + You may wish to investigate + :func:`iris.experimental.ugrid.recombine_submeshes`, which can be used + for a very specific type of :class:`~iris.experimental.ugrid.Mesh` + combination not detailed here. + +Arithmetic +---------- +.. |tagline: arithmetic| replace:: |pending| + +.. rubric:: |tagline: arithmetic| + +:class:`~iris.cube.Cube` Arithmetic (described in :doc:`/userguide/cube_maths`) +has not yet been adapted to handle :class:`~iris.cube.Cube`\s that include +:class:`~iris.experimental.ugrid.MeshCoord`\s. + + +.. todo: + Enumerate other popular operations that aren't yet possible + (and are they planned soon?) + +.. |new| replace:: ✨ New +.. |unchanged| replace:: ♻️ Unchanged +.. |different| replace:: ⚠️ Different +.. |pending| replace:: 🚧 Support Pending \ No newline at end of file diff --git a/docs/src/further_topics/ugrid/partner_packages.rst b/docs/src/further_topics/ugrid/partner_packages.rst new file mode 100644 index 0000000000..8e36f4ffc2 --- /dev/null +++ b/docs/src/further_topics/ugrid/partner_packages.rst @@ -0,0 +1,100 @@ +.. _ugrid partners: + +Iris' Mesh Partner Packages +**************************** +Python is an easy to use language and has formed a very strong collaborative +scientific community, which is why Iris is written in Python. *Performant* +Python relies on calls down to low level languages like C, which is ideal for +structured grid work since +they can be directly represented as NumPy arrays. This is more difficult when +working with unstructured meshes where extra steps are needed to determine data +position (:ref:`see the data model detail `), and we need +to find ways of again passing the operations down to more optimised languages. + +The Iris team are therefore developing 'wrapper' packages, which make it quick +and easy to analyse Iris mesh data via some popular Python packages that use +powerful tools under the hood, working in C and other languages. + +These solutions have been placed in their own 'partner packages' for several +reasons: + +* Can be useful to others who are not using Iris. + + * Everyone working with multi-dimensional geographic datasets shares common + problems that need solving. + * Wider user base = stronger community = better solutions. + +* Only some Iris users will need them - they are **optional** Iris dependencies. + + * They introduce a lot of new API. + * They introduce new large dependencies that take time to install and need + disk space. + +Below you can learn more about the partner packages and how they are useful. +Specifics of what operations would require their installation can be found in: +:doc:`operations`. + +.. important:: **Experimental** + + As with Iris' mesh support, these packages are still in the + experimental stages. They would love your feedback, but as immature + packages their API, documentation, test coverage and CI are still + 'under construction'. + + +.. _`ugrid geovista`: + +`GeoVista`_ +=========== +.. image:: images/geovistalogo.svg + :width: 300 + :class: no-scaled-link + +.. rubric:: "Cartographic rendering and mesh analytics powered by `PyVista`_" + +PyVista is described as "VTK for humans" - VTK is a very powerful toolkit for +working with meshes, and PyVista brings that power into the Python ecosystem. +GeoVista in turn makes it easy to use PyVista specifically for cartographic +work, designed from the start with the Iris +:class:`~iris.experimental.ugrid.Mesh` in mind. + +Applications +------------ +* Interactively plot mesh data: + + * On a 3D globe. + * On your favourite projection. + +* Extract a specific region from a mesh. +* Combine multiple meshes into one. + +.. _`ugrid iris-esmf-regrid`: + +`iris-esmf-regrid`_ +=================== +.. image:: images/iris-esmf-regrid.svg + :width: 300 + :class: no-scaled-link + +.. rubric:: "A collection of structured and unstructured ESMF regridding schemes for Iris" + +ESMF provide a sophisticated, performant regridding utility that supports a +variety of regridding types with both structured grids and unstructured meshes, +and this also has a flexible Python interface - ESMPy. iris-esmf-regrid takes +advantage of having a specific use-case - regridding Iris +:class:`~iris.cube.Cube`\s - to provide ESMPy-Iris wrappers that make the +process as easy as possible, with highly optimised performance. + +Applications +------------ +* Regrid structured to unstructured. +* Regrid unstructured to structured. +* Regrid with dask integration, computing in parallel and maintaining data + laziness. +* | Save a prepared regridder for re-use in subsequent runs. + | Regridders can even be re-used on sources with different masks - a + significant efficiency gain. + +.. _GeoVista: https://github.com/bjlittle/geovista +.. _PyVista: https://docs.pyvista.org/index.html +.. _iris-esmf-regrid: https://github.com/SciTools-incubator/iris-esmf-regrid diff --git a/docs/src/index.rst b/docs/src/index.rst index e59a6f0527..d6fc5f2f7e 100644 --- a/docs/src/index.rst +++ b/docs/src/index.rst @@ -131,6 +131,7 @@ For **Iris 2.4** and earlier documentation please see the further_topics/metadata further_topics/lenient_metadata further_topics/lenient_maths + further_topics/ugrid/index .. toctree:: diff --git a/docs/src/whatsnew/latest.rst b/docs/src/whatsnew/latest.rst index 35628c4355..6ba122e4d7 100644 --- a/docs/src/whatsnew/latest.rst +++ b/docs/src/whatsnew/latest.rst @@ -16,8 +16,10 @@ This document explains the changes made to Iris for this release The highlights for this minor release of Iris include: - * We've added experimental support for `UGRID`_ meshes which can now be loaded - and attached to a cube. + * We've added experimental support for + :ref:`Meshes `, which can now be loaded and + attached to a cube. Mesh support is based on the based on `CF-UGRID`_ + model. * We've also dropped support for ``Python 3.7``. And finally, get in touch with us on `GitHub`_ if you have any issues or @@ -37,7 +39,7 @@ This document explains the changes made to Iris for this release =========== #. `@bjlittle`_, `@pp-mo`_, `@trexfeathers`_ and `@stephenworsley`_ added - support for unstructured meshes, as described by `UGRID`_. This involved + support for :ref:`unstructured meshes `. This involved adding a data model (:pull:`3968`, :pull:`4014`, :pull:`4027`, :pull:`4036`, :pull:`4053`, :pull:`4439`) and API (:pull:`4063`, :pull:`4064`), and supporting representation (:pull:`4033`, :pull:`4054`) of data on meshes. @@ -54,14 +56,14 @@ This document explains the changes made to Iris for this release :class:`~iris.cube.Cube` via a :class:`~iris.experimental.ugrid.mesh.MeshCoord`. #. `@trexfeathers`_ added support for loading unstructured mesh data from netcdf data, - for files using the `UGRID`_ conventions. + for files using the `CF-UGRID`_ conventions. The context manager :obj:`~iris.experimental.ugrid.load.PARSE_UGRID_ON_LOAD` provides a way to load UGRID files so that :class:`~iris.cube.Cube`\ s can be returned with a :class:`~iris.experimental.ugrid.mesh.Mesh` attached. (:pull:`4058`). -#. `@pp-mo`_ added support to save cubes with meshes to netcdf files, using the - `UGRID`_ conventions. +#. `@pp-mo`_ added support to save cubes with :ref:`meshes ` to netcdf + files, using the `CF-UGRID`_ conventions. The existing :meth:`iris.save` function now does this, when saving cubes with meshes. A routine :meth:`iris.experimental.ugrid.save.save_mesh` allows saving :class:`~iris.experimental.ugrid.mesh.Mesh` objects to netcdf *without* any associated data @@ -82,7 +84,7 @@ This document explains the changes made to Iris for this release :class:`~iris.coords.AuxCoord` :attr:`~iris.coords.AuxCoord.points` and :class:`~iris.experimental.ugrid.mesh.Connectivity` :attr:`~iris.experimental.ugrid.mesh.Connectivity.indices` under the - `UGRID`_ model. (:pull:`4375`) + :ref:`mesh model `. (:pull:`4375`) #. `@bsherratt`_ added a `threshold` parameter to :meth:`~iris.cube.Cube.intersection` (:pull:`4363`) @@ -309,7 +311,7 @@ This document explains the changes made to Iris for this release :func:`~iris.analysis.cartography.wrap_lons` and updated affected tests using assertArrayAllClose following :issue:`3993`. (:pull:`4421`) - + #. `@rcomer`_ updated some tests to work with Matplotlib v3.5. (:pull:`4428`) #. `@rcomer`_ applied minor fixes to some regridding tests. (:pull:`4432`) @@ -343,7 +345,6 @@ This document explains the changes made to Iris for this release .. _GitHub: https://github.com/SciTools/iris/issues/new/choose .. _NEP-29: https://numpy.org/neps/nep-0029-deprecation_policy.html -.. _UGRID: http://ugrid-conventions.github.io/ugrid-conventions/ .. _sort-all: https://github.com/aio-libs/sort-all .. _faster documentation building: https://docs.readthedocs.io/en/stable/guides/conda.html#making-builds-faster-with-mamba .. _Metarelate: http://www.metarelate.net/