to_facet method (Deltares#350)

* Create to_facet method, closes Deltares#173

* Add to_node, to_edge, to_face methods; API, changelog

* Make sure source dim isn't chunked before indexing

* Defalt to nmax instead of contributors for the new dim name

* Drop python 3.9

* Update lockfile

* Remove py 3.9 from CI

* Pin minimum dask version to 2025.1.0 and update pixi dependencies

---------

Co-authored-by: Huite <huite.bootsma@gmail.com>
Co-authored-by: JoerivanEngelen <joerivanengelen@hotmail.com>

Author: 3 people

.github/workflows/ci.yml

@@ -37,7 +37,6 @@ jobs:
           - py312
           - py311
           - py310 
-          - py309 
     steps:
       - name: Check out repo
         uses: actions/checkout@v4

docs/api.rst

@@ -123,6 +123,9 @@ UgridDataArray or UgridDataset.
     UgridDataArrayAccessor.rasterize
     UgridDataArrayAccessor.rasterize_like
     UgridDataArrayAccessor.reindex_like
+    UgridDataArrayAccessor.to_node
+    UgridDataArrayAccessor.to_edge
+    UgridDataArrayAccessor.to_face
     UgridDataArrayAccessor.to_periodic
     UgridDataArrayAccessor.to_nonperiodic
     UgridDataArrayAccessor.to_geodataframe

docs/changelog.rst

@@ -31,8 +31,12 @@ Added
   inserted vertices.
 - Added :class:`xugrid.NetworkGridder` to grid networks (Ugrid1d) to a 2D grid.
   Currently only support gridding to a Ugrid2d grid.
-- :meth:`xugrid.UgridDataArray.interpolate_na` will now also work for Ugrid1d
+- :meth:`xugrid.UgridDataArrayAccessor.interpolate_na` will now also work for Ugrid1d
   topologies.
+- :meth:`xugrid.UgridDataArrayAccessor.to_node`,
+  :meth:`xugrid.UgridDataArrayAccessor.to_edge`, and
+  :meth:`xugrid.UgridDataArrayAccessor.to_face` have been added to map data
+  from one facet of a grid to another.
 - Added ``tolerance`` argument to :meth:`xugrid.Ugrid1d.sel_points`,
   :meth:`xugrid.Ugrid1d.refine_by_vertices` and
   :meth:`xugrid.Ugrid2d.sel_points`, :meth:`xugrid.Ugrid2d.locate_centroids`.

pixi.lock

@@ -7,11 +7,12 @@ name = "xugrid"
 description = "Xarray extension for unstructured grids"
 readme = { file = "README.rst", content-type = "text/x-rst" }
 maintainers = [{ name = "Huite Bootsma", email = "huite.bootsma@deltares.nl" }]
-requires-python = ">=3.9"
+requires-python = ">=3.10"
 dependencies = [
     'pandas',
     'numba',
     'numba_celltree>=0.4.1',
+    'dask>=2025.1.0',
     'numpy',
     'pooch',
     'scipy',
@@ -25,7 +26,6 @@ classifiers = [
     'Programming Language :: Python',
     'Operating System :: OS Independent',
     'Programming Language :: Python :: 3',
-    'Programming Language :: Python :: 3.9',
     'Programming Language :: Python :: 3.10',
     'Programming Language :: Python :: 3.11',
     'Programming Language :: Python :: 3.12',
@@ -75,15 +75,14 @@ platforms = ["win-64", "linux-64", "osx-64", "osx-arm64"]
 
 [tool.pixi.pypi-dependencies]
 xugrid = { path = ".", editable = true }
-numba_celltree = ">=0.4.1"
 
 [tool.pixi.dependencies]
-dask = "*"
+dask = ">=2025.1.0"
 geopandas = "*"
 mapbox_earcut = "*"
 matplotlib-base = "*"
 netcdf4 = "*"
-# numba_celltree = ">=0.4.0"
+numba_celltree = ">=0.4.1"
 numpy = "*"
 pip = "*"
 pooch = "*"
@@ -94,7 +93,7 @@ pyproj = "*"
 pytest = "*"
 pytest-cases = "*"
 pytest-cov = "*"
-python = ">=3.9"
+python = ">=3.10"
 ruff = "*"
 shapely = ">=2.0"
 scipy = "*"
@@ -126,16 +125,12 @@ python = "3.11.*"
 [tool.pixi.feature.py310.dependencies]
 python = "3.10.*"
 
-[tool.pixi.feature.py309.dependencies]
-python = "3.9.*"
-
 [tool.pixi.environments]
 default = { features = ["py312"], solve-group = "py312" }
 py313 = { features = ["py313"], solve-group = "py313" }
 py312 = { features = ["py312"], solve-group = "py312" }
 py311 = ["py311"]
 py310 = ["py310"]
-py309 = ["py309"]
 
 [tool.ruff.lint]
 # See https://docs.astral.sh/ruff/rules/

pyproject.toml

@@ -87,6 +87,7 @@ def test_ugrid1d_properties():
     assert grid.edge_dimension == f"{NAME}_nEdges"
     assert grid.n_node == 3
     assert grid.n_edge == 2
+    assert grid.facets == {"node": grid.node_dimension, "edge": grid.edge_dimension}
     expected_coords = np.array([[0.0, 0.0], [1.0, 1.0], [2.0, 2.0]])
     assert np.allclose(grid.node_coordinates, expected_coords)
     assert np.allclose(grid.edge_x, [0.5, 1.5])

tests/test_ugrid1d.py

@@ -146,6 +146,11 @@ def test_ugrid2d_properties():
     assert grid.n_edge == 10
     assert grid.n_face == 4
     assert grid.n_max_node_per_face == 4
+    assert grid.facets == {
+        "node": grid.node_dimension,
+        "edge": grid.edge_dimension,
+        "face": grid.face_dimension,
+    }
     assert np.array_equal(grid.n_node_per_face, [4, 4, 3, 3])
     assert np.allclose(grid.node_coordinates, VERTICES)
     assert grid.bounds == (0.0, 0.0, 2.0, 2.0)

tests/test_ugrid2d.py

@@ -379,6 +379,38 @@ def test_broadcasted_laplace_interpolate(self):
         assert set(actual.dims) == set(nd_uda2.dims)
         assert isinstance(actual.data, dask.array.Array)
 
+    def test_to_facets(self):
+        face_da = self.uda
+        grid = face_da.ugrid.grid
+
+        with pytest.raises(ValueError, match="No conversion needed"):
+            face_da.ugrid.to_face()
+        with pytest.raises(ValueError, match="already exists"):
+            face_da.ugrid.to_face(grid.face_dimension)
+
+        node_da = face_da.ugrid.to_node()
+        edge_da = face_da.ugrid.to_edge()
+        assert isinstance(node_da, xugrid.UgridDataArray)
+        assert isinstance(edge_da, xugrid.UgridDataArray)
+
+        back1 = node_da.mean("nmax").ugrid.to_face()
+        back2 = edge_da.mean("nmax").ugrid.to_face()
+        cross1 = node_da.mean("nmax").ugrid.to_edge()
+        cross2 = edge_da.mean("nmax").ugrid.to_node()
+
+        assert isinstance(back1, xugrid.UgridDataArray)
+        assert isinstance(back2, xugrid.UgridDataArray)
+        assert isinstance(cross1, xugrid.UgridDataArray)
+        assert isinstance(cross2, xugrid.UgridDataArray)
+
+        assert node_da.dims == (grid.node_dimension, "nmax")
+        assert edge_da.dims == (grid.edge_dimension, "nmax")
+        assert back1.dims == (grid.face_dimension, "nmax")
+
+        # Check fill values, should contain two NaNs for the fill value.
+        assert back1[2:, -1].isnull().all()
+        assert back1.isnull().sum() == 2
+
     def test_to_dataset(self):
         uda2 = self.uda.copy()
         uda2.ugrid.obj.name = "test"
@@ -1524,6 +1556,24 @@ def test_laplace_interpolate_facets():
         assert np.allclose(actual, 1.0)
 
 
+def test_to_facets_1d():
+    uds = ugrid1d_ds()
+    with pytest.raises(ValueError, match="Cannot map to face"):
+        uds["a1d"].ugrid.to_face()
+    with pytest.raises(ValueError, match="No conversion needed"):
+        uds["a1d"].ugrid.to_node()
+    with pytest.raises(ValueError, match="No conversion needed"):
+        uds["b1d"].ugrid.to_edge()
+
+    grid = uds.ugrid.grid
+    to_edge = uds["a1d"].ugrid.to_edge()
+    to_node = uds["b1d"].ugrid.to_node()
+    assert isinstance(to_edge, xugrid.UgridDataArray)
+    assert isinstance(to_node, xugrid.UgridDataArray)
+    assert to_edge.dims == (grid.edge_dimension, "nmax")
+    assert to_node.dims == (grid.node_dimension, "nmax")
+
+
 def test_laplace_interpolate_1d():
     uda = ugrid1d_ds()["a1d"]
     uda[:] = 1.0

tests/test_ugrid_dataset.py

@@ -276,6 +276,123 @@ def to_nonperiodic(self, xmax: float):
         grid, obj = self.grid.to_nonperiodic(xmax=xmax, obj=self.obj)
         return UgridDataArray(obj, grid)
 
+    def _to_facet(self, facet: str, newdim: str):
+        """
+        Map the data from one facet to another.
+
+        Parameters
+        ----------
+        facet: str
+            node, edge, face
+        newdim: str
+            how to name the new dimension, e.g. three nodes
+            per triangle face, two nodes per edge, etc.
+        """
+        grid = self.grid
+        obj = self.obj
+
+        gridfacets = grid.facets
+        if facet not in gridfacets:
+            raise ValueError(
+                f"Cannot map to {facet} for a {type(grid).__name__} topology."
+            )
+
+        if newdim in obj.dims:
+            raise ValueError(
+                f"Dimension {newdim} already exists. Please provide a new dimension name."
+            )
+
+        source_dim = set(grid.dimensions).intersection(obj.dims).pop()
+        target_dim = getattr(grid, f"{facet}_dimension")
+        if source_dim == target_dim:
+            raise ValueError(
+                f"No conversion needed, data is already {facet}-associated."
+            )
+
+        # Find out on which facet we're currently located
+        source = {v: k for k, v in gridfacets.items()}[source_dim]
+        connectivity = grid.format_connectivity_as_dense(
+            getattr(grid, f"{facet}_{source}_connectivity")
+        )
+        indexer = xr.DataArray(connectivity, dims=(target_dim, newdim))
+        # Ensure the source dimension is not chunked for efficient indexing.
+        obj = obj.chunk({source_dim: -1})
+        # Set the fill values (-1) to NaN
+        mapped = obj.isel({source_dim: indexer}).where(connectivity != -1)
+        return UgridDataArray(mapped, grid)
+
+    def to_node(self, dim: str = "nmax"):
+        """
+        Map data to nodes.
+
+        Creates a new dimension representing the contributing source elements
+        for each node, as multiple faces/edges can connect to a single node.
+
+        Parameters
+        ----------
+        dim : str, optional
+
+        Returns
+        -------
+        mapped: UgridDataArray
+            A new UgridDataArray with data mapped to the nodes of the grid.
+
+        Examples
+        --------
+        Compute the mean elevation based on the surrounding faces for each node:
+
+        >>> node_elevation = face_elevation.to_node().mean("contributors")
+        """
+        return self._to_facet("node", dim)
+
+    def to_edge(self, dim: str = "nmax"):
+        """
+        Map data to edges.
+
+        Creates a new dimension representing the contributing source elements
+        for each node, as two nodes or two faces are connected to an edge.
+
+        Parameters
+        ----------
+        dim : str, optional
+
+        Returns
+        -------
+        mapped: UgridDataArray
+            A new UgridDataArray with data mapped to the edges of the grid.
+
+        Examples
+        --------
+        Compute the mean elevation based on the nodes for each edge:
+
+        >>> edge_elevation = node_elevation.to_edge().mean("contributors")
+        """
+        return self._to_facet("edge", dim)
+
+    def to_face(self, dim: str = "nmax"):
+        """
+        Map data to faces.
+
+        Creates a new dimension representing the contributing source elements
+        for each node, as two edges or multiple nodes are connected to a face.
+
+        Parameters
+        ----------
+        dim : str, optional
+
+        Returns
+        -------
+        mapped: UgridDataArray
+            A new UgridDataArray with data mapped to the faces of the grid.
+
+        Examples
+        --------
+        Compute the mean elevation based on the nodes for each face:
+
+        >>> face_elevation = node_elevation.to_face().mean("contributors")
+        """
+        return self._to_facet("face", dim)
+
     def intersect_line(
         self, start: Sequence[float], end: Sequence[float]
     ) -> xr.DataArray:

xugrid/core/dataarray_accessor.py

@@ -311,6 +311,13 @@ def get_coordinates(self, dim: str) -> FloatArray:
                 f"Expected {self.node_dimension} or {self.edge_dimension}; got: {dim}"
             )
 
+    @property
+    def facets(self) -> dict[str, str]:
+        return {
+            "node": self.node_dimension,
+            "edge": self.edge_dimension,
+        }
+
     def get_connectivity_matrix(self, dim: str, xy_weights: bool):
         """Return the connectivity matrix for the specified UGRID dimension."""
         if dim == self.node_dimension: