Add function to import hotspot dataset

doc/api/index.rst

@@ -95,6 +95,7 @@ Operations on grids:
     grdlandmask
     grdgradient
     grdtrack
+    sphdistance
 
 Crossover analysis with x2sys:
 

pygmt/__init__.py

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

pygmt/datasets/__init__.py

@@ -5,6 +5,7 @@
 from pygmt.datasets.earth_relief import load_earth_relief
 from pygmt.datasets.samples import (
     load_fractures_compilation,
+    load_hotspots,
     load_japan_quakes,
     load_ocean_ridge_points,
     load_sample_bathymetry,

pygmt/datasets/samples.py

@@ -123,3 +123,38 @@ def load_fractures_compilation():
     fname = which("@fractures_06.txt", download="c")
     data = pd.read_csv(fname, header=None, sep=r"\s+", names=["azimuth", "length"])
     return data[["length", "azimuth"]]
+
+
+def load_hotspots():
+    """
+    Load a table with the locations, names, and suggested icon sizes of
+    hotspots.
+
+    This is the ``@hotspots.txt`` 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. Use ``print(data.describe())`` to see the available
+        columns.
+    """
+    fname = which("@hotspots.txt", download="c")
+    with open(fname) as hotspot_text:
+        hotspot_text.readline()
+        hotspot_text.readline()
+        hotspot_text.readline()
+        hotspots = []
+        for line in hotspot_text:
+            line_split = line.strip().split("\t")
+            # Add coordinates and icon_size of hotspot
+            hotspot = [float(item.strip()) for item in line_split[0].split()]
+            hotspot.append(line_split[1].title())  # Add name of hotspot
+            hotspots.append(hotspot)
+    data = pd.DataFrame(
+        hotspots, columns=["longitude", "latitude", "icon_size", "name"]
+    )
+    return data

pygmt/src/__init__.py

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

@@ -0,0 +1,79 @@
+"""
+sphdistance - Create Voronoi distance, node,
+or natural nearest-neighbor grid on a sphere
+"""
+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(
+    G="outgrid",
+    I="increment",
+    R="region",
+)
+@kwargs_to_strings(I="sequence", R="sequence")
+def sphdistance(table, **kwargs):
+    r"""
+    Create Voroni polygons from lat/long coordinates.
+
+    Reads one or more ASCII [or binary] files (or standard
+    input) containing lon, lat and performs the construction of Voronoi
+    polygons. These polygons are then processed to calculate the nearest
+    distance to each node of the lattice and written to 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}
+
+    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``)
+    """
+    if "I" not in kwargs.keys() or "R" not in kwargs.keys():
+        raise GMTInvalidInput("Both 'region' and 'increment' must be specified.")
+    kind = data_kind(table)
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        with Session() as lib:
+            if kind == "file":
+                file_context = dummy_context(table)
+            elif kind == "matrix":
+                file_context = lib.virtualfile_from_matrix(matrix=table)
+            else:
+                raise GMTInvalidInput("Unrecognized data type: {}".format(type(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("sphdistance", 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_datasets_samples.py

@@ -1,8 +1,10 @@
 """
 Test basic functionality for loading sample datasets.
 """
+import pandas as pd
 from pygmt.datasets import (
     load_fractures_compilation,
+    load_hotspots,
     load_japan_quakes,
     load_ocean_ridge_points,
     load_sample_bathymetry,
@@ -72,3 +74,18 @@ def test_fractures_compilation():
     assert summary.loc["max", "length"] == 984.652
     assert summary.loc["min", "azimuth"] == 0.0
     assert summary.loc["max", "azimuth"] == 360.0
+
+
+def test_hotspots():
+    """
+    Check that the @hotspots.txt dataset loads without errors.
+    """
+    data = load_hotspots()
+    assert data.shape == (55, 4)
+    assert data.columns.values.tolist() == [
+        "longitude",
+        "latitude",
+        "icon_size",
+        "name",
+    ]
+    assert isinstance(data, pd.DataFrame)

pygmt/tests/test_sphdistance.py

@@ -0,0 +1,57 @@
+"""
+Tests for sphdistance.
+"""
+import os
+
+import numpy as np
+import pytest
+from pygmt import grdinfo, sphdistance
+from pygmt.exceptions import GMTInvalidInput
+from pygmt.helpers import GMTTempFile
+
+
+@pytest.fixture(scope="module", name="array")
+def fixture_table():
+    """
+    Load the table data.
+    """
+    coords_list = [[85.5, 22.3], [82.3, 22.6], [85.8, 22.4], [86.5, 23.3]]
+    return np.array(coords_list)
+
+
+def test_sphdistance_outgrid(array):
+    """
+    Test sphdistance with a set outgrid.
+    """
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        result = sphdistance(
+            table=array, outgrid=tmpfile.name, increment=1, region=[82, 87, 22, 24]
+        )
+        assert result is None  # return value is None
+        assert os.path.exists(path=tmpfile.name)  # check that outgrid exists
+
+
+def test_sphdistance_no_outgrid(array):
+    """
+    Test sphdistance with no set outgrid.
+    """
+    temp_grid = sphdistance(table=array, increment=[1, 2], region=[82, 87, 22, 24])
+    assert temp_grid.dims == ("lat", "lon")
+    assert temp_grid.gmt.gtype == 1  # Geographic grid
+    assert temp_grid.gmt.registration == 0  # Gridline registration
+    result = grdinfo(grid=temp_grid, force_scan="a", per_column="n").strip().split()
+    assert int(result[0]) == 82  # x minimum
+    assert int(result[1]) == 87  # x maximum
+    assert int(result[2]) == 22  # y minimum
+    assert int(result[3]) == 24  # y maximum
+    assert int(result[6]) == 1  # x increment
+    assert int(result[7]) == 2  # y increment
+
+
+def test_sphdistance_fails(array):
+    """
+    Check that sphdistance fails correctly when neither increment nor region is
+    given.
+    """
+    with pytest.raises(GMTInvalidInput):
+        sphdistance(table=array)