Wrap grdhisteq

doc/api/index.rst

@@ -92,6 +92,7 @@ Operations on grids:
     grdcut
     grdfill
     grdfilter
+    grdhisteq
     grdlandmask
     grdgradient
     grdsample

pygmt/__init__.py

@@ -38,6 +38,7 @@
     grdfill,
     grdfilter,
     grdgradient,
+    grdhisteq,
     grdinfo,
     grdlandmask,
     grdsample,

pygmt/src/__init__.py

@@ -16,6 +16,7 @@
 from pygmt.src.grdfill import grdfill
 from pygmt.src.grdfilter import grdfilter
 from pygmt.src.grdgradient import grdgradient
+from pygmt.src.grdhisteq import grdhisteq
 from pygmt.src.grdimage import grdimage
 from pygmt.src.grdinfo import grdinfo
 from pygmt.src.grdlandmask import grdlandmask

pygmt/src/grdhisteq.py

@@ -0,0 +1,74 @@
+"""
+grdhisteq - Compute directional gradients from a grid.
+"""
+
+import xarray as xr
+from pygmt.clib import Session
+from pygmt.helpers import (
+    GMTTempFile,
+    build_arg_string,
+    fmt_docstring,
+    kwargs_to_strings,
+    use_alias,
+)
+
+
+@fmt_docstring
+@use_alias(
+    C="divisions",
+    G="outgrid",
+    R="region",
+    N="gaussian",
+    Q="quadratic",
+    V="verbose",
+    n="interpolation",
+)
+@kwargs_to_strings(R="sequence")
+def grdhisteq(grid, **kwargs):
+    r"""
+    Perform histogram equalization for a grid.
+
+    Full option list at :gmt-docs:`grdhisteq.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.
+    divisions : int or str
+        The number of divisions of data range to make [Default is 16].
+
+    {R}
+    {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``)
+    """
+    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 "G" not in kwargs.keys():  # 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("grdhisteq", arg_str)
+
+        if outgrid == tmpfile.name:  # if user did not set outgrid, return DataArray
+            with xr.open_dataarray(outgrid) as dataarray:
+                result = dataarray.load()
+                _ = result.gmt  # load GMTDataArray accessor information
+        else:
+            result = None  # if user sets an outgrid, return None
+
+        return result

pygmt/tests/test_grdhisteq.py

@@ -0,0 +1,42 @@
+"""
+Tests for grdhisteq.
+"""
+import os
+
+import numpy.testing as npt
+import pytest
+from pygmt import grdhisteq, grdinfo
+from pygmt.datasets import load_earth_relief
+from pygmt.helpers import GMTTempFile
+
+
+@pytest.fixture(scope="module", name="grid")
+def fixture_grid():
+    """
+    Load the grid data from the sample earth_relief file.
+    """
+    return load_earth_relief(resolution="01d", region=[-5, 5, -5, 5])
+
+
+def test_grdhisteq_outgrid(grid):
+    """
+    Test the azimuth and direction parameters for grdhisteq with a set outgrid.
+    """
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        result = grdhisteq(grid=grid, outgrid=tmpfile.name)
+        assert result is None  # return value is None
+        assert os.path.exists(path=tmpfile.name)  # check that outgrid exists
+        result = (
+            grdinfo(grid=tmpfile.name, force_scan="a", per_column="n").strip().split()
+        )
+
+
+def test_grdhisteq_no_outgrid(grid):
+    """
+    Test the azimuth and direction parameters for grdhisteq with no set
+    outgrid.
+    """
+    temp_grid = grdhisteq(grid=grid)
+    assert temp_grid.dims == ("lat", "lon")
+    assert temp_grid.gmt.gtype == 1  # Geographic grid
+    assert temp_grid.gmt.registration == 1  # Pixel registration