Wrap dimfilter

doc/api/index.rst

@@ -135,6 +135,7 @@ Operations on raster data
 .. autosummary::
     :toctree: generated
 
+    dimfilter
     grd2xyz
     grdclip
     grdcut

pygmt/__init__.py

@@ -33,6 +33,7 @@
     blockmedian,
     blockmode,
     config,
+    dimfilter,
     grd2cpt,
     grd2xyz,
     grdclip,

pygmt/src/__init__.py

@@ -9,6 +9,7 @@
 from pygmt.src.colorbar import colorbar
 from pygmt.src.config import config
 from pygmt.src.contour import contour
+from pygmt.src.dimfilter import dimfilter
 from pygmt.src.grd2cpt import grd2cpt
 from pygmt.src.grd2xyz import grd2xyz
 from pygmt.src.grdclip import grdclip

pygmt/src/dimfilter.py

@@ -0,0 +1,144 @@
+"""
+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
+
+
+@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. **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 is 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 : str or xarray.DataArray
+        The file name of the input grid or the grid loaded as a DataArray.
+    outgrid : str or None
+        The name of the output netCDF file with extension .nc to store the grid
+        in.
+    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**].
+        Sets 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**]
+        Sets 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: Same
+        as input.]
+    region : str or list
+        [*xmin*, *xmax*, *ymin*, *ymax*].
+        Defines the region of the output points. [Default: Same as input.]
+    {V}
+
+    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``)
+    """
+    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=None, data=grid)
+            with file_context as infile:
+                if "G" not in kwargs:  # if outgrid is unset, output to tempfile
+                    kwargs.update({"G": tmpfile.name})
+                outgrid = kwargs["G"]
+                arg_str = " ".join([infile, build_arg_string(kwargs)])
+                lib.call_module("dimfilter", arg_str)
+
+        return load_dataarray(outgrid) if outgrid == tmpfile.name else None

pygmt/tests/test_dimfilter.py

@@ -0,0 +1,81 @@
+"""
+Tests for dimfilter.
+"""
+import os
+
+import pytest
+import xarray as xr
+from pygmt import dimfilter, load_dataarray
+from pygmt.exceptions import GMTInvalidInput
+from pygmt.helpers import GMTTempFile
+from pygmt.helpers.testing import load_static_earth_relief
+
+
+@pytest.fixture(scope="module", name="grid")
+def fixture_grid():
+    """
+    Load the grid data from the static_earth_relief file.
+    """
+    return load_static_earth_relief()
+
+
+@pytest.fixture(scope="module", name="expected_grid")
+def fixture_grid_result():
+    """
+    Load the expected dimfilter grid result.
+    """
+    return xr.DataArray(
+        data=[
+            [346.0, 344.5, 349.0, 349.0],
+            [344.5, 318.5, 344.5, 394.0],
+            [344.5, 356.5, 345.5, 352.5],
+            [367.5, 349.0, 385.5, 349.0],
+            [435.0, 385.5, 413.5, 481.5],
+        ],
+        coords=dict(
+            lon=[-54.5, -53.5, -52.5, -51.5],
+            lat=[-23.5, -22.5, -21.5, -20.5, -19.5],
+        ),
+        dims=["lat", "lon"],
+    )
+
+
+def test_dimfilter_outgrid(grid, expected_grid):
+    """
+    Test the required parameters for dimfilter with a set outgrid.
+    """
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        result = dimfilter(
+            grid=grid,
+            outgrid=tmpfile.name,
+            filter="m600",
+            distance=4,
+            sectors="l6",
+            region=[-55, -51, -24, -19],
+        )
+        assert result is None  # return value is None
+        assert os.path.exists(path=tmpfile.name)  # check that outgrid exists
+        temp_grid = load_dataarray(tmpfile.name)
+        xr.testing.assert_allclose(a=temp_grid, b=expected_grid)
+
+
+def test_dimfilter_no_outgrid(grid, expected_grid):
+    """
+    Test the required parameters for dimfilter with no set outgrid.
+    """
+    result = dimfilter(
+        grid=grid, filter="m600", distance=4, sectors="l6", region=[-55, -51, -24, -19]
+    )
+    assert result.dims == ("lat", "lon")
+    assert result.gmt.gtype == 1  # Geographic grid
+    assert result.gmt.registration == 1  # Pixel registration
+    xr.testing.assert_allclose(a=result, b=expected_grid)
+
+
+def test_dimfilter_fails(grid):
+    """
+    Check that dimfilter fails correctly when not all of sectors, filters, and
+    distance are specified.
+    """
+    with pytest.raises(GMTInvalidInput):
+        dimfilter(grid=grid, sectors="l6", distance=4)