Skip to content

pygmt.grdfilter: Add parameters filter_type/filter_width and more to control filter parameters #4378

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
seisman wants to merge 1 commit into
base: main
Choose a base branch
from
+220 −19
Draft

pygmt.grdfilter: Add parameters filter_type/filter_width and more to control filter parameters #4378

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
229 changes: 212 additions & 17 deletions pygmt/src/grdfilter.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,16 +9,139 @@
from pygmt._typing import PathLike
from pygmt.alias import Alias, AliasSystem
from pygmt.clib import Session
from pygmt.exceptions import GMTParameterError
from pygmt.helpers import build_arg_list, fmt_docstring, use_alias

__doctest_skip__ = ["grdfilter"]


def _alias_option_F( # noqa: N802
filter_type=None,
filter_width=None,
hist_bin_width=None,
highpass=False,
median_quantile=None,
hist_center_bins=False,
mode_extreme=None,
filter=None, # noqa: A002
):
"""
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_type="boxcar", filter_width=2.0)
'b2.0'
>>> parse(filter_type="cosine_arch", filter_width=(5, 10))
'c5/10'
>>> parse(filter_type="gaussian", filter_width=100, highpass=True)
'g100+h'
>>> parse(filter_type="median", median_quantile=0.25)
'm+q0.25'
>>> parse(
... filter_type="histogram",
... filter_width=100,
... hist_bin_width=1.0,
... hist_center_bins=True,
... mode_extreme="max",
... )
'h100/1.0+c+u'
"""
if filter is not None:
kwdict = {
"filter_type": filter_type,
"filter_width": filter_width,
"hist_bin_width": hist_bin_width,
"highpass": highpass,
"median_quantile": median_quantile,
"hist_center_bins": hist_center_bins,
"mode_extreme": mode_extreme,
}
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 string.

if median_quantile is not None and filter_type != "median":
raise GMTParameterError(
conflicts_with=("median_quantile", [f"filter_type={filter_type!r}"]),
reason="'median_quantile' is allowed only when 'filter_type' is 'median'.",
)
if hist_bin_width is not None and filter_type != "histogram":
raise GMTParameterError(
conflicts_with=("hist_bin_width", [f"filter_type={filter_type!r}"]),
reason="'hist_bin_width' is allowed only when 'filter_type' is 'histogram'.",
)
if hist_center_bins is not False and filter_type != "histogram":
raise GMTParameterError(
conflicts_with=("hist_center_bins", [f"filter_type={filter_type!r}"]),
reason="'hist_center_bins' is allowed only when 'filter_type' is 'histogram'.",
)
if mode_extreme is not None and filter_type not in {"mlprob", "histogram"}:
raise GMTParameterError(
conflicts_with=("mode_extreme", [f"filter_type={filter_type!r}"]),
reason="'mode_extreme' is allowed only when 'filter_type' is 'mlprob' or 'histogram'.",
)

return [
Alias(
filter_type,
name="filter_type",
mapping={
"boxcar": "b",
"cosine_arch": "c",
"gaussian": "g",
"custom": "f",
"operator": "o",
"median": "m",
"mlprob": "p",
"histogram": "h",
"minall": "l",
"minpos": "L",
"maxall": "u",
"maxneg": "U",
},
),
Alias(filter_width, name="filter_width", sep="/"),
Alias(hist_bin_width, name="hist_bin_width", prefix="/"),
Alias(hist_center_bins, name="hist_center_bins", prefix="+c"),
Alias(highpass, name="highpass", prefix="+h"),
Alias(median_quantile, name="median_quantile", prefix="+q"),
Alias(mode_extreme, name="mode_extreme", mapping={"min": "+l", "max": "+u"}),
]


