Wrap sphdistance

doc/api/index.rst

@@ -97,6 +97,7 @@ Operations on grids:
     grdproject
     grdsample
     grdtrack
+    sphdistance
     xyz2grd
 
 Crossover analysis with x2sys:

pygmt/__init__.py

@@ -46,6 +46,7 @@
     grdtrack,
     info,
     makecpt,
+    sphdistance,
     surface,
     which,
     x2sys_cross,

pygmt/src/__init__.py

@@ -35,6 +35,7 @@
 from pygmt.src.plot3d import plot3d
 from pygmt.src.rose import rose
 from pygmt.src.solar import solar
+from pygmt.src.sphdistance import sphdistance
 from pygmt.src.subplot import set_panel, subplot
 from pygmt.src.surface import surface
 from pygmt.src.text import text_ as text  # "text" is an argument within "text_"

pygmt/src/sphdistance.py

@@ -0,0 +1,66 @@
+"""
+sphdistance - Create Voronoi distance, node,
+or natural nearest-neighbor grid on a sphere
+"""
+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(
+    G="outgrid",
+    I="increment",
+    R="region",
+    V="verbose",
+)
+@kwargs_to_strings(I="sequence", R="sequence")
+def sphdistance(table, **kwargs):
+    r"""
+    Create Voroni polygons from lat/lon coordinates.
+
+    Reads one or more ASCII [or binary] files (or standard
+    input) containing lon, lat and performs the construction of Voronoi
+    polygons. These polygons are then processed to calculate the nearest
+    distance to each node of the lattice and written to the specified grid.
+
+    {aliases}
+
+    Parameters
+    ----------
+    outgrid : str or None
+        The name of the output netCDF file with extension .nc to store the grid
+        in.
+    {I}
+    {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``)
+    """
+    if "I" not in kwargs.keys() or "R" not in kwargs.keys():
+        raise GMTInvalidInput("Both 'region' and 'increment' must be specified.")
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        with Session() as lib:
+            file_context = lib.virtualfile_from_data(check_kind="vector", data=table)
+            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 = build_arg_string(kwargs)
+                arg_str = " ".join([infile, arg_str])
+                lib.call_module("sphdistance", arg_str)
+
+        return load_dataarray(outgrid) if outgrid == tmpfile.name else None

pygmt/tests/test_sphdistance.py

@@ -0,0 +1,55 @@
+"""
+Tests for sphdistance.
+"""
+import os
+
+import numpy as np
+import numpy.testing as npt
+import pytest
+from pygmt import sphdistance
+from pygmt.exceptions import GMTInvalidInput
+from pygmt.helpers import GMTTempFile
+
+
+@pytest.fixture(scope="module", name="array")
+def fixture_table():
+    """
+    Load the table data.
+    """
+    coords_list = [[85.5, 22.3], [82.3, 22.6], [85.8, 22.4], [86.5, 23.3]]
+    return np.array(coords_list)
+
+
+def test_sphdistance_outgrid(array):
+    """
+    Test sphdistance with a set outgrid.
+    """
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        result = sphdistance(
+            table=array, outgrid=tmpfile.name, increment=1, region=[82, 87, 22, 24]
+        )
+        assert result is None  # return value is None
+        assert os.path.exists(path=tmpfile.name)  # check that outgrid exists
+
+
+def test_sphdistance_no_outgrid(array):
+    """
+    Test sphdistance with no set outgrid.
+    """
+    temp_grid = sphdistance(table=array, increment=[1, 2], region=[82, 87, 22, 24])
+    assert temp_grid.dims == ("lat", "lon")
+    assert temp_grid.gmt.gtype == 1  # Geographic grid
+    assert temp_grid.gmt.registration == 0  # Gridline registration
+    npt.assert_allclose(temp_grid.max(), 232977.546875)
+    npt.assert_allclose(temp_grid.min(), 0)
+    npt.assert_allclose(temp_grid.median(), 0)
+    npt.assert_allclose(temp_grid.mean(), 62469.17)
+
+
+def test_sphdistance_fails(array):
+    """
+    Check that sphdistance fails correctly when neither increment nor region is
+    given.
+    """
+    with pytest.raises(GMTInvalidInput):
+        sphdistance(table=array)