From 2304053e4f297d71e6f59a7009743e16f1a78b2a Mon Sep 17 00:00:00 2001 From: Will Schlitzer Date: Mon, 11 Apr 2022 08:58:11 -0600 Subject: [PATCH] Wrap dimfilter (#1492) Co-authored-by: Meghan Jones Co-authored-by: Dongdong Tian Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com> --- doc/api/index.rst | 1 + pygmt/__init__.py | 1 + pygmt/src/__init__.py | 1 + pygmt/src/dimfilter.py | 144 ++++++++++++++++++++++++++++++++++ pygmt/tests/test_dimfilter.py | 81 +++++++++++++++++++ 5 files changed, 228 insertions(+) create mode 100644 pygmt/src/dimfilter.py create mode 100644 pygmt/tests/test_dimfilter.py diff --git a/doc/api/index.rst b/doc/api/index.rst index 5287a36836e..bce929e305f 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -135,6 +135,7 @@ Operations on raster data .. autosummary:: :toctree: generated + dimfilter grd2xyz grdclip grdcut diff --git a/pygmt/__init__.py b/pygmt/__init__.py index 9d0a7331743..ca1c804b11a 100644 --- a/pygmt/__init__.py +++ b/pygmt/__init__.py @@ -33,6 +33,7 @@ blockmedian, blockmode, config, + dimfilter, grd2cpt, grd2xyz, grdclip, diff --git a/pygmt/src/__init__.py b/pygmt/src/__init__.py index fc069a3dddb..1ba67a0f506 100644 --- a/pygmt/src/__init__.py +++ b/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 diff --git a/pygmt/src/dimfilter.py b/pygmt/src/dimfilter.py new file mode 100644 index 00000000000..72e7cd695fe --- /dev/null +++ b/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 diff --git a/pygmt/tests/test_dimfilter.py b/pygmt/tests/test_dimfilter.py new file mode 100644 index 00000000000..879a81d6b6c --- /dev/null +++ b/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)