Wrap surface

Initial commit for GenericMappingTools#243. Implement GMT surface function under gridding.py. Test cases checking file/matrix/vectors type inputs stored in test_surface.py. Sample dataset for tests uses newly created load_tut_ship function in datasets/tutorial.py.

Note that original GMT surface module https://gmt.soest.hawaii.edu/doc/latest/surface.html requires a -Goutputfile.nc parameter. Here, the implementation of surface does not respect this -G parameter properly. Instead, it stores the output file to a GMTTempFile and returns an xarray.Dataset in Python.

gmt/__init__.py

@@ -13,6 +13,7 @@
 # Import modules to make the high-level GMT Python API
 from .session_management import begin as _begin, end as _end
 from .figure import Figure
+from .gridding import surface
 from .modules import info, grdinfo, which
 from . import datasets
 

gmt/datasets/__init__.py

@@ -1,5 +1,5 @@
 """
 Load sample data included with GMT (downloaded from the GMT cache server).
 """
-from .tutorial import load_japan_quakes, load_usgs_quakes
+from .tutorial import load_japan_quakes, load_tut_ship, load_usgs_quakes
 from .earth_relief import load_earth_relief

gmt/datasets/tutorial.py

@@ -38,6 +38,27 @@ def load_japan_quakes():
     return data
 
 
+def load_tut_ship():
+    """
+    Load a table of ship observations of bathymetry off Baja California as a
+    pandas.DataFrame.
+
+    This is the ``@tut_ship.xyz`` dataset used in the GMT tutorials.
+
+    The data are downloaded to a cache directory (usually ``~/.gmt/cache``) the
+    first time you invoke this function. Afterwards, it will load the data from
+    the cache. So you'll need an internet connection the first time around.
+
+    Returns
+    -------
+    data :  pandas.Dataframe
+        The data table. Columns are x, y, and z.
+    """
+    fname = which("@tut_ship.xyz", download="c")
+    data = pd.read_csv(fname, sep="\t", header=None, names=["x", "y", "z"])
+    return data
+
+
 def load_usgs_quakes():
     """
     Load a table of global earthquakes form the USGS as a pandas.Dataframe.

gmt/gridding.py

@@ -0,0 +1,71 @@
+"""
+GMT modules for Gridding of Data Tables
+"""
+import xarray as xr
+
+from .clib import Session
+from .helpers import (
+    build_arg_string,
+    fmt_docstring,
+    GMTTempFile,
+    use_alias,
+    data_kind,
+    dummy_context,
+)
+from .exceptions import GMTInvalidInput
+
+
+@fmt_docstring
+@use_alias(I="spacing", R="region")
+def surface(x=None, y=None, z=None, data=None, **kwargs):
+    """
+    Grids table data using adjustable tension continuous curvature splines
+
+    Takes a matrix, (x,y,z) pairs, or a file name as input.
+
+    Must provide either *data* or *x*, *y*, and *z*.
+
+    {gmt_module_docs}
+
+    Parameters
+    ----------
+    x, y, z : 1d arrays
+        Arrays of x and y coordinates and values z of the data points.
+    data : str or 2d array
+        Either a data file name or a 2d numpy array with the tabular data.
+
+    spacing (I) :
+        ``'xinc[unit][+e|n][/yinc[unit][+e|n]]'``.
+        x_inc [and optionally y_inc] is the grid spacing.
+
+    region (R) : str or list
+        ``'xmin/xmax/ymin/ymax[+r][+uunit]'``.
+        Specify the region of interest.
+
+    {aliases}
+
+    Returns
+    -------
+    array: xarray.DataArray
+        The output grid as a DataArray
+    """
+    kind = data_kind(data, x, y, z)
+    if kind == "vectors" and z is None:
+        raise GMTInvalidInput("Must provide z with x and y.")
+
+    with GMTTempFile(suffix=".nc") as outfile:
+        with Session() as lib:
+            if kind == "file":
+                file_context = dummy_context(data)
+            elif kind == "matrix":
+                file_context = lib.virtualfile_from_matrix(data)
+            elif kind == "vectors":
+                file_context = lib.virtualfile_from_vectors(x, y, z)
+            else:
+                raise GMTInvalidInput("Unrecognized data type: {}".format(type(data)))
+            with file_context as infile:
+                kwargs.update({"G": outfile.name})
+                arg_str = " ".join([infile, build_arg_string(kwargs)])
+                lib.call_module(module="surface", args=arg_str)
+        result = xr.open_dataset(outfile)
+    return result

gmt/tests/test_surface.py

@@ -0,0 +1,36 @@
+"""
+Tests for surface
+"""
+from .. import surface
+from .. import which
+from ..datasets import load_tut_ship
+
+
+def test_surface_input_file():
+    """
+    Run surface by passing in a filename
+    """
+    fname = which("@tut_ship.xyz", download="c")
+    outputfile = surface(data=fname, I="5m", R="245/255/20/30")
+    return outputfile
+
+
+def test_surface_input_data_array():
+    """
+    Run surface by passing in a numpy array into data
+    """
+    ship_data = load_tut_ship()
+    data = ship_data.values  # convert pandas.DataFrame to numpy.ndarray
+    outputfile = surface(data=data, I="5m", R="245/255/20/30")
+    return outputfile
+
+
+def test_surface_input_xyz():
+    """
+    Run surface by passing in x, y, z numpy.ndarrays individually
+    """
+    ship_data = load_tut_ship()
+    outputfile = surface(
+        x=ship_data.x, y=ship_data.y, z=ship_data.z, I="5m", R="245/255/20/30"
+    )
+    return outputfile