Add logic for earcutting polygons, and weights dataframes

docs/api.rst

@@ -32,6 +32,7 @@ Top-level functions
     merge
     merge_partitions
     burn_vector_geometry
+    earcut_triangulate_polygons
     polygonize
 
 UgridDataArray
@@ -335,6 +336,7 @@ UGRID2D Topology
     Ugrid2d.from_structured_intervals2d
     Ugrid2d.from_shapely
     Ugrid2d.to_shapely
+    Ugrid2d.earcut_triangulate_polygons
 
     Ugrid2d.plot
 

docs/changelog.rst

@@ -41,7 +41,13 @@ Added
   maximum connectivity dimensions and their corresponding size.
 - :attr:`xugrid.AbstractUgrid.max_connectivity_dimensions` which returns all
   maximum connectivity dimensions.
-
+- :func:`xugrid.earcut_triangulate_polygons` and
+  :meth:`xugrid.Ugrid2d.earcut_triangulate_polygons` have been added to break
+  down polygon geodataframes into a triangular mesh for further processing.
+- :meth:`xugrid.OverlapRegridder.weights_as_dataframe` has been added to
+  extract regridding weights (overlaps) from the regridders. This method is
+  also available for :class:`BarycentricInterpolator`,
+  :class:`CentroidLocatorRegridder`, and :class:`RelativeOverlapRegridder`.
 
 Changed
 ~~~~~~~

tests/test_regrid/test_regridder.py

@@ -1,4 +1,5 @@
 import numpy as np
+import pandas as pd
 import pytest
 import xarray as xr
 
@@ -46,6 +47,32 @@ def test_structured_to_unstructured(
     assert isinstance(actual, xu.UgridDataArray)
 
 
+@pytest.mark.parametrize(
+    "regridder_class",
+    [
+        CentroidLocatorRegridder,
+        OverlapRegridder,
+        RelativeOverlapRegridder,
+        BarycentricInterpolator,
+    ],
+)
+def test_weights_as_dataframe(
+    regridder_class,
+    disk,
+    quads_structured,
+):
+    regridder = regridder_class(quads_structured, disk)
+    df = regridder.weights_as_dataframe()
+    assert isinstance(df, pd.DataFrame)
+    assert "source_index" in df
+    assert "target_index" in df
+    assert "weight" in df
+
+    regridder._weights = None
+    with pytest.raises(ValueError):
+        regridder.weights_as_dataframe()
+
+
 def test_centroid_locator_regridder_structured(
     grid_data_a, grid_data_a_layered, grid_data_b, expected_results_centroid
 ):

tests/test_ugrid2d.py

@@ -1387,3 +1387,58 @@ def test_equals():
     assert not grid.equals(xr_grid)
     grid_copy.attrs["attr"] = "something_else"
     assert not grid.equals(grid_copy)
+
+
+def test_earcut_triangulate_polygons():
+    with pytest.raises(TypeError):
+        xugrid.Ugrid2d.earcut_triangulate_polygons("abc")
+
+    x = np.array([0.0, 1.0, 2.0])
+    y = np.array([0.0, 0.0, 0.0])
+    gdf = gpd.GeoDataFrame(geometry=[shapely.linestrings(x, y)])
+    with pytest.raises(TypeError):
+        xugrid.Ugrid2d.earcut_triangulate_polygons(gdf)
+
+    xy = np.array(
+        [
+            [0.0, 0.0],
+            [1.0, 0.0],
+            [1.0, 1.0],
+            [0.0, 1.0],
+            [0.0, 0.0],
+        ]
+    )
+    hole = np.array(
+        [
+            [
+                [0.25, 0.25],
+                [0.75, 0.25],
+                [0.75, 0.75],
+                [0.25, 0.25],
+            ]
+        ]
+    )
+    polygon = shapely.polygons(xy, holes=hole)
+    gdf = gpd.GeoDataFrame(geometry=[polygon])
+
+    grid = xugrid.Ugrid2d.earcut_triangulate_polygons(polygons=gdf)
+    assert isinstance(grid, xugrid.Ugrid2d)
+    assert grid.n_face == 7
+
+    grid, index = xugrid.Ugrid2d.earcut_triangulate_polygons(
+        polygons=gdf, return_index=True
+    )
+    assert isinstance(grid, xugrid.Ugrid2d)
+    assert isinstance(index, np.ndarray)
+    assert np.array_equal(index, np.zeros(7, dtype=int))
+
+    gdf = gpd.GeoDataFrame(data={"a": [10.0], "b": [20.0]}, geometry=[polygon])
+    uda = xugrid.earcut_triangulate_polygons(polygons=gdf)
+    assert isinstance(uda, xugrid.UgridDataArray)
+    assert np.allclose(uda.to_numpy(), 0)
+
+    uda = xugrid.earcut_triangulate_polygons(polygons=gdf, column="a")
+    assert np.allclose(uda.to_numpy(), 10.0)
+
+    uda = xugrid.earcut_triangulate_polygons(polygons=gdf, column="b")
+    assert np.allclose(uda.to_numpy(), 20.0)

xugrid/__init__.py

@@ -20,7 +20,7 @@
     OverlapRegridder,
     RelativeOverlapRegridder,
 )
