Issue #1372 test regridder

docs/api/changelog.rst

@@ -9,6 +9,17 @@ The format is based on `Keep a Changelog`_, and this project adheres to
 [Unreleased]
 ------------
 
+From this release on, we recommend using `xugrid's regridding utilities
+<https://deltares.github.io/xugrid/examples/regridder_overview.html>`_ for
+regridding individual grids instead of :class:`imod.prepare.Regridder`. Xugrid's
+regridders are tested to be about 10 times faster than
+:class:`imod.prepare.Regridder`. There is one small difference: xugrid's
+``xugrid.BaryCentricInterpolator`` considers sample points of the destination
+grid that lie on the source grid's cell edges to be inside, whereas
+:class:`imod.prepare.Regridder` considers them to be outside. This difference is
+negligible for most applications, but might create slightly fewer ``np.nan``
+values than before.
+
 Removed
 ~~~~~~~
 - ``imod.flow`` module has been removed for generating iMODFLOW models. Use

docs/faq/how-do-i/modification.rst

@@ -99,19 +99,25 @@ Make sure the grids have the same spatial coordinates.
 Change cellsize (and extent)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+`xugrid's regridding functionality
+<https://deltares.github.io/xugrid/examples/regridder_overview.html>`_ allows
+regridding structured grids (next to regridding unstructured grids).
+
 Nearest neighbor:
 
 .. code-block:: python
 
-    regridder = imod.prepare.Regridder(source, destination, method="nearest")
-    out = regridder.regrid(source)
+    import xugrid as xu
+    regridder = xu.CentroidLocatorRegridder(source=source, target=like)
+    result = regridder.regrid(source)
 
 Area weighted mean:
 
 .. code-block:: python
 
-    regridder = imod.prepare.Regridder(source, destination, method="mean")
-    out = regridder.regrid(source)
+    import xugrid as xu
+    regridder = xu.OverlapRegridder(source=source, target=like, method="mean")
+    result = regridder.regrid(source)
 
 Change time resolution
 ~~~~~~~~~~~~~~~~~~~~~~

imod/prepare/regrid.py

