Wrap sph2grd

doc/api/index.rst

@@ -96,6 +96,7 @@ Operations on grids:
     grdgradient
     grdsample
     grdtrack
+    sph2grd
     xyz2grd
 
 Crossover analysis with x2sys:

pygmt/__init__.py

@@ -44,6 +44,7 @@
     grdtrack,
     info,
     makecpt,
+    sph2grd,
     surface,
     which,
     x2sys_cross,

pygmt/src/__init__.py

@@ -34,6 +34,7 @@
 from pygmt.src.plot3d import plot3d
 from pygmt.src.rose import rose
 from pygmt.src.solar import solar
+from pygmt.src.sph2grd import sph2grd
 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/sph2grd.py

@@ -0,0 +1,68 @@
+"""
+sph2grd - Compute grid from spherical harmonic coefficients
+"""
+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(
+    G="outgrid",
+    I="increment",
+    R="region",
+    V="verbose",
+)
+@kwargs_to_strings(I="sequence", R="sequence")
+def sph2grd(table, **kwargs):
+    r"""
+    Create spherical grid files in tension of data.
+
+    Reads a spherical harmonics coefficient table with records of L, M,
+    C[L,M], S[L,M] and evaluates the spherical harmonic model on 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``)
+    """
+    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("sph2grd", 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_sph2grd.py

@@ -0,0 +1,34 @@
+"""
+Tests for sph2grd.
+"""
+import os
+
+import numpy.testing as npt
+from pygmt import sph2grd
+from pygmt.helpers import GMTTempFile
+
+
+def test_sph2grd_outgrid():
+    """
+    Test sph2grd with a set outgrid.
+    """
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        result = sph2grd(
+            table="@EGM96_to_36.txt", outgrid=tmpfile.name, increment=1, region="g"
+        )
+        assert result is None  # return value is None
+        assert os.path.exists(path=tmpfile.name)  # check that outgrid exists
+
+
+def test_sph2grd_no_outgrid():
+    """
+    Test sph2grd with no set outgrid.
+    """
+    temp_grid = sph2grd(table="@EGM96_to_36.txt", increment=1, region="g")
+    assert temp_grid.dims == ("y", "x")
+    assert temp_grid.gmt.gtype == 0  # Cartesian grid
+    assert temp_grid.gmt.registration == 0  # Gridline registration
+    npt.assert_allclose(temp_grid.max(), 0.00021961, rtol=1e-4)
+    npt.assert_allclose(temp_grid.min(), -0.0004326, rtol=1e-4)
+    npt.assert_allclose(temp_grid.median(), -0.00010894, rtol=1e-4)
+    npt.assert_allclose(temp_grid.mean(), -0.00010968, rtol=1e-4)