@fmt_docstring
@use_alias(D="distance", F="filter", f="coltypes")
def grdfilter(
@use_alias(D="distance", f="coltypes")
def grdfilter( # noqa: PLR0913
grid: PathLike | xr.DataArray,
outgrid: PathLike | None = None,
filter_type: Literal[
"boxcar",
"cosine_arch",
"gaussian",
"custom",
"operator",
"median",
"mlprob",
"histogram",
"minall",
"minpos",
"maxall",
"maxneg",
]
| None = None,
filter_width: Sequence[float] | None = None,
highpass: bool = False,
median_quantile: float | None = None,
hist_bin_width: float | None = None,
hist_center_bins: bool = False,
mode_extreme: Literal["min", "max"] | None = None,
filter: str | None = None, # noqa: A002
spacing: Sequence[float | str] | None = None,
nans: Literal["ignore", "replace", "preserve"] | None = None,
toggle: bool = False,
Expand All @@ -29,7 +152,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
Expand Down Expand Up @@ -58,19 +181,78 @@ 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_type
The filter type. Choose among convolution and non-convolution filters.

Convolution filters include:

- ``"boxcar"``: All weights are equal.
- ``"cosine_arch"``: Weights follow a cosine arch curve.
- ``"gaussian"``: Weights are given by the Gaussian function, where filter width
is 6 times the conventional Gaussian sigma.
- ``"custom"``: Weights are given by the precomputed values in the filter weight
grid file *weight*, which must have odd dimensions; also requires ``distance=0``
and output spacing must match input spacing or be integer multiples.
- ``"operator"``: Weights are given by the precomputed values in the filter weight
grid file *weight*, which must have odd dimensions; also requires ``distance=0``
and output spacing must match input spacing or be integer multiples. Weights
are assumed to sum to zero so no accumulation of weight sums and normalization
will be done.

Non-convolution filters include:

- ``"median"``: Returns median value. To select another quantile, use the
parameter ``median_quantile`` in the 0-1 range [Default is 0.5, i.e., median].
- ``"mlprob"``: Maximum likelihood probability (a mode estimator). Return modal
value. If more than one mode is found we return their average value. Set
``mode_extreme`` to ``"min"`` or ``"max"`` to return the lowermost or uppermost
of the modal values.
- ``"histogram"``: Histogram mode (another mode estimator). Return the modal value
as the center of the dominant peak in a histogram. Use parameter
``histogram_center_bins`` to center the bins on multiples of bin width [Default
has bin edges that are multiples of bin width]. Use parameter
``histogram_bin_width`` to set the bin width. If more than one mode is found we
return their average value. Set ``mode_extreme`` to ``"min"`` or ``"max"`` to
return the lowermost or uppermost of the modal values.

- ``"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.
filter_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 ``"p"`` or ``0``). For isotropic
filters, ``width`` can also be a path to a grid file for variable filter width,
in which case the grid must have the same registration and dimensions as the
output filtered grid.
highpass
By default, the filter is a low-pass filter. If True, then the filter is a
high-pass filter. [Default is ``False``].
median_quantile
Quantile to use when ``filter_type="median"``. Must be a float in the range 0-1.
[Default is 0.5 (median)].
hist_bin_width
Bin width to use when ``filter_type="histogram"``.
hist_center_bins
Center the histogram bins on multiples of *histogram_bin_width* when
``filter_type="histogram"``. By default, the bins are aligned such that
their edges are on multiples of *hist_bin_width*.
mode_extreme
Choose which extreme to return when ``filter_type="mlprob"`` or
``filter_type="histogram"`` and multiple modes are found. Options are: ``"min"``
to return the lowermost mode, or ``"max"`` to return the uppermost mode. By
default, the average of all modes is returned.
filter
Set the filter type.

.. deprecated:: v0.19.0

This parameter is deprecated. Use the parameters ``filter_type``,
``filter_width``, ``hist_bin_width``, ``highpass``, ``median_quantile``,
``hist_center_bins``, and ``mode_extreme`` instead. This parameter still
accepts raw GMT CLI strings for the ``-F`` option of the ``grdfilter``
module for backward compatibility.
distance : str
State how the grid (x,y) relates to the filter *width*:

Expand Down Expand Up @@ -130,7 +312,8 @@ def grdfilter(
>>> # and return a filtered grid (saved as netCDF file).
>>> pygmt.grdfilter(
... grid="@earth_relief_30m_g",
... filter="m600",
... filter_type="median",
... filter_width=600,
... distance="4",
... region=[150, 250, 10, 40],
... spacing=0.5,
Expand All @@ -140,9 +323,21 @@ 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()
>>> smooth_field = pygmt.grdfilter(grid=grid, filter="g600", distance="4")
>>> smooth_field = pygmt.grdfilter(
... grid=grid, filter_type="gaussian", filter_width=600, distance="4"
... )
"""
aliasdict = AliasSystem(
F=_alias_option_F(
filter_type=filter_type,
filter_width=filter_width,
hist_bin_width=hist_bin_width,
highpass=highpass,
median_quantile=median_quantile,
hist_center_bins=hist_center_bins,
mode_extreme=mode_extreme,
filter=filter,
),
I=Alias(spacing, name="spacing", sep="/", size=2),
N=Alias(
nans, name="nans", mapping={"ignore": "i", "replace": "r", "preserve": "p"}
Expand Down
10 changes: 8 additions & 2 deletions pygmt/tests/test_grdfilter.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,12 @@ def test_grdfilter_dataarray_in_dataarray_out(grid, expected_grid):
Test grdfilter with an input DataArray, and output as DataArray.
"""
result = grdfilter(
grid=grid, filter="g600", distance="4", region=[-53, -49, -20, -17], cores=2
grid=grid,
filter_type="gaussian",
filter_width=600,
distance="4",
region=[-53, -49, -20, -17],
cores=2,
)
# check information of the output grid
assert isinstance(result, xr.DataArray)
Expand All @@ -65,7 +70,8 @@ def test_grdfilter_dataarray_in_file_out(grid, expected_grid):
result = grdfilter(
grid,
outgrid=tmpfile.name,
filter="g600",
filter_type="gaussian",
filter_width=600,
distance="4",
region=[-53, -49, -20, -17],
)
Expand Down