pygmt.sphdistance: Create Voronoi distance, node, or natural nearest-neighbor grid on a sphere

[seisman]

This issue serves as the central place for discussing and tracking the implementation of the pygmt.sphdistance function in PyGMT. The issue will be closed when the initial implementation is complete. Progress is tracked at PyGMT: Wrapping GMT modules.

Documentation

• GMT: https://docs.generic-mapping-tools.org/dev/sphdistance.html
• GMT.jl: https://www.generic-mapping-tools.org/GMTjl_doc/documentation/modules/sphdistance
• PyGMT: https://www.pygmt.org/dev/api/generated/pygmt.sphdistance.html

GMT Option Flags and Modifiers

☑️: Implemented; ⬜: To be implemented/discussed; Strikethrough: Won't implement.

• ☑️ -C (single_form): Store only one form of location coordinates at a time to save memory (geographic or Cartesian 3-D vectors).
• ☑️ -D (duplicate): Skip duplicate points.
• ☑️ -E (quantity): d|n|z[dist]. Specify the quantity assigned to grid nodes: distance (d), Voronoi polygon ID (n), or natural nearest-neighbor z-value (z).
• ☑️ -G (outgrid): Output grid file name.
• ☑️ -I (spacing): Grid spacing.
• ☑️ -L (unit): Unit used for distance calculations.
• ☑️ -N (node_table): Read Voronoi polygon node information from a separate file.
• ☑️ -Q (voronoi): Use pre-calculated Voronoi polygons from a file instead of computing them.
• ☑️ -R (region): Output grid region.
• ☑️ -V (verbose): Verbosity level.
• -X/-Y: Use Figure.shift_origin instead.
• ⬜ -b: Binary input/output.
• ⬜ -d: Replace NaN with a specified nodata value on input/output.
• ⬜ -e: Pattern matching to select input rows.
• ⬜ -f (coltypes): Column data types.
• ⬜ -g: Gap detection.
• ⬜ -h: Read/write header records.
• ⬜ -i: Select input columns.
• ⬜ -j: Distance calculation mode (great circle, rhumb, Cartesian).
• ⬜ -q: Select rows by row number or range.
• ⬜ -r: Set grid node registration to gridline or pixel.
• ⬜ -s: Skip rows containing NaN values.
• ⬜ -w: Wrap repeated cycles.
• --PAR=value: Use pygmt.config instead.

Notes on Input Formats

• data: Accepts a file path, 2-D numpy.ndarray, or pandas.DataFrame with (lon, lat) in the first two columns. Alternatively, provide x and y as separate 1-D arrays.
• outgrid: If not set, returns an xarray.DataArray; if set to a file path, writes the grid to disk and returns None.
• Both spacing and region are required; a GMTParameterError is raised if either is missing.
• The quantity="z" mode implements natural nearest-neighbour interpolation on the sphere, which is distinct from the planar pygmt.nearneighbor.

Linked Pull Requests

• Initial feature implementation – Wrap sphdistance #1383
• Add inline docstring example
• Add common data I/O options (-b, -d, -e, -f, -g, -h, -i, -j, -q, -r, -s, -w)
• Add a gallery or tutorial example (e.g., computing the distance to the nearest ocean island from any location on Earth, or filling a global grid by natural nearest-neighbour interpolation from sparse station data)

Related Issues and Discussions

• pygmt.sphdistance is the spherical counterpart to Cartesian Voronoi gridding; for Cartesian data or small-area geographic data, use pygmt.triangulate.regular_grid with the Voronoi option instead.
• The voronoi parameter (-Q) allows reuse of a pre-computed Voronoi tessellation (e.g., from pygmt.sphtriangulate) to speed up repeated gridding of the same point set.
• Setting quantity="n" produces a grid of polygon IDs useful for nearest-neighbour assignment tasks (e.g., associating grid cells to the closest weather station).