-from xugrid.ugrid.burn import burn_vector_geometry
+from xugrid.ugrid.burn import burn_vector_geometry, earcut_triangulate_polygons
 from xugrid.ugrid.conventions import UgridRolesAccessor
 from xugrid.ugrid.partitioning import merge_partitions
 from xugrid.ugrid.polygonize import polygonize
@@ -51,6 +51,7 @@
     "OverlapRegridder",
     "RelativeOverlapRegridder",
     "burn_vector_geometry",
+    "earcut_triangulate_polygons",
     "UgridRolesAccessor",
     "merge_partitions",
     "polygonize",

xugrid/regrid/regridder.py

@@ -4,6 +4,7 @@
 
 import numba
 import numpy as np
+import pandas as pd
 import xarray as xr
 
 import xugrid as xu
@@ -85,10 +86,12 @@ def __init__(
     ):
         self._source = setup_grid(source)
         self._target = setup_grid(target)
+        self._weights = None
         self._compute_weights(self._source, self._target)
         return
 
-    @abc.abstractproperty
+    @property
+    @abc.abstractmethod
     def weights(self):
         pass
 
@@ -241,6 +244,31 @@ def to_dataset(self) -> xr.Dataset:
         target_ds = self._target.to_dataset("__target")
         return xr.merge((weights_ds, source_ds, target_ds))
 
+    def weights_as_dataframe(self) -> pd.DataFrame:
+        """
+        Return the weights as a three column dataframe:
+
+        * source index
+        * target index
+        * weight
+
+        Returns
+        -------
+        weights: pd.DataFrame
+        """
+        matrix = self._weights
+        if matrix is None:
+            raise ValueError("Weights have not been computed yet.")
+        if isinstance(matrix, MatrixCSR):
+            matrix = matrix.to_coo()
+        return pd.DataFrame(
+            {
+                "target_index": matrix.row,
+                "source_index": matrix.col,
+                "weight": matrix.data,
+            }
+        )
+
     @staticmethod
     def _csr_from_dataset(dataset: xr.Dataset) -> MatrixCSR:
         """

xugrid/ugrid/burn.py

@@ -22,6 +22,12 @@
 except ImportError:
     shapely = MissingOptionalModule("shapely")
 
+try:
+    import mapbox_earcut
+
+except ImportError:
+    mapbox_earcut = MissingOptionalModule("mapbox_earcut")
+
 
 @nb.njit(inline="always")
 def in_bounds(p: Point, a: Point, b: Point) -> bool:
@@ -227,7 +233,7 @@ def burn_vector_geometry(
     column: str | None = None,
     fill: Union[int, float] = np.nan,
     all_touched: bool = False,
-) -> None:
+) -> xugrid.UgridDataArray:
     """
     Burn vector geometries (points, lines, polygons) into a Ugrid2d mesh.
 
