Wrap sphinterpolate (#1418)

Wrapping the sphinterpolate function which creates spherical
grid files in tension of data. Official GMT documentation is at 
https://docs.generic-mapping-tools.org/6.2/sphdistance.html.
Aliased outgrid (G), spacing (I), region (R), verbose (V).

Co-authored-by: Wei Ji <[email protected]>
Co-authored-by: Dongdong Tian <[email protected]>

Author: 3 people

doc/api/index.rst

@@ -84,6 +84,7 @@ Operations on tabular data:
     blockmode
     nearneighbor
     sph2grd
+    sphinterpolate
     surface
 
 Operations on grids:

pygmt/__init__.py

@@ -51,6 +51,7 @@
     nearneighbor,
     sph2grd,
     sphdistance,
+    sphinterpolate,
     surface,
     which,
     x2sys_cross,

pygmt/src/__init__.py

@@ -39,6 +39,7 @@
 from pygmt.src.solar import solar
 from pygmt.src.sph2grd import sph2grd
 from pygmt.src.sphdistance import sphdistance
+from pygmt.src.sphinterpolate import sphinterpolate
 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/sphinterpolate.py

@@ -0,0 +1,69 @@
+"""
+sphinterpolate - Spherical gridding in tension of data on a sphere
+"""
+from pygmt.clib import Session
+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="spacing",
+    R="region",
+    V="verbose",
+)
+@kwargs_to_strings(I="sequence", R="sequence")
+def sphinterpolate(data, **kwargs):
+    r"""
+    Create spherical grid files in tension of data.
+
+    Reads a table containing *lon, lat, z* columns and performs a Delaunay
+    triangulation to set up a spherical interpolation in tension. Several
+    options may be used to affect the outcome, such as choosing local versus
+    global gradient estimation or optimize the tension selection to satisfy one
+    of four criteria.
+
+    Full option list at :gmt-docs:`sphinterpolate.html`
+
+    {aliases}
+
+    Parameters
+    ----------
+    data : str or {table-like}
+        Pass in (x, y, z) or (longitude, latitude, elevation) values by
+        providing a file name to an ASCII data table, a 2D
+        {table-classes}.
+    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``)
+    """
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        with Session() as lib:
+            file_context = lib.virtualfile_from_data(check_kind="vector", data=data)
+            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("sphinterpolate", arg_str)
+
+        return load_dataarray(outgrid) if outgrid == tmpfile.name else None

pygmt/tests/test_sphinterpolate.py

@@ -0,0 +1,42 @@
+"""
+Tests for sphinterpolate.
+"""
+import os
+
+import numpy.testing as npt
+import pytest
+from pygmt import sphinterpolate
+from pygmt.datasets import load_mars_shape
+from pygmt.helpers import GMTTempFile
+
+
+@pytest.fixture(scope="module", name="mars")
+def fixture_mars_shape():
+    """
+    Load the data from the sample bathymetry dataset.
+    """
+    return load_mars_shape()
+
+
+def test_sphinterpolate_outgrid(mars):
+    """
+    Test sphinterpolate with a set outgrid.
+    """
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        result = sphinterpolate(data=mars, outgrid=tmpfile.name, spacing=1, region="g")
+        assert result is None  # return value is None
+        assert os.path.exists(path=tmpfile.name)  # check that outgrid exists
+
+
+def test_sphinterpolate_no_outgrid(mars):
+    """
+    Test sphinterpolate with no set outgrid.
+    """
+    temp_grid = sphinterpolate(data=mars, spacing=1, region="g")
+    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(), 14628.144)
+    npt.assert_allclose(temp_grid.min(), -6908.1987)
+    npt.assert_allclose(temp_grid.median(), 118.96849)
+    npt.assert_allclose(temp_grid.mean(), 272.60593)