Iris to GeoVista conversion

docs/src/further_topics/ugrid/operations.rst

@@ -488,7 +488,7 @@ earlier:
 
 
         # We'll re-use this to plot some real global data later.
-        >>> def cube_faces_to_polydata(cube):
+        >>> def cube_to_polydata(cube):
         ...     lons, lats = cube.mesh.node_coords
         ...     face_node = cube.mesh.face_node_connectivity
         ...     indices = face_node.indices_by_location()
@@ -515,7 +515,7 @@ earlier:
 
         # Convert our mesh+data to a PolyData object.
         # Just plotting a single height level.
-        >>> face_polydata = cube_faces_to_polydata(face_cube[:, 0])
+        >>> face_polydata = cube_to_polydata(face_cube[:, 0])
         >>> print(face_polydata)
         PolyData (0x7ff4861ff4c0)
           N Cells:	2
@@ -590,7 +590,7 @@ Here's another example using a global cubed-sphere data set:
                 title                       Created by xios
                 uuid                        489bcef5-3d1c-4529-be42-4ab5f8c8497b
 
-        >>> global_polydata = cube_faces_to_polydata(global_cube)
+        >>> global_polydata = cube_to_polydata(global_cube)
         >>> print(global_polydata)
         PolyData (0x7f761b536160)
           N Cells:	13824

lib/iris/experimental/geovista.py

@@ -0,0 +1,121 @@
+# 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.
+"""Experimental module for using some GeoVista operations with Iris cubes."""
+
+from geovista import Transform
+from geovista.common import VTK_CELL_IDS, VTK_POINT_IDS
+
+from iris.exceptions import CoordinateNotFoundError
+from iris.experimental.ugrid import Mesh
+
+
+def _get_coord(cube, axis):
+    try:
+        coord = cube.coord(axis=axis, dim_coords=True)
+    except CoordinateNotFoundError:
+        coord = cube.coord(axis=axis)
+    return coord
+
+
+def cube_to_polydata(cube, **kwargs):
+    if cube.mesh:
+        if cube.ndim != 1:
+            raise NotImplementedError("Cubes with a mesh must be one dimensional")
+        lons, lats = cube.mesh.node_coords
+        face_node = cube.mesh.face_node_connectivity
+        indices = face_node.indices_by_location()
+
+        polydata = Transform.from_unstructured(
+            xs=lons.points,
+            ys=lats.points,
+            connectivity=indices,
+            data=cube.data,
+            name=f"{cube.name()} / {cube.units}",
+            start_index=face_node.start_index,
+            **kwargs,
+        )
+    # TODO: Add support for point clouds
+    elif cube.ndim == 2:
+        x_coord = _get_coord(cube, "X")
+        y_coord = _get_coord(cube, "Y")
+        transform_kwargs = dict(
+            xs=x_coord.contiguous_bounds(),
+            ys=y_coord.contiguous_bounds(),
+            data=cube.data,
+            name=f"{cube.name()} / {cube.units}",
+            **kwargs,
+        )
+        coord_system = cube.coord_system()
+        if coord_system:
+            transform_kwargs["crs"] = coord_system.as_cartopy_crs().proj4_init
+
+        if x_coord.ndim == 2 and y_coord.ndim == 2:
+            polydata = Transform.from_2d(**transform_kwargs)
+
+        elif x_coord.ndim == 1 and y_coord.ndim == 1:
+            polydata = Transform.from_1d(**transform_kwargs)
+
+        else:
+            raise NotImplementedError("Only 1D and 2D coordinates are supported")
+    else:
+        raise NotImplementedError("Cube must have a mesh or have 2 dimensions")
+
+    return polydata
+
+
+def extract_unstructured_region(cube, polydata, region, **kwargs):
+    if cube.mesh:
+        # Find what dimension the mesh is in on the cube
+        mesh_dim = cube.mesh_dim()
+        recreate_mesh = False
+
+        if cube.location == "face":
+            polydata_length = polydata.GetNumberOfCells()
+            indices_key = VTK_CELL_IDS
+            recreate_mesh = True
+        elif cube.location == "node":
+            polydata_length = polydata.GetNumberOfPoints()
+            indices_key = VTK_POINT_IDS
+        else:
+            raise NotImplementedError(
+                f"Must be on face or node. Found: {cube.location}."
+            )
+
+        if cube.shape[mesh_dim] != polydata_length:
+            raise ValueError(
+                f"The mesh on the cube and the polydata"
+                f"must have the same shape."
+                f" Found Mesh: {cube.shape[mesh_dim]},"
+                f" Polydata: {polydata_length}."
+            )
+
+        region_polydata = region.enclosed(polydata, **kwargs)
+        indices = region_polydata[indices_key]
+        if len(indices) == 0:
+            raise IndexError("No part of `polydata` falls within `region`.")
+
+        my_tuple = tuple(
+            [slice(None) if i != mesh_dim else indices for i in range(cube.ndim)]
+        )
+
+        region_cube = cube[my_tuple]
+
+        if recreate_mesh:
+            coords_on_mesh_dim = region_cube.coords(dimensions=mesh_dim)
+            new_mesh = Mesh.from_coords(
+                *[c for c in coords_on_mesh_dim if c.has_bounds()]
+            )
+
+            new_mesh_coords = new_mesh.to_MeshCoords(cube.location)
+
+            for coord in new_mesh_coords:
+                region_cube.remove_coord(coord.name())
+                region_cube.add_aux_coord(coord, mesh_dim)
+
+    # TODO: Support unstructured point based data without a mesh
+    else:
+        raise ValueError("Cube must have a mesh")
+
+    return region_cube