@@ -262,7 +268,7 @@ def burn_vector_geometry(
     GEOM_NAMES = {v: k for k, v in shapely.GeometryType.__members__.items()}
 
     if not isinstance(gdf, gpd.GeoDataFrame):
-        raise TypeError(f"gdf must be GeoDataFrame, received: {type(like).__name__}")
+        raise TypeError(f"gdf must be GeoDataFrame, received: {type(gdf).__name__}")
     if isinstance(like, (xugrid.UgridDataArray, xugrid.UgridDataset)):
         like = like.ugrid.grid
     if not isinstance(like, xugrid.Ugrid2d):
@@ -307,3 +313,103 @@ def burn_vector_geometry(
         obj=xr.DataArray(output, dims=[like.face_dimension], name=column),
         grid=like,
     )
+
+
+def grid_from_earcut_polygons(
+    polygons: "geopandas.GeoDataFrame",  # type: ignore # noqa
+    return_index: bool = False,
+):
+    import geopandas as gpd
+
+    if not isinstance(polygons, gpd.GeoDataFrame):
+        raise TypeError(f"Expected GeoDataFrame, received: {type(polygons).__name__}")
+
+    geometry = polygons.geometry
+    POLYGON = shapely.GeometryType.POLYGON
+    GEOM_NAMES = {v: k for k, v in shapely.GeometryType.__members__.items()}
+    geometry_id = shapely.get_type_id(geometry)
+
+    if not (geometry_id == POLYGON).all():
+        GEOM_NAMES = {v: k for k, v in shapely.GeometryType.__members__.items()}
+        received = ", ".join(
+            [GEOM_NAMES[geom_id] for geom_id in np.unique(geometry_id)]
+        )
+        raise TypeError(
+            "geometry contains unsupported geometry types. Can only triangulate "
+            f"Polygon geometries. Received: {received}"
+        )
+
+    # Shapely does not provide a vectorized manner to get the interior rings
+    # easily, an index argument is always required, which is a poor fit for
+    # heterogeneous polygons (where some have no holes, and some may have
+    # hundreds).
+    # map_box_earcut is also not vectorized for polygons, so we simply loop
+    # over the geometries here.
+    exterior_coordinates = [
+        shapely.get_coordinates(exterior) for exterior in geometry.exterior
+    ]
+    interior_coordinates = [
+        [shapely.get_coordinates(p_interior) for p_interior in p_interiors]
+        for p_interiors in geometry.interiors
+    ]
+
+    all_triangles = []
+    offset = 0
+    for exterior, interiors in zip(exterior_coordinates, interior_coordinates):
+        rings = np.cumsum([len(exterior)] + [len(interior) for interior in interiors])
+        vertices = np.vstack([exterior] + interiors).astype(np.float64)
+        triangles = mapbox_earcut.triangulate_float64(vertices, rings).reshape((-1, 3))
+        all_triangles.append(triangles + offset)
+        offset += len(vertices)
+
+    face_nodes = np.concatenate(all_triangles).reshape((-1, 3))
+    all_vertices = shapely.get_coordinates(geometry)
+    node_x = all_vertices[:, 0]
+    node_y = all_vertices[:, 1]
+    grid = xugrid.Ugrid2d(node_x, node_y, -1, face_nodes)
+
+    if return_index:
+        n_triangles = [len(triangles) for triangles in all_triangles]
+        index = np.repeat(np.arange(len(geometry)), n_triangles)
+        return grid, index
+    else:
+        return grid
+
+
+def earcut_triangulate_polygons(
+    polygons: "geopandas.GeoDataframe",  # type: ignore # noqa
+    column: str | None = None,
+) -> xugrid.UgridDataArray:
+    """
+    Break down polygons using mapbox_earcut, and create a mesh from the
+    resulting triangles.
+
+    If no ``column`` argument is provided, the polygon index will be assigned
+    to the grid faces.
+
+    Parameters
+    ----------
+    polygons: geopandas.GeoDataFrame
+        Polygons to convert to triangles.
+    column: str, optional
+        Name of the geodataframe column of which to the values to assign
+        to the grid faces.
+
+    Returns
+    -------
+    triangulated: UgridDataArray
+    """
+    grid, index = grid_from_earcut_polygons(polygons, return_index=True)
+
+    if column is not None:
+        da = (
+            polygons[column]
+            .reset_index(drop=True)
+            .to_xarray()
+            .isel(index=index)
+            .rename({"index": grid.face_dimension})
+        )
+    else:
+        da = xr.DataArray(data=index, dims=(grid.face_dimension,))
+
+    return xugrid.UgridDataArray(da, grid)

xugrid/ugrid/ugrid2d.py

@@ -1936,6 +1936,28 @@ def from_polygon(
         ugrid._meshkernel = _mesh_kernel
         return ugrid
 
+    @staticmethod
+    def earcut_triangulate_polygons(polygons, return_index: bool = False):
+        """
+        Break down polygons using mapbox_earcut, and create a mesh from the
+        resulting triangles.
+
+        Parameters
+        ----------
+        polygons: ndarray of shapely polygons
+        return_index: bool, default is False.
+
+        Returns
+        -------
+        grid: xugrid.Ugrid2d
+        index: ndarray of integer, optional
+            The polygon index for each triangle. Only provided if ``return_index``
+            is True.
+        """
+        return xugrid.ugrid.burn.grid_from_earcut_polygons(
+            polygons, return_index=return_index
+        )
+
     @classmethod
     def from_geodataframe(cls, geodataframe: "geopandas.GeoDataFrame") -> "Ugrid2d":  # type: ignore # noqa
         """