Iris to GeoVista conversion

lib/iris/experimental/geovista.py

@@ -0,0 +1,133 @@
+# 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):
+    """Helper function to get the coordinates from the cube."""
+    try:
+        coord = cube.coord(axis=axis, dim_coords=True)
+    except CoordinateNotFoundError:
+        coord = cube.coord(axis=axis)
+    return coord
+
+
+def cube_faces_to_polydata(cube, **kwargs):
+    """Function to convert a cube or the mesh attached to a cube into a polydata
+    object that can be used by GeoVista to generate plots.
+
+    Parameters
+    ----------
+    cube : :class:`~iris.cube.Cube`
+        Incoming cube containing the arrays or the mesh to be converted into the
+        polydata object.
+
+    **kwargs : dict
+        Additional keyword arguments to be passed to the Transform method.
+
+    """
+    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,
+        )
+    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 region_extraction(cube, polydata, region, **kwargs):
+    """Function to extract a region from a cube and its associated mesh and return
+    a new cube containing the region.
+
+    """
+    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("Must be on face or node.")
+
+        if cube.shape[mesh_dim] != polydata_length:
+            raise ValueError(
+                "The mesh on the cube and the polydata must have the" " same shape."
+            )
+
+        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/test_cube_to_poly.py

@@ -0,0 +1,47 @@
+import numpy as np
+
+from iris import load_cube
+from iris.experimental.geovista import cube_faces_to_polydata
+from iris.experimental.ugrid import PARSE_UGRID_ON_LOAD
+from iris.tests import get_data_path
+
+
+def test_integration_2d():
+    pass
+
+
+def test_integration_1d():
+    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")
+
+    polydata = cube_faces_to_polydata(global_cube)
+
+    assert polydata.GetNumberOfCells() == 13824
+    assert polydata.GetNumberOfPoints() == 13826
+    np.testing.assert_array_equal(polydata.active_scalars, global_cube.data)
+
+
+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_faces_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/unit/experimental/geobridge/__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."""

lib/iris/tests/unit/experimental/geobridge/test_cube_faces_to_polydata.py

@@ -0,0 +1,192 @@
+# 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 `iris.experimental.geobridge.cube_faces_to_polydata` function."""
+
+
+from unittest.mock import Mock
+
+from geovista import Transform
+import numpy as np
+import pytest
+
+import iris.analysis.cartography
+import iris.coord_systems
+from iris.experimental.geovista import cube_faces_to_polydata
+from iris.tests.stock import lat_lon_cube, sample_2d_latlons
+from iris.tests.stock.mesh import sample_mesh_cube
+
+
+@pytest.fixture()
+def cube_mesh():
+    return sample_mesh_cube()
+
+
+@pytest.fixture()
+def cube_1d():
+    sample_1d_cube = lat_lon_cube()
+    for coord in sample_1d_cube.dim_coords:
+        coord.coord_system = None
+    return sample_1d_cube
+
+
+@pytest.fixture()
+def cube_2d():
+    return sample_2d_latlons()
+
+
+@pytest.fixture()
+def default_cs():
+    return iris.coord_systems.GeogCS(
+        iris.analysis.cartography.DEFAULT_SPHERICAL_EARTH_RADIUS
+    )
+
+
+class ParentClass:
+    MOCKED_OPERATION = NotImplemented
+
+    @pytest.fixture()
+    def expected(self):
+        pass
+
+    @pytest.fixture()
+    def operation(self):
+        pass
+
+    @pytest.fixture()
+    def cube(self):
+        pass
+
+    @pytest.fixture()
+    def cube_with_crs(self, default_cs, cube):
+        cube_crs = cube.copy()
+        for coord in cube_crs.coords():
+            coord.coord_system = default_cs
+        return cube_crs
+
+    @pytest.fixture()
+    def mocked_operation(self):
+        mocking = Mock()
+        setattr(Transform, self.MOCKED_OPERATION, mocking)
+        return mocking
+
+    @staticmethod
+    def test_to_poly(expected, mocked_operation, cube):
+        cube_faces_to_polydata(cube)
+        actual = mocked_operation.call_args.kwargs
+        for key, expected_value in expected.items():
+            if hasattr(expected_value, "shape"):
+                np.testing.assert_array_equal(actual[key], expected_value)
+            else:
+                assert actual[key] == expected_value
+
+    @staticmethod
+    def test_to_poly_crs(mocked_operation, default_cs, cube_with_crs):
+        cube_faces_to_polydata(cube_with_crs)
+        actual = mocked_operation.call_args.kwargs
+        assert actual["crs"] == default_cs.as_cartopy_crs().proj4_init
+
+    @staticmethod
+    def test_to_poly_kwargs(mocked_operation, cube):
+        kwargs = {"test": "test"}
+        cube_faces_to_polydata(cube, **kwargs)
+        actual = mocked_operation.call_args.kwargs
+        assert actual["test"] == "test"
+
+
+class Test2dToPoly(ParentClass):
+    MOCKED_OPERATION = "from_2d"
+
+    @pytest.fixture()
+    def expected(self, cube_2d):
+        return {
+            "xs": cube_2d.coord(axis="X").contiguous_bounds(),
+            "ys": cube_2d.coord(axis="Y").contiguous_bounds(),
+            "data": cube_2d.data,
+            "name": cube_2d.name() + " / " + str(cube_2d.units),
+        }
+
+    @pytest.fixture()
+    def operation(self):
+        return "from_2d"
+
+    @pytest.fixture()
+    def cube(self, cube_2d):
+        return cube_2d
+
+
+class Test1dToPoly(ParentClass):
+    MOCKED_OPERATION = "from_1d"
+
+    @pytest.fixture()
+    def expected(self, cube_1d):
+        return {
+            "xs": cube_1d.coord(axis="X").contiguous_bounds(),
+            "ys": cube_1d.coord(axis="Y").contiguous_bounds(),
+            "data": cube_1d.data,
+            "name": cube_1d.name() + " / " + str(cube_1d.units),
+        }
+
+    @pytest.fixture()
+    def operation(self):
+        return "from_1d"
+
+    @pytest.fixture()
+    def cube(self, cube_1d):
+        return cube_1d
+
+
+class TestMeshToPoly(ParentClass):
+    MOCKED_OPERATION = "from_unstructured"
+
+    @pytest.fixture()
+    def expected(self, cube_mesh):
+        return {
+            "xs": cube_mesh.mesh.node_coords[0].points,
+            "ys": cube_mesh.mesh.node_coords[1].points,
+            "connectivity": cube_mesh.mesh.face_node_connectivity.indices_by_location(),
+            "data": cube_mesh.data,
+            "name": cube_mesh.name() + " / " + str(cube_mesh.units),
+            "start_index": 0,
+        }
+
+    @pytest.fixture()
+    def operation(self):
+        return "from_unstructured"
+
+    @pytest.fixture()
+    def cube(self, cube_mesh):
+        return cube_mesh
+
+    @pytest.mark.skip(reason="Meshes do not support crs currently")
+    def test_to_poly_crs(self, expected, actual):
+        return NotImplemented
+
+
+class TestExtras:
+    @pytest.fixture()
+    def cube_1d_2d(self, cube_2d):
+        my_cube = cube_2d.copy()
+        lat_coord = my_cube.aux_coords[0]
+        lat_coord_small = lat_coord[0]
+        lat_coord_small.bounds = None
+        lat_coord_small.points = np.arange(len(lat_coord_small.points))
+        my_cube.remove_coord(lat_coord)
+        my_cube.add_aux_coord(lat_coord_small, 1)
+        return my_cube
+
+    def test_not_1d_or_2d(self, cube_1d_2d):
+        with pytest.raises(
+            NotImplementedError,
+            match=r"Only 1D and 2D coordinates are supported",
+        ):
+            cube_faces_to_polydata(cube_1d_2d)
+
+    def test_no_mesh_or_2d(self, cube_1d):
+        cube = cube_1d[0]
+        with pytest.raises(
+            NotImplementedError,
+            match=r"Cube must have a mesh or have 2 dimensions",
+        ):
+            cube_faces_to_polydata(cube)