pygmt.grdfill: Deprecate parameter 'mode' and add parameters 'constantfill/gridfill/neighborfill/splinefill' (deprecated since v0.15.0 and will be removed in v0.19.0) (#3855)

Co-authored-by: Michael Grund <[email protected]>
Co-authored-by: Yvonne Fröhlich <[email protected]>

Author: 3 people

pygmt/src/grdfill.py

@@ -2,6 +2,8 @@
 grdfill - Interpolate across holes in a grid.
 """
 
+import warnings
+
 import xarray as xr
 from pygmt.clib import Session
 from pygmt.exceptions import GMTInvalidInput
@@ -16,12 +18,66 @@
 __doctest_skip__ = ["grdfill"]
 
 
+def _parse_fill_mode(
+    constantfill=None, gridfill=None, neighborfill=None, splinefill=None
+) -> str | None:
+    """
+    Parse the fill parameters and return the appropriate string for the -A option.
+
+    >>> import numpy as np
+    >>> import xarray as xr
+    >>> _parse_fill_mode(constantfill=20.0)
+    'c20.0'
+    >>> _parse_fill_mode(gridfill="bggrid.nc")
+    'g'
+    >>> _parse_fill_mode(gridfill=xr.DataArray(np.zeros((10, 10))))
+    'g'
+    >>> _parse_fill_mode(neighborfill=20)
+    'n20'
+    >>> _parse_fill_mode(neighborfill=True)
+    'n'
+    >>> _parse_fill_mode(splinefill=0.5)
+    's0.5'
+    >>> _parse_fill_mode(splinefill=True)
+    's'
+    >>> _parse_fill_mode(constantfill=20, gridfill="bggrid.nc")
+    Traceback (most recent call last):
+    ...
+    pygmt.exceptions.GMTInvalidInput: The ... parameters are mutually exclusive.
+    """
+    fill_params = [constantfill, gridfill, neighborfill, splinefill]
+    if sum(param is not None for param in fill_params) > 1:
+        msg = (
+            "The 'constantfill', 'gridfill', 'neighborfill', and 'splinefill' "
+            "parameters are mutually exclusive."
+        )
+        raise GMTInvalidInput(msg)
+
+    if constantfill is not None:
+        return f"c{constantfill}"
+    if gridfill is not None:
+        return "g"  # Append grid file name later to support xarray.DataArray.
+    if neighborfill is not None and neighborfill is not False:
+        return "n" if neighborfill is True else f"n{neighborfill}"
+    if splinefill is not None and splinefill is not False:
+        return "s" if splinefill is True else f"s{splinefill}"
+    return None
+
+
 @fmt_docstring
 # TODO(PyGMT>=0.19.0): Remove the deprecated 'no_data' parameter.
 @deprecate_parameter("no_data", "hole", "v0.15.0", remove_version="v0.19.0")
 @use_alias(A="mode", N="hole", R="region", V="verbose")
 @kwargs_to_strings(R="sequence")
-def grdfill(grid, outgrid: str | None = None, **kwargs) -> xr.DataArray | None:
+def grdfill(
+    grid: str | xr.DataArray,
+    outgrid: str | None = None,
+    constantfill: float | None = None,
+    gridfill: str | xr.DataArray | None = None,
+    neighborfill: float | bool | None = None,
+    splinefill: float | bool | None = None,
+    **kwargs,
+) -> xr.DataArray | None:
     r"""
     Interpolate across holes in a grid.
 
@@ -31,25 +87,40 @@ def grdfill(grid, outgrid: str | None = None, **kwargs) -> xr.DataArray | None:
     replace the hole values. If no holes are found the original unchanged grid is
     returned.
 
-    Full option list at :gmt-docs:`grdfill.html`
+    Full option list at :gmt-docs:`grdfill.html`.
 
     {aliases}
 
     Parameters
     ----------
     {grid}
     {outgrid}
-    mode : str
-        Specify the hole-filling algorithm to use.  Choose from **c** for
-        constant fill and append the constant value, **n** for nearest
-        neighbor (and optionally append a search radius in
-        pixels [default radius is :math:`r^2 = \sqrt{{ X^2 + Y^2 }}`,
-        where (*X,Y*) are the node dimensions of the grid]), or
-        **s** for bicubic spline (optionally append a *tension*
-        parameter [Default is no tension]).
+    constantfill
+        Fill the holes with a constant value. Specify the constant value to use.
+    gridfill
+        Fill the holes with values sampled from another (possibly coarser) grid. Specify
+        the grid (a file name or an :class:`xarray.DataArray`) to use for the fill.
+    neighborfill
+        Fill the holes with the nearest neighbor. Specify the search radius in pixels.
+        If set to ``True``, the default search radius will be used
+        (:math:`r^2 = \sqrt{{n^2 + m^2}}`, where (*n,m*) are the node dimensions of the
+        grid).
+    splinefill
+        Fill the holes with a bicubic spline. Specify the tension value to use. If set
+        to ``True``, no tension will be used.
     hole : float
         Set the node value used to identify a point as a member of a hole [Default is
         NaN].
+    mode : str
+        Specify the hole-filling algorithm to use. Choose from **c** for constant fill
+        and append the constant value, **n** for nearest neighbor (and optionally append
+        a search radius in pixels [default radius is :math:`r^2 = \sqrt{{ X^2 + Y^2 }}`,
+        where (*X,Y*) are the node dimensions of the grid]), or **s** for bicubic spline
+        (optionally append a *tension* parameter [Default is no tension]).
+
+        .. deprecated:: 0.15.0
+            Use ``constantfill``, ``gridfill``, ``neighborfill``, or ``splinefill``
+            instead. The parameter will be removed in v0.19.0.
     {region}
     {verbose}
 
@@ -67,19 +138,37 @@ def grdfill(grid, outgrid: str | None = None, **kwargs) -> xr.DataArray | None:
     >>> import pygmt
     >>> # Load a bathymetric grid with missing data
     >>> earth_relief_holes = pygmt.datasets.load_sample_data(name="earth_relief_holes")
-    >>> # Perform grid filling operations on the sample grid
-    >>> # Set all empty values to "20"
-    >>> filled_grid = pygmt.grdfill(grid=earth_relief_holes, mode="c20")
+    >>> # Fill the holes with a constant value of 20
+    >>> filled_grid = pygmt.grdfill(grid=earth_relief_holes, constantfill=20)
     """
+    # TODO(PyGMT>=0.19.0): Remove the deprecated 'mode' parameter.
+    if kwargs.get("A") is not None:  # The deprecated 'mode' parameter is given.
+        warnings.warn(
+            "The 'mode' parameter is deprecated since v0.15.0 and will be removed in "
+            "v0.19.0. Use 'constantfill'/'gridfill'/'neighborfill'/'splinefill' "
+            "instead.",
+            FutureWarning,
+            stacklevel=1,
+        )
+    else:
+        # Determine the -A option from the fill parameters.
+        kwargs["A"] = _parse_fill_mode(constantfill, gridfill, neighborfill, splinefill)
+
     if kwargs.get("A") is None and kwargs.get("L") is None:
         msg = "At least parameter 'mode' or 'L' must be specified."
         raise GMTInvalidInput(msg)
 
     with Session() as lib:
         with (
             lib.virtualfile_in(check_kind="raster", data=grid) as vingrd,
+            lib.virtualfile_in(
+                check_kind="raster", data=gridfill, required_data=False
+            ) as vbggrd,
             lib.virtualfile_out(kind="grid", fname=outgrid) as voutgrd,
         ):
+            if gridfill is not None:
+                # Fill by a grid. Append the actual or virtual grid file name.
+                kwargs["A"] = f"g{vbggrd}"
             kwargs["G"] = voutgrd
             lib.call_module(
                 module="grdfill", args=build_arg_list(kwargs, infile=vingrd)

pygmt/tests/test_grdfill.py

@@ -5,6 +5,7 @@
 from pathlib import Path
 
 import numpy as np
+import numpy.testing as npt
 import pytest
 import xarray as xr
 from pygmt import grdfill, load_dataarray
@@ -49,23 +50,8 @@ def fixture_expected_grid():
             [347.5, 331.5, 309.0, 282.0, 190.0, 208.0, 299.5, 348.0],
         ],
         coords={
-            "lon": [-54.5, -53.5, -52.5, -51.5, -50.5, -49.5, -48.5, -47.5],
-            "lat": [
-                -23.5,
-                -22.5,
-                -21.5,
-                -20.5,
-                -19.5,
-                -18.5,
-                -17.5,
-                -16.5,
-                -15.5,
-                -14.5,
-                -13.5,
-                -12.5,
-                -11.5,
-                -10.5,
-            ],
+            "lon": np.arange(-54.5, -46.5, 1),
+            "lat": np.arange(-23.5, -9.5, 1),
         },
         dims=["lat", "lon"],
     )
@@ -76,7 +62,7 @@ def test_grdfill_dataarray_out(grid, expected_grid):
     """
     Test grdfill with a DataArray output.
     """
-    result = grdfill(grid=grid, mode="c20")
+    result = grdfill(grid=grid, constantfill=20)
     # check information of the output grid
     assert isinstance(result, xr.DataArray)
     assert result.gmt.gtype == GridType.GEOGRAPHIC
@@ -91,7 +77,7 @@ def test_grdfill_asymmetric_pad(grid, expected_grid):
 
     Regression test for https://github.com/GenericMappingTools/pygmt/issues/1745.
     """
-    result = grdfill(grid=grid, mode="c20", region=[-55, -50, -24, -16])
+    result = grdfill(grid=grid, constantfill=20, region=[-55, -50, -24, -16])
     # check information of the output grid
     assert isinstance(result, xr.DataArray)
     assert result.gmt.gtype == GridType.GEOGRAPHIC
@@ -107,16 +93,40 @@ def test_grdfill_file_out(grid, expected_grid):
     Test grdfill with an outgrid set.
     """
     with GMTTempFile(suffix=".nc") as tmpfile:
-        result = grdfill(grid=grid, mode="c20", outgrid=tmpfile.name)
+        result = grdfill(grid=grid, constantfill=20, outgrid=tmpfile.name)
         assert result is None  # return value is None
         assert Path(tmpfile.name).stat().st_size > 0  # check that outfile exists
         temp_grid = load_dataarray(tmpfile.name)
         xr.testing.assert_allclose(a=temp_grid, b=expected_grid)
 
 
+def test_grdfill_gridfill_dataarray(grid):
+    """
+    Test grdfill with a DataArray input.
+    """
+    bggrid = xr.DataArray(
+        np.arange(grid.size).reshape(grid.shape),
+        dims=grid.dims,
+        coords={"lon": grid.lon, "lat": grid.lat},
+    )
+    result = grdfill(grid=grid, gridfill=bggrid)
+    assert not result.isnull().any()
+    npt.assert_array_equal(result[3:6, 3:5], bggrid[3:6, 3:5])
+
+
 def test_grdfill_required_args(grid):
     """
     Test that grdfill fails without arguments for `mode` and `L`.
     """
     with pytest.raises(GMTInvalidInput):
         grdfill(grid=grid)
+
+
+# TODO(PyGMT>=0.19.0): Remove this test.
+def test_grdfill_deprecated_mode(grid, expected_grid):
+    """
+    Test that grdfill fails with deprecated `mode` argument.
+    """
+    with pytest.warns(FutureWarning):
+        result = grdfill(grid=grid, mode="c20")
+        xr.testing.assert_allclose(a=result, b=expected_grid)