Add implementation of nearneighbor

The implementation mirrors that of the already implemented surface.

The documentation has been cobbled together from similar portions in
surface and from the [official GMT documentation on
nearneighbor](https://docs.generic-mapping-tools.org/latest/nearneighbor.html).
Additional long-form flags have been added, reflecting the flags
particularly useful for nearneighbor, mainly `-Ssearch_radius`,
`-Eempty` and `-Nsectors`.

Tests have been blatently copied from those testing surface, since the
two functions should operate broadly similarly.

pygmt/__init__.py

@@ -43,6 +43,7 @@
  grdtrack,
  info,
  makecpt,
+ nearneighbor,
  surface,
  which,
  x2sys_cross,

pygmt/src/__init__.py

@@ -29,6 +29,7 @@
 from pygmt.src.logo import logo
 from pygmt.src.makecpt import makecpt
 from pygmt.src.meca import meca
+from pygmt.src.nearneighbor import nearneighbor
 from pygmt.src.plot import plot
 from pygmt.src.plot3d import plot3d
 from pygmt.src.rose import rose

pygmt/src/nearneighbor.py

@@ -0,0 +1,131 @@
+"""
+nearneighbor - Grid table data using a "Nearest neighbor" algorithm
+"""
+
+import xarray as xr
+from pygmt.clib import Session
+from pygmt.exceptions import GMTInvalidInput
+from pygmt.helpers import (
+ GMTTempFile,
+ build_arg_string,
+ data_kind,
+ dummy_context,
+ fmt_docstring,
+ kwargs_to_strings,
+ use_alias,
+)
+
+
+@fmt_docstring
+@use_alias(
+ I="spacing",
+ R="region",
+ V="verbose",
+ G="outfile",
+ a="aspatial",
+ f="coltypes",
+ r="registration",
+ S="search_radius",
+ E="empty",
+ N="sectors",
+)
+@kwargs_to_strings(R="sequence")
+def nearneighbor(x=None, y=None, z=None, data=None, **kwargs):
+ r"""
+ Grid table data using a "Nearest neighbor" algorithm
+
+ **nearneighbor** reads arbitrarily located (*x,y,z*\ [,\ *w*]) triples
+ [quadruplets] and uses a nearest neighbor algorithm to assign a weighted
+ average value to each node that has one or more data points within a search
+ radius centered on the node with adequate coverage across a subset of the
+ chosen sectors. The node value is computed as a weighted mean of the
+ nearest point from each sector inside the search radius. 
+
+ Takes a matrix, xyz triples, or a file name as input.
+
+ Must provide either ``data`` or ``x``, ``y``, and ``z``.
+
+ Full option list at :gmt-docs:`nearneighbor.html`
+
+ {aliases}
+
+ 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.
+
+ {I}
+
+ region : str or list
+ *xmin/xmax/ymin/ymax*\[**+r**][**+u**\ *unit*].
+ Specify the region of interest.
+
+ search_radius : str
+ Sets the search radius that determines which data points are considered close to a node.
+
+ outfile : str
+ Optional. The file name for the output netcdf file with extension .nc
+ to store the grid in.
+
+ {V}
+ {a}
+ {f}
+ {r}
+
+ empty : str
+ Optional. Set the value assigned to empty nodes. Defaults to NaN.
+
+ sectors : str
+ Optional. **\ *sectors*\ [**+m**\ *min_sectors*]\|\ **n**
+ The circular search area centered on each node is divided into sectors
+ sectors. Average values will only be computed if there is at least one
+ value inside each of at least min_sectors of the sectors for a given
+ node. Nodes that fail this test are assigned the value NaN (but see
+ -E). If +m is omitted then min_sectors is set to be at least 50% of
+ sectors (i.e., rounded up to next integer) [Default is a quadrant
+ search with 100% coverage, i.e., sectors = min_sectors = 4]. Note that
+ only the nearest value per sector enters into the averaging; the more
+ distant points are ignored. Alternatively, use -Nn to call GDALʻs
+ nearest neighbor algorithm instead.
+
+ Returns
+ -------
+ ret: xarray.DataArray or None
+ Return type depends on whether the ``outfile`` parameter is set:
+
+ - :class:`xarray.DataArray`: if ``outfile`` is not set
+ - None if ``outfile`` is set (grid output will be stored in file set by
+ ``outfile``)
+ """
+
+ 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 tmpfile:
+ 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:
+ if "G" not in kwargs.keys(): # if outfile is unset, output to tmpfile
+ kwargs.update({"G": tmpfile.name})
+ outfile = kwargs["G"]
+ arg_str = " ".join([infile, build_arg_string(kwargs)])
+ lib.call_module(module="nearneighbor", args=arg_str)
+
+ if outfile == tmpfile.name: # if user did not set outfile, return DataArray
+ with xr.open_dataarray(outfile) as dataarray:
+ result = dataarray.load()
+ _ = result.gmt # load GMTDataArray accessor information
+ elif outfile != tmpfile.name: # if user sets an outfile, return None
+ result = None
+
+ return result

pygmt/tests/test_nearneighbor.py

@@ -0,0 +1,119 @@
+"""
+Tests for nearneighbor.
+"""
+import os
+
+import pytest
+import xarray as xr
+from pygmt import nearneighbor, which
+from pygmt.datasets import load_sample_bathymetry
+from pygmt.exceptions import GMTInvalidInput
+from pygmt.helpers import data_kind
+
+TEST_DATA_DIR = os.path.join(os.path.dirname(__file__), "data")
+TEMP_GRID = os.path.join(TEST_DATA_DIR, "tmp_grid.nc")
+
+
+@pytest.fixture(scope="module", name="ship_data")
+def fixture_ship_data():
+ """
+ Load the data from the sample bathymetry dataset.
+ """
+ return load_sample_bathymetry()
+
+
+def test_nearneighbor_input_file():
+ """
+ Run nearneighbor by passing in a filename.
+ """
+ fname = which("@tut_ship.xyz", download="c")
+ output = nearneighbor(data=fname, spacing="5m", region=[245, 255, 20, 30], search_radius="10m")
+ assert isinstance(output, xr.DataArray)
+ assert output.gmt.registration == 0 # Gridline registration
+ assert output.gmt.gtype == 0 # Cartesian type
+ return output
+
+
+def test_nearneighbor_input_data_array(ship_data):
+ """
+ Run nearneighbor by passing in a numpy array into data.
+ """
+ data = ship_data.values # convert pandas.DataFrame to numpy.ndarray
+ output = nearneighbor(data=data, spacing="5m", region=[245, 255, 20, 30], search_radius="10m")
+ assert isinstance(output, xr.DataArray)
+ return output
+
+
+def test_nearneighbor_input_xyz(ship_data):
+ """
+ Run nearneighbor by passing in x, y, z numpy.ndarrays individually.
+ """
+ output = nearneighbor(
+ x=ship_data.longitude,
+ y=ship_data.latitude,
+ z=ship_data.bathymetry,
+ spacing="5m",
+ region=[245, 255, 20, 30],
+ search_radius="10m",
+ )
+ assert isinstance(output, xr.DataArray)
+ return output
+
+
+def test_nearneighbor_input_xy_no_z(ship_data):
+ """
+ Run nearneighbor by passing in x and y, but no z.
+ """
+ with pytest.raises(GMTInvalidInput):
+ nearneighbor(
+ x=ship_data.longitude,
+ y=ship_data.latitude,
+ spacing="5m",
+ region=[245, 255, 20, 30],
+ search_radius="10m",
+ )
+
+
+def test_nearneighbor_wrong_kind_of_input(ship_data):
+ """
+ Run nearneighbor using grid input that is not file/matrix/vectors.
+ """
+ data = ship_data.bathymetry.to_xarray() # convert pandas.Series to xarray.DataArray
+ assert data_kind(data) == "grid"
+ with pytest.raises(GMTInvalidInput):
+ nearneighbor(data=data, spacing="5m", region=[245, 255, 20, 30], search_radius="10m")
+
+
+def test_nearneighbor_with_outfile_param(ship_data):
+ """
+ Run nearneighbor with the -Goutputfile.nc parameter.
+ """
+ data = ship_data.values # convert pandas.DataFrame to numpy.ndarray
+ try:
+ output = nearneighbor(
+ data=data, spacing="5m", region=[245, 255, 20, 30], outfile=TEMP_GRID, search_radius="10m"
+ )
+ assert output is None # check that output is None since outfile is set
+ assert os.path.exists(path=TEMP_GRID) # check that outfile exists at path
+ with xr.open_dataarray(TEMP_GRID) as grid:
+ assert isinstance(grid, xr.DataArray) # ensure netcdf grid loads ok
+ finally:
+ os.remove(path=TEMP_GRID)
+ return output
+
+
+def test_nearneighbor_short_aliases(ship_data):
+ """
+ Run nearneighbor using short aliases -I for spacing, -R for region, -G for
+ outfile, -S for search radius.
+ """
+ data = ship_data.values # convert pandas.DataFrame to numpy.ndarray
+ try:
+ output = nearneighbor(data=data, I="5m", R=[245, 255, 20, 30], G=TEMP_GRID, S="10m")
+ assert output is None # check that output is None since outfile is set
+ assert os.path.exists(path=TEMP_GRID) # check that outfile exists at path
+ with xr.open_dataarray(TEMP_GRID) as grid:
+ assert isinstance(grid, xr.DataArray) # ensure netcdf grid loads ok
+ finally:
+ os.remove(path=TEMP_GRID)
+ return output