diff --git a/pygmt/src/grdfill.py b/pygmt/src/grdfill.py index 7188ea14497..a6f489c8121 100644 --- a/pygmt/src/grdfill.py +++ b/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,7 +87,7 @@ 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} @@ -39,17 +95,32 @@ def grdfill(grid, outgrid: str | None = None, **kwargs) -> xr.DataArray | None: ---------- {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,10 +138,22 @@ 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) @@ -78,8 +161,14 @@ def grdfill(grid, outgrid: str | None = None, **kwargs) -> xr.DataArray | None: 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) diff --git a/pygmt/tests/test_grdfill.py b/pygmt/tests/test_grdfill.py index 8d36125c279..8763a482383 100644 --- a/pygmt/tests/test_grdfill.py +++ b/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)