@@ -299,12 +299,14 @@ def _nd_regrid(src, dst, src_coords, dst_coords, iter_regrid, use_relative_weigh
 
 
 WARNING_MSG = textwrap.dedent(
-    """{name} is deprecated and we plan to remove it in the
-    final v1.0 release. Use the regridder in xugrid instead. To regrid a single
-    array, see:
+    """{name} is deprecated so we will remove it eventually. Use the
+    the regridder in ``xugrid`` instead, which is around 10 times faster. To
+    regrid a single array, see:
     https://deltares.github.io/xugrid/examples/regridder_overview.html. To
-    regrid Modflow6 packages or entire simulations, see the user guide:
-    https://deltares.github.io/imod-python/user-guide/08-regridding.html."""
+    regrid Modflow6 packages or entire simulations, see the iMOD Python user
+    guide:
+    https://deltares.github.io/imod-python/user-guide/08-regridding.html.
+    """
 )
 
 
@@ -313,6 +315,22 @@ class Regridder(object):
     Object to repeatedly regrid similar objects. Compiles once on first call,
     can then be repeatedly called without JIT compilation overhead.
 
+    .. caution::
+
+        ``imod.prepare.Regridder`` is deprecated so we will remove it
+        eventually. Use the regridder in ``xugrid`` instead, which is around 10
+        times faster. It is as simple as:
+
+        >>> import xugrid as xu
+        >>> regridder = xu.OverlapRegridder(source=source, target=like, method="mean")
+        >>> result = regridder.regrid(source)
+
+        For more information, see:
+        https://deltares.github.io/xugrid/examples/regridder_overview.html. To
+        regrid MODFLOW6 packages or entire MODFLOW6 simulations, see the iMOD
+        Python user guide:
+        https://deltares.github.io/imod-python/user-guide/08-regridding.html.
+
     Attributes
     ----------
     method : str, function

imod/tests/test_regrid.py

@@ -5,6 +5,7 @@
 import pandas as pd
 import pytest
 import xarray as xr
+import xugrid as xu
 
 import imod
 
@@ -34,6 +35,17 @@ def weightedmean(values, weights):
     return vsum / wsum
 
 
+def weightedmean_xu(values, weights, workspace):
+    vsum = 0.0
+    wsum = 0.0
+    for i in range(values.size):
+        v = values[i]
+        w = weights[i]
+        vsum += w * v
+        wsum += w
+    return vsum / wsum
+
+
 def conductance(values, weights):
     v_agg = 0.0
     w_sum = 0.0
@@ -405,17 +417,21 @@ def test_regrid_mean2d(chunksize):
         ]
     )
     assert np.allclose(out.values, compare)
+    # Verify xugrid returns same result
+    out_xu = xu.OverlapRegridder(
+        source=source, target=like, method=weightedmean_xu
+    ).regrid(source)
+    assert np.allclose(out_xu.values, compare)
 
     # Now with chunking
     source = source.chunk({"x": chunksize, "y": chunksize})
     out = imod.prepare.Regridder(method=weightedmean).regrid(source, like)
     assert np.allclose(out.values, compare)
-
-    # Now with orthogonal chunks
-    source = source.chunk({"x": 2})
-    likecoords = {"y": dst_x}
-    like = xr.DataArray(np.empty(2), likecoords, ["y"])
-    out = imod.prepare.Regridder(method=weightedmean).regrid(source, like)
+    # Verify xugrid returns same result
+    out_xu = xu.OverlapRegridder(
+        source=source, target=like, method=weightedmean_xu
+    ).regrid(source)
+    assert np.allclose(out_xu.values, compare)
 
 
 @pytest.mark.parametrize("chunksize", [1, 2, 3])
@@ -446,13 +462,22 @@ def test_regrid_mean2d_over3darray(chunksize):
     )
     compare = np.empty((5, 2, 2))
     compare[:, ...] = compare_values
-
     assert np.allclose(out.values, compare)
+    # Verify xugrid returns same result
+    out_xu = xu.OverlapRegridder(
+        source=source, target=like, method=weightedmean_xu
+    ).regrid(source)
+    assert np.allclose(out_xu.values, compare)
 
     # Now with chunking
     source = source.chunk({"x": chunksize, "y": chunksize})
     out = imod.prepare.Regridder(method=weightedmean).regrid(source, like)
     assert np.allclose(out.values, compare)
+    # Verify xugrid returns same result
+    out_xu = xu.OverlapRegridder(
+        source=source, target=like, method=weightedmean_xu
+    ).regrid(source)
+    assert np.allclose(out_xu.values, compare)
 
 
 def test_regrid_conductance2d():
@@ -465,10 +490,16 @@ def test_regrid_conductance2d():
     src_da = xr.DataArray(
         [[10.0, 10.0], [10.0, 10.0]], {"y": [7.5, 2.5], "x": [2.5, 7.5]}, dims
     )
+    src_da[0, 0] = np.nan
 
     regridder = imod.prepare.Regridder(method=conductance, use_relative_weights=True)
     dst_da = regridder.regrid(src_da, like_da)
     assert float(src_da.sum()) == float(dst_da.sum())
+    # Verify xugrid returns same result
+    dst_da_xu = xu.RelativeOverlapRegridder(
+        source=src_da, target=like_da, method="conductance"
+    ).regrid(src_da)
+    assert float(src_da.sum()) == float(dst_da_xu.sum())
 
     # Second case, different domain, smaller cellsizes
     dx = np.array([2.5, 2.5, 2.5, 3.5])
@@ -478,19 +509,30 @@ def test_regrid_conductance2d():
     like_da = xr.DataArray(np.empty((4, 4)), coords, dims)
     dst_da = regridder.regrid(src_da, like_da)
     assert float(src_da.sum()) == float(dst_da.sum())
+    # Verify xugrid returns same result
+    dst_da_xu = xu.RelativeOverlapRegridder(
+        source=src_da, target=like_da, method="conductance"
+    ).regrid(src_da)
+    assert float(src_da.sum()) == float(dst_da_xu.sum())
 
     # Third case, same domain, small to large cellsizes
     y = np.arange(10.0, 0.0, -2.5) - 1.25
     x = np.arange(0.0, 10.0, 2.5) + 1.25
     coords = {"y": y, "x": x}
     dims = ("y", "x")
     src_da = xr.DataArray(np.full((4, 4), 10.0), coords, dims)
+    src_da[0, 0] = np.nan
     like_da = xr.DataArray(
         [[10.0, 10.0], [10.0, 10.0]], {"y": [7.5, 2.5], "x": [2.5, 7.5]}, dims
     )
 
     dst_da = regridder.regrid(src_da, like_da)
     assert float(src_da.sum()) == float(dst_da.sum())
+    # Verify xugrid returns same result
+    dst_da_xu = xu.RelativeOverlapRegridder(
+        source=src_da, target=like_da, method="conductance"
+    ).regrid(src_da)
+    assert float(src_da.sum()) == float(dst_da_xu.sum())
 
     # Fourth case, larger domain, small to large cellsizes
     like_da = xr.DataArray(
@@ -499,6 +541,11 @@ def test_regrid_conductance2d():
 
     dst_da = regridder.regrid(src_da, like_da)
     assert float(src_da.sum()) == float(dst_da.sum())
+    # Verify xugrid returns same result
+    dst_da_xu = xu.RelativeOverlapRegridder(
+        source=src_da, target=like_da, method="conductance"
+    ).regrid(src_da)
+    assert float(src_da.sum()) == float(dst_da_xu.sum())
 
 
 def test_regrid_errors():
@@ -541,7 +588,8 @@ def test_str_method():
     out = imod.prepare.Regridder(method="mean").regrid(source, like)
     assert np.allclose(out.values, compare)
 
-    out = imod.prepare.Regridder(method="nearest").regrid(source, like)
+    # Test if "nearest" name is still supported
+    imod.prepare.Regridder(method="nearest").regrid(source, like)
 
 
 def test_no_overlap():