Wrap sphinterpolate

doc/api/index.rst

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

pygmt/__init__.py

@@ -50,6 +50,7 @@
     makecpt,
     nearneighbor,
     sphdistance,
+    sphinterpolate,
     surface,
     which,
     x2sys_cross,

pygmt/src/__init__.py

@@ -38,6 +38,7 @@
 from pygmt.src.rose import rose
 from pygmt.src.solar import solar
 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,68 @@
+"""
+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 one or more ASCII [or binary] files (or standard input) containing
+    *lon, lat, z* and performs a Delaunay triangulation to set up a spherical
+    interpolation in tension. The final grid is saved to the specified file.
+    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.
+
+    {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)