"""
dimfilter - Directional filtering of grids in the space domain.
"""
from pygmt.clib import Session
from pygmt.exceptions import GMTInvalidInput
from pygmt.helpers import (
GMTTempFile,
build_arg_string,
fmt_docstring,
kwargs_to_strings,
use_alias,
)
from pygmt.io import load_dataarray
__doctest_skip__ = ["dimfilter"]
[docs]
@fmt_docstring
@use_alias(
D="distance",
F="filter",
G="outgrid",
I="spacing",
N="sectors",
R="region",
V="verbose",
)
@kwargs_to_strings(I="sequence", R="sequence")
def dimfilter(grid, **kwargs):
r"""
Filter a grid by dividing the filter circle.
Filter a grid in the space (or time) domain by
dividing the given filter circle into the given number of sectors,
applying one of the selected primary convolution or non-convolution
filters to each sector, and choosing the final outcome according to the
selected secondary filter. It computes distances using Cartesian or
Spherical geometries. The output grid can optionally be generated as a
subregion of the input and/or with a new increment using ``spacing``,
which may add an "extra space" in the input data to prevent edge
effects for the output grid. If the filter is low-pass, then the output
may be less frequently sampled than the input. :func:`pygmt.dimfilter`
will not produce a smooth output as other spatial filters
do because it returns a minimum median out of *N* medians of *N*
sectors. The output can be rough unless the input data are noise-free.
Thus, an additional filtering (e.g., Gaussian via :func:`pygmt.grdfilter`)
of the DiM-filtered data is generally recommended.
Full option list at :gmt-docs:`dimfilter.html`
{aliases}
Parameters
----------
{grid}
{outgrid}
distance : int or str
Distance flag tells how grid (x,y) relates to filter width, as follows:
- **0**\ : grid (x,y) in same units as *width*, Cartesian distances.
- **1**\ : grid (x,y) in degrees, *width* in kilometers, Cartesian
distances.
- **2**\ : grid (x,y) in degrees, *width* in km, dx scaled by
cos(middle y), Cartesian distances.
The above options are fastest because they allow weight matrix to be
computed only once. The next two options are slower because they
recompute weights for each latitude.
- **3**\ : grid (x,y) in degrees, *width* in km, dx scaled by
cosine(y), Cartesian distance calculation.
- **4**\ : grid (x,y) in degrees, *width* in km, Spherical distance
calculation.
filter : str
**x**\ *width*\ [**+l**\|\ **u**].
Set the primary filter type. Choose among convolution and
non-convolution filters. Use the filter code **x** followed by
the full diameter *width*. Available convolution filters are:
- (**b**) Boxcar: All weights are equal.
- (**c**) Cosine Arch: Weights follow a cosine arch curve.
- (**g**) Gaussian: Weights are given by the Gaussian function.
Non-convolution filters are:
- (**m**) Median: Returns median value.
- (**p**) Maximum likelihood probability (a mode estimator): Return
modal value. If more than one mode is found we return their average
value. Append **+l** or **+h** to the filter width if you want
to return the smallest or largest of each sector's modal values.
sectors : str
**x**\ *sectors*\ [**+l**\|\ **u**]
Set the secondary filter type **x** and the number of bow-tie sectors.
*sectors* must be integer and larger than 0. When *sectors* is
set to 1, the secondary filter is not effective. Available secondary
filters **x** are:
- (**l**) Lower: Return the minimum of all filtered values.
- (**u**) Upper: Return the maximum of all filtered values.
- (**a**) Average: Return the mean of all filtered values.
- (**m**) Median: Return the median of all filtered values.
- (**p**) Mode: Return the mode of all filtered values:
If more than one mode is found we return their average
value. Append **+l** or **+h** to the sectors if you rather want to
return the smallest or largest of the modal values.
spacing : str or list
*x_inc* [and optionally *y_inc*] is the output increment. Append
**m** to indicate minutes, or **c** to indicate seconds. If the new
*x_inc*, *y_inc* are NOT integer multiples of the old ones (in the
input data), filtering will be considerably slower. [Default is same
as input.]
region : str or list
[*xmin*, *xmax*, *ymin*, *ymax*].
Define the region of the output points [Default is same as input].
{verbose}
Returns
-------
ret: xarray.DataArray or None
Return type depends on whether the ``outgrid`` parameter is set:
- :class:`xarray.DataArray` if ``outgrid`` is not set
- None if ``outgrid`` is set (grid output will be stored in file set by
``outgrid``)
Example
-------
>>> import pygmt
>>> # Load a grid of Earth relief data
>>> grid = pygmt.datasets.load_earth_relief()
>>> # Create a filtered grid from an input grid.
>>> filtered_grid = pygmt.dimfilter(
... grid=grid,
... # Set filter type to "median" and the diameter width to 600 km
... filter="m600",
... # Set grid in degrees, width in km
... distance=4,
... # Create 6 sectors and return the lowest values in the sector
... sectors="l6",
... # Set the region longitude range from 55W to 51W, and the
... # latitude range from 24S to 19S
... region=[-55, -51, -24, -19],
... )
"""
if not all(arg in kwargs for arg in ["D", "F", "N"]) and "Q" not in kwargs:
raise GMTInvalidInput(
"""At least one of the following parameters must be specified:
distance, filters, or sectors."""
)
with GMTTempFile(suffix=".nc") as tmpfile:
with Session() as lib:
file_context = lib.virtualfile_from_data(check_kind="raster", data=grid)
with file_context as infile:
if (outgrid := kwargs.get("G")) is None:
kwargs["G"] = outgrid = tmpfile.name # output to tmpfile
lib.call_module(
module="dimfilter", args=build_arg_string(kwargs, infile=infile)
)
return load_dataarray(outgrid) if outgrid == tmpfile.name else None