pygmt.grdfilter: Add parameters filter/width/highpass to set the filter [Part 1] (#4401)

Author: seisman

pygmt/src/grdfilter.py

@@ -15,11 +15,68 @@
 __doctest_skip__ = ["grdfilter"]
 
 
+def _alias_option_F(  # noqa: N802
+    filter=None,  # noqa: A002
+    width=None,
+    highpass=False,
+):
+    """
+    Helper function to create the alias list for the -F option.
+
+    Examples
+    --------
+    >>> def parse(**kwargs):
+    ...     return AliasSystem(F=_alias_option_F(**kwargs)).get("F")
+    >>> parse(filter="boxcar", width=2.0)
+    'b2.0'
+    >>> parse(filter="cosarch", width=(5, 10))
+    'c5/10'
+    >>> parse(filter="gaussian", width=100, highpass=True)
+    'g100+h'
+    """
+    _filter_mapping = {
+        "boxcar": "b",
+        "cosarch": "c",
+        "gaussian": "g",
+        "minall": "l",
+        "minpos": "L",
+        "maxall": "u",
+        "maxneg": "U",
+    }
+    # Check if the 'filter' parameter is using the old GMT command string syntax.
+    _old_filter_syntax = isinstance(filter, str) and filter not in _filter_mapping
+
+    if _old_filter_syntax:
+        kwdict = {"width": width, "highpass": highpass}
+        if any(v is not None and v is not False for v in kwdict.values()):
+            raise GMTParameterError(
+                conflicts_with=("filter", kwdict.keys()),
+                reason="'filter' is specified using the unrecommended GMT command string syntax.",
+            )
+        return Alias(filter, name="filter")  # Deprecated raw GMT command string.
+
+    if filter is None or width is None:
+        raise GMTParameterError(required=["filter", "width"])
+
+    return [
+        Alias(filter, name="filter", mapping=_filter_mapping),
+        Alias(width, name="width", sep="/", size=2),
+        Alias(highpass, name="highpass", prefix="+h"),
+    ]
+
+
 @fmt_docstring
-@use_alias(F="filter", f="coltypes")
-def grdfilter(
+@use_alias(f="coltypes")
+def grdfilter(  # noqa: PLR0913
     grid: PathLike | xr.DataArray,
     outgrid: PathLike | None = None,
+    filter: Literal[  # noqa: A002
+        "boxcar", "cosarch", "gaussian", "minall", "minpos", "maxall", "maxneg"
+    ]
+    | str
+    | None = None,
+    width: float | Sequence[float] | None = None,
+    highpass: bool = False,
     distance: Literal[
         "pixel",
         "cartesian",
@@ -40,7 +97,7 @@ def grdfilter(
     cores: int | bool = False,
     **kwargs,
 ) -> xr.DataArray | None:
-    r"""
+    """
     Filter a grid in the space (or time) domain.
 
     Filter a grid file in the space (or time) domain using one of the selected
@@ -57,6 +114,7 @@ def grdfilter(
 
     $aliases
        - D = distance
+       - F = filter, width, **+h**: highpass
        - G = outgrid
        - I = spacing
        - N = nans
@@ -70,19 +128,35 @@ def grdfilter(
     ----------
     $grid
     $outgrid
-    filter : str
-        **b**\|\ **c**\|\ **g**\|\ **o**\|\ **m**\|\ **p**\|\ **h**\ *width*\
-        [/*width2*\][*modifiers*].
-        Name of the filter type you wish to apply, followed by the *width*:
-
-        - **b**: Box Car
-        - **c**: Cosine Arch
-        - **g**: Gaussian
-        - **o**: Operator
-        - **m**: Median
-        - **p**: Maximum Likelihood probability
-        - **h**: Histogram
+    filter
+        The filter type. Choose among convolution and non-convolution filters.
+
+        Convolution filters include:
 
+        - ``"boxcar"``: All weights are equal.
+        - ``"cosarch"``: Weights follow a cosine arch curve.
+        - ``"gaussian"``: Weights are given by the Gaussian function, where filter width
+          is 6 times the conventional Gaussian sigma.
+
+        Non-convolution filters include:
+
+        - ``"minall"``: Return minimum of all values.
+        - ``"minpos"``: Return minimum of all positive values only.
+        - ``"maxall"``: Return maximum of all values.
+        - ``"maxneg"``: Return maximum of all negative values only.
+
+        **Note**: There are still a few other filter types available in GMT (e.g.,
+        histogram and mode filters), but they are not implemented in PyGMT yet. As a
+        workaround, pass the raw GMT command string to this parameter to use these other
+        filter types. Refer to :gmt-docs:`grdfilter.html#f` for the full syntax of this
+        parameter.
+    width
+        The full diameter width of the filter. It can be a single value for an isotropic
+        filter, or a pair of values for a rectangular filter (width in x- and
+        y-directions, requiring ``distance`` be either ``"pixel"`` or ``"cartesian"``).
+    highpass
+        By default, the filter is a low-pass filter. If True, then the filter is a
+        high-pass filter. [Default is ``False``].
     distance
         Determine how grid (*x, y*) relates to filter *width* and how distances are
         calculated. Valid values are list below. The first four options are fastest
@@ -159,11 +233,12 @@ def grdfilter(
     --------
     >>> from pathlib import Path
     >>> import pygmt
-    >>> # Apply a median filter of 600 km (full width) to the @earth_relief_30m_g grid
+    >>> # Apply a gaussian filter of 600 km (full width) to the @earth_relief_30m_g grid
     >>> # and return a filtered grid (saved as netCDF file).
     >>> pygmt.grdfilter(
     ...     grid="@earth_relief_30m_g",
-    ...     filter="m600",
+    ...     filter="gaussian",
+    ...     width=600,
     ...     distance="geo_spherical",
     ...     region=[150, 250, 10, 40],
     ...     spacing=0.5,
@@ -173,10 +248,13 @@ def grdfilter(
     >>> # Apply a Gaussian smoothing filter of 600 km to the input DataArray and return
     >>> # a filtered DataArray with the smoothed grid.
     >>> grid = pygmt.datasets.load_earth_relief()
-    >>> smoothed = pygmt.grdfilter(grid=grid, filter="g600", distance="geo_spherical")
+    >>> smoothed = pygmt.grdfilter(
+    ...     grid=grid, filter="gaussian", width=600, distance="geo_spherical"
+    ... )
     """
     if kwargs.get("D", distance) is None:
         raise GMTParameterError(required="distance")
+    # filter is also required but will be checked in _alias_option_F.
 
     aliasdict = AliasSystem(
         D=Alias(
@@ -192,6 +270,11 @@ def grdfilter(
                 "geo_mercator": 5,
             },
         ),
+        F=_alias_option_F(
+            filter=filter,
+            width=width,
+            highpass=highpass,
+        ),
         I=Alias(spacing, name="spacing", sep="/", size=2),
         N=Alias(
             nans, name="nans", mapping={"ignore": "i", "replace": "r", "preserve": "p"}

pygmt/tests/test_grdfilter.py

@@ -48,7 +48,8 @@ def test_grdfilter_dataarray_in_dataarray_out(grid, expected_grid):
     """
     result = grdfilter(
         grid=grid,
-        filter="g600",
+        filter="gaussian",
+        width=600,
         distance="geo_spherical",
         region=[-53, -49, -20, -17],
         cores=2,
@@ -69,7 +70,8 @@ def test_grdfilter_dataarray_in_file_out(grid, expected_grid):
         result = grdfilter(
             grid,
             outgrid=tmpfile.name,
-            filter="g600",
+            filter="gaussian",
+            width=600,
             distance="geo_spherical",
             region=[-53, -49, -20, -17],
         )
@@ -85,13 +87,35 @@ def test_grdfilter_fails():
     """
     with pytest.raises(GMTTypeError):
         grdfilter(
-            np.arange(10).reshape((5, 2)), filter="g600", distance="geo_spherical"
+            np.arange(10).reshape((5, 2)),
+            filter="gaussian",
+            width=600,
+            distance="geo_spherical",
         )
 
 
 def test_grdfilter_required(grid):
     """
     Test that grdfilter raises an exception when required parameters are missing.
     """
+    with pytest.raises(GMTParameterError, match="distance"):
+        grdfilter(grid=grid)
+    with pytest.raises(GMTParameterError, match="distance"):
+        grdfilter(grid=grid, filter="gaussian")
+    with pytest.raises(GMTParameterError, match="width"):
+        grdfilter(grid=grid, filter="gaussian", distance="geo_spherical")
+    with pytest.raises(GMTParameterError, match="filter"):
+        grdfilter(grid=grid, width=600, distance="geo_spherical")
     with pytest.raises(GMTParameterError, match="distance"):
         grdfilter(grid=grid, filter="g600")
+
+
+def test_grdfilter_mixed_syntax(grid):
+    """
+    Test grdfilter's filter parameter with mixed syntax.
+    """
+    kwargs = {"grid": grid, "filter": "g600", "distance": 4}
+    with pytest.raises(GMTParameterError):
+        grdfilter(width=600, **kwargs)
+    with pytest.raises(GMTParameterError):
+        grdfilter(highpass=True, **kwargs)