Skip to content

Commit f910a2a

Browse files
seismanseisman
committed
pygmt.grdfilter: Add parameters filter_type/filter_width and more to set the filter
1 parent 83f06c3 commit f910a2a

2 files changed

Lines changed: 220 additions & 19 deletions

File tree

‎pygmt/src/grdfilter.py‎

Lines changed: 212 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -9,16 +9,139 @@
99
from pygmt._typing import PathLike
1010
from pygmt.alias import Alias, AliasSystem
1111
from pygmt.clib import Session
12+
from pygmt.exceptions import GMTParameterError
1213
from pygmt.helpers import build_arg_list, fmt_docstring, use_alias
1314

1415
__doctest_skip__ = ["grdfilter"]
1516

1617

18+
def _alias_option_F( # noqa: N802
19+
filter_type=None,
20+
filter_width=None,
21+
hist_bin_width=None,
22+
highpass=False,
23+
median_quantile=None,
24+
hist_center_bins=False,
25+
mode_extreme=None,
26+
filter=None, # noqa: A002
27+
):
28+
"""
29+
Helper function to create the alias list for the -F option.
30+
31+
Examples
32+
--------
33+
>>> def parse(**kwargs):
34+
... return AliasSystem(F=_alias_option_F(**kwargs)).get("F")
35+
>>> parse(filter_type="boxcar", filter_width=2.0)
36+
'b2.0'
37+
>>> parse(filter_type="cosine_arch", filter_width=(5, 10))
38+
'c5/10'
39+
>>> parse(filter_type="gaussian", filter_width=100, highpass=True)
40+
'g100+h'
41+
>>> parse(filter_type="median", median_quantile=0.25)
42+
'm+q0.25'
43+
>>> parse(
44+
... filter_type="histogram",
45+
... filter_width=100,
46+
... hist_bin_width=1.0,
47+
... hist_center_bins=True,
48+
... mode_extreme="max",
49+
... )
50+
'h100/1.0+c+u'
51+
"""
52+
if filter is not None:
53+
kwdict = {
54+
"filter_type": filter_type,
55+
"filter_width": filter_width,
56+
"hist_bin_width": hist_bin_width,
57+
"highpass": highpass,
58+
"median_quantile": median_quantile,
59+
"hist_center_bins": hist_center_bins,
60+
"mode_extreme": mode_extreme,
61+
}
62+
if any(v is not None and v is not False for v in kwdict.values()):
63+
raise GMTParameterError(
64+
conflicts_with=("filter", kwdict.keys()),
65+
reason="'filter' is specified using the unrecommended GMT command string syntax.",
66+
)
67+
return Alias(filter, name="filter") # Deprecated raw GMT string.
68+
69+
if median_quantile is not None and filter_type != "median":
70+
raise GMTParameterError(
71+
conflicts_with=("median_quantile", [f"filter_type={filter_type!r}"]),
72+
reason="'median_quantile' is allowed only when 'filter_type' is 'median'.",
73+
)
74+
if hist_bin_width is not None and filter_type != "histogram":
75+
raise GMTParameterError(
76+
conflicts_with=("hist_bin_width", [f"filter_type={filter_type!r}"]),
77+
reason="'hist_bin_width' is allowed only when 'filter_type' is 'histogram'.",
78+
)
79+
if hist_center_bins is not False and filter_type != "histogram":
80+
raise GMTParameterError(
81+
conflicts_with=("hist_center_bins", [f"filter_type={filter_type!r}"]),
82+
reason="'hist_center_bins' is allowed only when 'filter_type' is 'histogram'.",
83+
)
84+
if mode_extreme is not None and filter_type not in {"mlprob", "histogram"}:
85+
raise GMTParameterError(
86+
conflicts_with=("mode_extreme", [f"filter_type={filter_type!r}"]),
87+
reason="'mode_extreme' is allowed only when 'filter_type' is 'mlprob' or 'histogram'.",
88+
)
89+
90+
return [
91+
Alias(
92+
filter_type,
93+
name="filter_type",
94+
mapping={
95+
"boxcar": "b",
96+
"cosine_arch": "c",
97+
"gaussian": "g",
98+
"custom": "f",
99+
"operator": "o",
100+
"median": "m",
101+
"mlprob": "p",
102+
"histogram": "h",
103+
"minall": "l",
104+
"minpos": "L",
105+
"maxall": "u",
106+
"maxneg": "U",
107+
},
108+
),
109+
Alias(filter_width, name="filter_width", sep="/"),
110+
Alias(hist_bin_width, name="hist_bin_width", prefix="/"),
111+
Alias(hist_center_bins, name="hist_center_bins", prefix="+c"),
112+
Alias(highpass, name="highpass", prefix="+h"),
113+
Alias(median_quantile, name="median_quantile", prefix="+q"),
114+
Alias(mode_extreme, name="mode_extreme", mapping={"min": "+l", "max": "+u"}),
115+
]
116+
117+
17118
@fmt_docstring
18-
@use_alias(D="distance", F="filter", f="coltypes")
19-
def grdfilter(
119+
@use_alias(D="distance", f="coltypes")
120+
def grdfilter( # noqa: PLR0913
20121
grid: PathLike | xr.DataArray,
21122
outgrid: PathLike | None = None,
123+
filter_type: Literal[
124+
"boxcar",
125+
"cosine_arch",
126+
"gaussian",
127+
"custom",
128+
"operator",
129+
"median",
130+
"mlprob",
131+
"histogram",
132+
"minall",
133+
"minpos",
134+
"maxall",
135+
"maxneg",
136+
]
137+
| None = None,
138+
filter_width: Sequence[float] | None = None,
139+
highpass: bool = False,
140+
median_quantile: float | None = None,
141+
hist_bin_width: float | None = None,
142+
hist_center_bins: bool = False,
143+
mode_extreme: Literal["min", "max"] | None = None,
144+
filter: str | None = None, # noqa: A002
22145
spacing: Sequence[float | str] | None = None,
23146
nans: Literal["ignore", "replace", "preserve"] | None = None,
24147
toggle: bool = False,
@@ -29,7 +152,7 @@ def grdfilter(
29152
cores: int | bool = False,
30153
**kwargs,
31154
) -> xr.DataArray | None:
32-
r"""
155+
"""
33156
Filter a grid in the space (or time) domain.
34157
35158
Filter a grid file in the space (or time) domain using one of the selected
@@ -58,19 +181,78 @@ def grdfilter(
58181
----------
59182
$grid
60183
$outgrid
61-
filter : str
62-
**b**\|\ **c**\|\ **g**\|\ **o**\|\ **m**\|\ **p**\|\ **h**\ *width*\
63-
[/*width2*\][*modifiers*].
64-
Name of the filter type you wish to apply, followed by the *width*:
65-
66-
- **b**: Box Car
67-
- **c**: Cosine Arch
68-
- **g**: Gaussian
69-
- **o**: Operator
70-
- **m**: Median
71-
- **p**: Maximum Likelihood probability
72-
- **h**: Histogram
184+
filter_type
185+
The filter type. Choose among convolution and non-convolution filters.
186+
187+
Convolution filters include:
188+
189+
- ``"boxcar"``: All weights are equal.
190+
- ``"cosine_arch"``: Weights follow a cosine arch curve.
191+
- ``"gaussian"``: Weights are given by the Gaussian function, where filter width
192+
is 6 times the conventional Gaussian sigma.
193+
- ``"custom"``: Weights are given by the precomputed values in the filter weight
194+
grid file *weight*, which must have odd dimensions; also requires ``distance=0``
195+
and output spacing must match input spacing or be integer multiples.
196+
- ``"operator"``: Weights are given by the precomputed values in the filter weight
197+
grid file *weight*, which must have odd dimensions; also requires ``distance=0``
198+
and output spacing must match input spacing or be integer multiples. Weights
199+
are assumed to sum to zero so no accumulation of weight sums and normalization
200+
will be done.
201+
202+
Non-convolution filters include:
73203
204+
- ``"median"``: Returns median value. To select another quantile, use the
205+
parameter ``median_quantile`` in the 0-1 range [Default is 0.5, i.e., median].
206+
- ``"mlprob"``: Maximum likelihood probability (a mode estimator). Return modal
207+
value. If more than one mode is found we return their average value. Set
208+
``mode_extreme`` to ``"min"`` or ``"max"`` to return the lowermost or uppermost
209+
of the modal values.
210+
- ``"histogram"``: Histogram mode (another mode estimator). Return the modal value
211+
as the center of the dominant peak in a histogram. Use parameter
212+
``histogram_center_bins`` to center the bins on multiples of bin width [Default
213+
has bin edges that are multiples of bin width]. Use parameter
214+
``histogram_bin_width`` to set the bin width. If more than one mode is found we
215+
return their average value. Set ``mode_extreme`` to ``"min"`` or ``"max"`` to
216+
return the lowermost or uppermost of the modal values.
217+
218+
- ``"minall"``: Return minimum of all values.
219+
- ``"minpos"``: Return minimum of all positive values only.
220+
- ``"maxall"``: Return maximum of all values.
221+
- ``"maxneg"``: Return maximum of all negative values only.
222+
filter_width
223+
The full diameter width of the filter. It can be a single value for an isotropic
224+
filter, or a pair of values for a rectangular filter (width in x- and
225+
y-directions, requiring ``distance`` be either ``"p"`` or ``0``). For isotropic
226+
filters, ``width`` can also be a path to a grid file for variable filter width,
227+
in which case the grid must have the same registration and dimensions as the
228+
output filtered grid.
229+
highpass
230+
By default, the filter is a low-pass filter. If True, then the filter is a
231+
high-pass filter. [Default is ``False``].
232+
median_quantile
233+
Quantile to use when ``filter_type="median"``. Must be a float in the range 0-1.
234+
[Default is 0.5 (median)].
235+
hist_bin_width
236+
Bin width to use when ``filter_type="histogram"``.
237+
hist_center_bins
238+
Center the histogram bins on multiples of *histogram_bin_width* when
239+
``filter_type="histogram"``. By default, the bins are aligned such that
240+
their edges are on multiples of *hist_bin_width*.
241+
mode_extreme
242+
Choose which extreme to return when ``filter_type="mlprob"`` or
243+
``filter_type="histogram"`` and multiple modes are found. Options are: ``"min"``
244+
to return the lowermost mode, or ``"max"`` to return the uppermost mode. By
245+
default, the average of all modes is returned.
246+
filter
247+
Set the filter type.
248+
249+
.. deprecated:: v0.19.0
250+
251+
This parameter is deprecated. Use the parameters ``filter_type``,
252+
``filter_width``, ``hist_bin_width``, ``highpass``, ``median_quantile``,
253+
``hist_center_bins``, and ``mode_extreme`` instead. This parameter still
254+
accepts raw GMT CLI strings for the ``-F`` option of the ``grdfilter``
255+
module for backward compatibility.
74256
distance : str
75257
State how the grid (x,y) relates to the filter *width*:
76258
@@ -130,7 +312,8 @@ def grdfilter(
130312
>>> # and return a filtered grid (saved as netCDF file).
131313
>>> pygmt.grdfilter(
132314
... grid="@earth_relief_30m_g",
133-
... filter="m600",
315+
... filter_type="median",
316+
... filter_width=600,
134317
... distance="4",
135318
... region=[150, 250, 10, 40],
136319
... spacing=0.5,
@@ -140,9 +323,21 @@ def grdfilter(
140323
>>> # Apply a Gaussian smoothing filter of 600 km to the input DataArray and return
141324
>>> # a filtered DataArray with the smoothed grid.
142325
>>> grid = pygmt.datasets.load_earth_relief()
143-
>>> smooth_field = pygmt.grdfilter(grid=grid, filter="g600", distance="4")
326+
>>> smooth_field = pygmt.grdfilter(
327+
... grid=grid, filter_type="gaussian", filter_width=600, distance="4"
328+
... )
144329
"""
145330
aliasdict = AliasSystem(
331+
F=_alias_option_F(
332+
filter_type=filter_type,
333+
filter_width=filter_width,
334+
hist_bin_width=hist_bin_width,
335+
highpass=highpass,
336+
median_quantile=median_quantile,
337+
hist_center_bins=hist_center_bins,
338+
mode_extreme=mode_extreme,
339+
filter=filter,
340+
),
146341
I=Alias(spacing, name="spacing", sep="/", size=2),
147342
N=Alias(
148343
nans, name="nans", mapping={"ignore": "i", "replace": "r", "preserve": "p"}

‎pygmt/tests/test_grdfilter.py‎

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,12 @@ def test_grdfilter_dataarray_in_dataarray_out(grid, expected_grid):
4747
Test grdfilter with an input DataArray, and output as DataArray.
4848
"""
4949
result = grdfilter(
50-
grid=grid, filter="g600", distance="4", region=[-53, -49, -20, -17], cores=2
50+
grid=grid,
51+
filter_type="gaussian",
52+
filter_width=600,
53+
distance="4",
54+
region=[-53, -49, -20, -17],
55+
cores=2,
5156
)
5257
# check information of the output grid
5358
assert isinstance(result, xr.DataArray)
@@ -65,7 +70,8 @@ def test_grdfilter_dataarray_in_file_out(grid, expected_grid):
6570
result = grdfilter(
6671
grid,
6772
outgrid=tmpfile.name,
68-
filter="g600",
73+
filter_type="gaussian",
74+
filter_width=600,
6975
distance="4",
7076
region=[-53, -49, -20, -17],
7177
)

0 commit comments

Comments
 (0)