lib/iris/experimental/ugrid/mesh.py

@@ -862,6 +862,7 @@ def check_shape(array_name):
 
         #####
         # TODO: remove axis assignment once Mesh supports arbitrary coords.
+        # TODO: consider filtering coords as the first action in this method.
         axes_present = [guess_coord_axis(coord) for coord in coords]
         axes_required = ("X", "Y")
         if all([req in axes_present for req in axes_required]):

lib/iris/tests/integration/experimental/geovista/__init__.py

@@ -0,0 +1,5 @@
+# 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.
+"""Integration tests for the :mod:`iris.experimental.geovista` module."""

lib/iris/tests/integration/experimental/geovista/test_cube_to_poly.py

@@ -0,0 +1,66 @@
+# 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.
+"""Integration tests for the `iris.experimental.geovista.cube_to_polydata` function."""
+
+import numpy as np
+
+from iris import load_cube
+from iris.experimental.geovista import cube_to_polydata
+from iris.experimental.ugrid import PARSE_UGRID_ON_LOAD
+from iris.tests import get_data_path
+
+
+def test_integration_2d():
+    file_path = get_data_path(
+        [
+            "NetCDF",
+            "ORCA2",
+            "votemper.nc",
+        ]
+    )
+    with PARSE_UGRID_ON_LOAD.context():
+        global_cube = load_cube(file_path, "votemper")
+
+    polydata = cube_to_polydata(global_cube[0, 1, :])
+
+    assert polydata.GetNumberOfCells() == 26640
+    assert polydata.GetNumberOfPoints() == 26969
+
+
+def test_integration_1d():
+    file_path = get_data_path(
+        [
+            "NetCDF",
+            "global",
+            "xyt",
+            "SMALL_hires_wind_u_for_ipcc4.nc",
+        ]
+    )
+    with PARSE_UGRID_ON_LOAD.context():
+        global_cube = load_cube(file_path)
+
+    polydata = cube_to_polydata(global_cube[0, :])
+
+    assert polydata.GetNumberOfCells() == 51200
+    assert polydata.GetNumberOfPoints() == 51681
+
+
+def test_integration_mesh():
+    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")
+
+    polydata = cube_to_polydata(global_cube[0, :])
+
+    assert polydata.GetNumberOfCells() == 864
+    assert polydata.GetNumberOfPoints() == 866
+    np.testing.assert_array_equal(polydata.active_scalars, global_cube[0, :].data)

lib/iris/tests/integration/experimental/geovista/test_extract_unstructured_region.py

@@ -0,0 +1,34 @@
+# 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.
+"""Integration tests for the `iris.experimental.geovista.extract_unstructured_region` function."""
+
+from geovista.geodesic import BBox
+
+from iris import load_cube
+from iris.experimental.geovista import cube_to_polydata, extract_unstructured_region
+from iris.experimental.ugrid import PARSE_UGRID_ON_LOAD
+from iris.tests import get_data_path
+
+
+def test_face_region_extraction():
+    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")
+    polydata = cube_to_polydata(global_cube[0, :])
+    region = BBox(lons=[0, 70, 70, 0], lats=[-25, -25, 45, 45])
+
+    extracted_cube = extract_unstructured_region(
+        global_cube, polydata, region, preference="center"
+    )
+
+    assert extracted_cube.ndim == 2
+    assert extracted_cube.shape == (72, 101)

lib/iris/tests/unit/experimental/geovista/__init__.py

@@ -0,0 +1,5 @@
+# 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.
+"""Unit tests for the :mod:`iris.experimental.geovista` module."""