Wrap triangulate

doc/api/index.rst

@@ -82,6 +82,7 @@ Operations on tabular data:
     blockmean
     blockmedian
     surface
+    triangulate
 
 Operations on grids:
 

pygmt/__init__.py

@@ -47,6 +47,7 @@
     info,
     makecpt,
     surface,
+    triangulate,
     which,
     x2sys_cross,
     x2sys_init,

pygmt/src/__init__.py

@@ -38,6 +38,7 @@
 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_"
+from pygmt.src.triangulate import triangulate
 from pygmt.src.velo import velo
 from pygmt.src.which import which
 from pygmt.src.wiggle import wiggle

pygmt/src/triangulate.py

@@ -0,0 +1,116 @@
+"""
+triangulate - Delaunay triangulation or Voronoi partitioning and gridding of
+Cartesian data.
+"""
+import pandas as pd
+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",
+    J="projection",
+    R="region",
+    V="verbose",
+    b="binary",
+    d="nodata",
+    e="find",
+    f="coltypes",
+    h="header",
+    i="incols",
+    r="registration",
+)
+@kwargs_to_strings(R="sequence")
+def triangulate(table=None, x=None, y=None, z=None, **kwargs):
+    """
+    Delaunay triangulation or Voronoi partitioning and gridding of Cartesian
+    data.
+
+    Triangulate reads in x,y[,z] data and performs Delaunay triangulation,
+    i.e., it find how the points should be connected to give the most
+    equilateral triangulation possible. If a map projection (give *region*
+    and *projection*) is chosen then it is applied before the triangulation
+    is calculated.
+
+    Must provide either *table* or *x*, *y*, and *z*.
+
+    Full option list at :gmt-docs:`triangulate.html`
+
+    {aliases}
+
+    Parameters
+    ----------
+    x/y/z : np.ndarray
+        Arrays of x and y coordinates and values z of the data points.
+    table : 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}.
+    projection : str
+        Select map projection.
+    region
+        ``'xmin/xmax/ymin/ymax[+r][+uunit]'``.
+        Specify the region of interest.
+    spacing : str
+        ``'xinc[unit][+e|n][/yinc[unit][+e|n]]'``.
+        x_inc [and optionally y_inc] is the grid spacing.
+    outgrid : bool or str
+        Use triangulation to grid the data onto an even grid (specified with
+        *region* and *spacing*). Set to True, or pass in the name of the output
+        grid file. The interpolation is performed in the original coordinates,
+        so if your triangles are close to the poles you are better off
+        projecting all data to a local coordinate system before using
+        *triangulate* (this is true of all gridding routines) or instead
+        select *sphtriangulate*.
+    {V}
+    {b}
+    {d}
+    {e}
+    {f}
+    {h}
+    {i}
+    {r}
+        Only valid with *outgrid*.
+
+    Returns
+    -------
+    ret: xarray.DataArray or None
+        Return type depends on whether the outgrid parameter is set:
+
+        - pandas.DataFrame if outgrid is None (default)
+        - xarray.DataArray if outgrid is True
+        - None if outgrid is a str (grid output will be stored in outgrid)
+    """
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        with Session() as lib:
+            # Choose how data will be passed into the module
+            table_context = lib.virtualfile_from_data(
+                check_kind="vector", data=table, x=x, y=y, z=z
+            )
+            with table_context as infile:
+                if "G" not in kwargs:  # table output if outgrid is unset
+                    kwargs.update({">": tmpfile.name})
+                else:  # NetCDF or xarray.DataArray output if outgrid is set
+                    if (
+                        kwargs["G"] is True
+                    ):  # xarray.DataArray output if outgrid is True
+                        kwargs.update({"G": tmpfile.name})
+                    outgrid = kwargs["G"]
+                arg_str = " ".join([infile, build_arg_string(kwargs)])
+                lib.call_module(module="triangulate", args=arg_str)
+
+        try:
+            result = load_dataarray(outgrid) if outgrid == tmpfile.name else None
+        except UnboundLocalError:  # if outgrid unset, return pd.DataFrame
+            result = pd.read_csv(tmpfile.name, sep="\t", header=None)
+
+    return result

pygmt/tests/test_triangulate.py

@@ -0,0 +1,99 @@
+"""
+Tests for triangulate.
+"""
+import os
+
+import pandas as pd
+import pytest
+import xarray as xr
+from pygmt import triangulate
+from pygmt.datasets import load_sample_bathymetry
+from pygmt.exceptions import GMTInvalidInput
+from pygmt.helpers import GMTTempFile, data_kind
+
+
+@pytest.fixture(scope="module", name="dataframe")
+def fixture_dataframe():
+    """
+    Load the table data from the sample bathymetry dataset.
+    """
+    return load_sample_bathymetry()
+
+
+def test_triangulate_input_file():
+    """
+    Run triangulate by passing in a filename.
+    """
+    output = triangulate(table="@tut_ship.xyz")
+    assert isinstance(output, pd.DataFrame)
+    assert output.shape == (161935, 3)
+
+
+def test_triangulate_input_data_array(dataframe):
+    """
+    Run triangulate by passing in a numpy array into data.
+    """
+    data = dataframe.to_numpy()
+    output = triangulate(table=data)
+    assert isinstance(output, pd.DataFrame)
+    assert output.shape == (161935, 3)
+
+
+def test_triangulate_input_xyz(dataframe):
+    """
+    Run triangulate by passing in x, y, z numpy.ndarrays individually.
+    """
+    output = triangulate(
+        x=dataframe.longitude,
+        y=dataframe.latitude,
+        z=dataframe.bathymetry,
+    )
+    assert isinstance(output, pd.DataFrame)
+    assert output.shape == (161935, 3)
+
+
+def test_triangulate_input_xy_no_z(dataframe):
+    """
+    Run triangulate by passing in x and y, but no z.
+    """
+    output = triangulate(x=dataframe.longitude, y=dataframe.latitude)
+    assert isinstance(output, pd.DataFrame)
+    assert output.shape == (161935, 3)
+
+
+def test_triangulate_wrong_kind_of_input(dataframe):
+    """
+    Run triangulate using grid input that is not file/matrix/vectors.
+    """
+    data = dataframe.bathymetry.to_xarray()  # convert pandas.Series to xarray.DataArray
+    assert data_kind(data) == "grid"
+    with pytest.raises(GMTInvalidInput):
+        triangulate(table=data)
+
+
+def test_triangulate_with_outgrid_true(dataframe):
+    """
+    Run triangulate with outgrid=True and see it load into an xarray.DataArray.
+    """
+    data = dataframe.to_numpy()
+    output = triangulate(
+        table=data, spacing="5m", region=[245, 255, 20, 30], outgrid=True
+    )
+    assert isinstance(output, xr.DataArray)
+    assert output.shape == (121, 121)
+
+
+def test_triangulate_with_outgrid_param(dataframe):
+    """
+    Run triangulate with the -Goutputfile.nc parameter.
+    """
+    data = dataframe.to_numpy()
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        output = triangulate(
+            table=data, spacing="5m", region=[245, 255, 20, 30], outgrid=tmpfile.name
+        )
+        assert output is None  # check that output is None since outgrid is set
+        assert os.path.exists(path=tmpfile.name)  # check that outgrid exists
+        with xr.open_dataarray(tmpfile.name) as grid:
+            assert isinstance(grid, xr.DataArray)  # ensure netcdf grid loads ok
+            assert grid.shape == (121, 121)