diff --git a/pygmt/src/grdfilter.py b/pygmt/src/grdfilter.py index 285e4848375..352e40e0673 100644 --- a/pygmt/src/grdfilter.py +++ b/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"} diff --git a/pygmt/tests/test_grdfilter.py b/pygmt/tests/test_grdfilter.py index 15b9a977b3e..b318da9ba8e 100644 --- a/pygmt/tests/test_grdfilter.py +++ b/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,7 +87,10 @@ 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", ) @@ -93,5 +98,24 @@ 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)