Add Convention.triangulate() function

Author: mx-moth

src/emsarray/conventions/_base.py

@@ -19,7 +19,7 @@
 from emsarray import plot as _plot
 from emsarray import utils
 from emsarray.exceptions import InvalidGeometryWarning, NoSuchCoordinateError
-from emsarray.operations import depth, point_extraction
+from emsarray.operations import depth, point_extraction, triangulate
 from emsarray.operations.cache import hash_attributes, hash_int, hash_string
 from emsarray.state import State
 from emsarray.types import Bounds, DataArrayOrName, Pathish
@@ -1949,6 +1949,51 @@ def hash_geometry(self, hash: "hashlib._Hash") -> None:
             # Hash dataset attributes
             hash_attributes(hash, data_array.attrs)
 
+    def make_triangulation(self) -> triangulate.Triangulation[GridKind]:
+        """
+        Triangulates the polygons in the dataset.
+        Subclasses may have improved implementations.
+
+        This requires the dataset to have a grid with polygonal geometry.
+        If there is an additional grid corresponding to the vertices of the polygons
+        then subclasses should use this information to inform the triangulation.
+
+        Returns
+        -------
+        tuple of vertices, triangles, and `cell_indexes`
+            A tuple of three numpy arrays is returned,
+            containing vertices, triangles, and cell indexes respectively.
+
+            `vertices` is a numpy array of shape (V, 2)
+            where V is the number of unique vertices in the dataset.
+            The vertex coordinates are in (x, y) or (lon, lat) order.
+            If the dataset has a vertex grid associated with the polygons
+            then this array replicates that data.
+
+            `triangles` is a numpy array of shape (T, 3)
+            where T is the number of triangles in the dataset.
+            Each triangle is a set of three vertex indexes.
+
+            `cell_indexes` is a numpy list of length T.
+            Each entry indicates which polygon from the dataset a triangle is a part of.
+
+        See also
+        --------
+        :func:`~emsarray.operations.triangulate.triangulate`
+        """
+        grid = self.default_grid
+        if not issubclass(grid.geometry_type, shapely.Polygon):
+            raise ValueError("Can not triangulate a dataset that does not have polygonal geometry")
+        polygons = grid.geometry
+        vertices = triangulate.find_unique_vertices(polygons)
+        polygon_vertex_indexes = triangulate.polygons_to_vertex_indexes(polygons, vertices)
+        vertex_coordinates, triangles, face_indexes = triangulate.triangulate(vertices, polygons, polygon_vertex_indexes)
+        return triangulate.Triangulation(
+            vertices=vertex_coordinates,
+            triangles=triangles,
+            face_indexes=face_indexes,
+            face_grid_kind=grid.grid_kind)
+
 
 type DimensionIndex[GridKind] = tuple[GridKind, *tuple[int, ...]]
 

src/emsarray/conventions/arakawa_c.py

@@ -19,6 +19,7 @@
 from xarray.core.dataset import DatasetCoordinates
 
 from emsarray import masking, utils
+from emsarray.operations import triangulate
 from emsarray.types import Pathish
 
 from ._base import DimensionConvention, Specificity
@@ -405,6 +406,24 @@ def make_clip_mask(
     def apply_clip_mask(self, clip_mask: xarray.Dataset, work_dir: Pathish) -> xarray.Dataset:
         return masking.mask_grid_dataset(self.dataset, clip_mask, work_dir)
 
+    def make_triangulation(self) -> triangulate.Triangulation:
+        vertices = self.grids[ArakawaCGridKind.node].geometry
+        polygons = self.grids[ArakawaCGridKind.face].geometry
+
+        j_size, i_size = self.node.shape
+        xx, yy = numpy.meshgrid(numpy.arange(i_size - 1), numpy.arange(j_size - 1) * i_size)
+        bottom_left_indexes = (xx + yy).flatten()
+        offsets = numpy.array([0, 1, i_size + 1, i_size])
+        polygon_vertex_indexes = bottom_left_indexes[:, None].repeat(4, axis=1) + offsets
+
+        vertex_coordinates, triangles, face_indexes = triangulate.triangulate(vertices, polygons, polygon_vertex_indexes)
+        return triangulate.Triangulation[ArakawaCGridKind](
+            vertices=vertex_coordinates,
+            triangles=triangles,
+            face_indexes=face_indexes,
+            face_grid_kind=ArakawaCGridKind.face,
+            vertex_grid_kind=ArakawaCGridKind.node)
+
 
 def c_mask_from_centres(
     face_mask: numpy.ndarray,

src/emsarray/conventions/ugrid.py

@@ -26,6 +26,7 @@
 from emsarray.exceptions import (
     ConventionViolationError, ConventionViolationWarning
 )
+from emsarray.operations import triangulate
 from emsarray.types import Bounds, Pathish
 
 from ._base import DimensionConvention, Specificity
@@ -1407,3 +1408,15 @@ def drop_geometry(self) -> xarray.Dataset:
         dataset = super().drop_geometry()
         dataset.attrs.pop('Conventions', None)
         return dataset
+
+    def make_triangulation(self) -> triangulate.Triangulation:
+        vertices = self.grids[UGridKind.node].geometry
+        polygons = self.grids[UGridKind.face].geometry
+        polygon_vertex_indexes = self.topology.face_node_array
+        vertex_coordinates, triangles, face_indexes = triangulate.triangulate(vertices, polygons, polygon_vertex_indexes)
+        return triangulate.Triangulation[UGridKind](
+            vertices=vertex_coordinates,
+            triangles=triangles,
+            face_indexes=face_indexes,
+            face_grid_kind=UGridKind.face,
+            vertex_grid_kind=UGridKind.node)