From a6cef02b485b042c50a3614800ca582c50dc4b33 Mon Sep 17 00:00:00 2001
From: Wei Ji <weiji.leong@vuw.ac.nz>
Date: Wed, 22 May 2019 14:28:14 +0200
Subject: [PATCH 01/11] Wrap grdtrack

Initial commit for #307. Implement GMT grdtrack function under sampling.py. Test cases checking proper pandas.DataFrame and xarray.DataArray inputs stored in test_grdtrack.py. Sample datasets for tests uses newly created load_east_pacific_rise_grid and load_ocean_ridge_points functions in datasets/tutorial.py.

GMT grdtrack documentation can be found at https://gmt.soest.hawaii.edu/doc/latest/grdtrack. Originally, grdtrack should take in an xyfile and -Ggridfile as parameters, and pass the output table to stdout. Here, the implementation (currently) takes in a pandas.DataFrame table and xarray.DataArray grid instead as input, and returns a pandas.DataFrame with an extra column for the sampled grid data.
---
 pygmt/__init__.py            |  1 +
 pygmt/datasets/__init__.py   |  8 ++++-
 pygmt/datasets/tutorial.py   | 45 +++++++++++++++++++++++
 pygmt/sampling.py            | 70 ++++++++++++++++++++++++++++++++++++
 pygmt/tests/test_grdtrack.py | 66 ++++++++++++++++++++++++++++++++++
 5 files changed, 189 insertions(+), 1 deletion(-)
 create mode 100644 pygmt/sampling.py
 create mode 100644 pygmt/tests/test_grdtrack.py

diff --git a/pygmt/__init__.py b/pygmt/__init__.py
index 9c40f485c36..031521d6f5e 100644
--- a/pygmt/__init__.py
+++ b/pygmt/__init__.py
@@ -15,6 +15,7 @@
 from .session_management import begin as _begin, end as _end
 from .figure import Figure
 from .gridding import surface
+from .sampling import grdtrack
 from .modules import info, grdinfo, which
 from . import datasets
 
diff --git a/pygmt/datasets/__init__.py b/pygmt/datasets/__init__.py
index 44d6229fcee..17b4020b900 100644
--- a/pygmt/datasets/__init__.py
+++ b/pygmt/datasets/__init__.py
@@ -2,5 +2,11 @@
 #
 # Load sample data included with GMT (downloaded from the GMT cache server).
 
-from .tutorial import load_japan_quakes, load_sample_bathymetry, load_usgs_quakes
+from .tutorial import (
+    load_east_pacific_rise_grid,
+    load_japan_quakes,
+    load_ocean_ridge_points,
+    load_sample_bathymetry,
+    load_usgs_quakes,
+)
 from .earth_relief import load_earth_relief
diff --git a/pygmt/datasets/tutorial.py b/pygmt/datasets/tutorial.py
index bee3c18609d..23e0aebe8f9 100644
--- a/pygmt/datasets/tutorial.py
+++ b/pygmt/datasets/tutorial.py
@@ -2,10 +2,33 @@
 Functions to load sample data from the GMT tutorials.
 """
 import pandas as pd
+import xarray as xr
 
 from .. import which
 
 
+def load_east_pacific_rise_grid():
+    """
+    Load a grid of bathymetry over part of the East Pacific Rise as a xarray.DataArray.
+
+    This is the ``@spac_33.nc`` 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 : xarray.DataArray
+        The data grid. Coordinates in longitude (lon) and latitude (lat).
+        Data attributes: bathymetry (z) in metres.
+    """
+    fname = which("@spac_33.nc", download="c")
+    with xr.open_dataarray(fname) as dataarray:
+        data = dataarray.load()
+    return data
+
+
 def load_japan_quakes():
     """
     Load a table of earthquakes around Japan as a pandas.Dataframe.
@@ -38,6 +61,28 @@ def load_japan_quakes():
     return data
 
 
+def load_ocean_ridge_points():
+    """
+    Load a table of ocean ridge points for the entire world as a pandas.DataFrame.
+
+    This is the ``@ridge.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. Columns are longitude and latitude.
+    """
+    fname = which("@ridge.txt", download="c")
+    data = pd.read_csv(
+        fname, sep=r"\s+", names=["longitude", "latitude"], skiprows=1, comment=">"
+    )
+    return data
+
+
 def load_sample_bathymetry():
     """
     Load a table of ship observations of bathymetry off Baja California as a
diff --git a/pygmt/sampling.py b/pygmt/sampling.py
new file mode 100644
index 00000000000..81d91f7ffa2
--- /dev/null
+++ b/pygmt/sampling.py
@@ -0,0 +1,70 @@
+"""
+GMT modules for Sampling of 1-D and 2-D Data
+"""
+import pandas as pd
+import xarray as xr
+
+from .clib import Session
+from .helpers import build_arg_string, fmt_docstring, GMTTempFile, data_kind
+from .exceptions import GMTInvalidInput
+
+
+@fmt_docstring
+def grdtrack(table: pd.DataFrame, grid: xr.DataArray, newcolname: str = "z_", **kwargs):
+    """
+    Sample grids at specified (x,y) locations.
+
+    Grdtrack reads one or more grid files and a table with (x,y) [or (lon,lat)]
+    positions in the first two columns (more columns may be present). It interpolates
+    the grid(s) at the positions in the table and writes out the table with the
+    interpolated values added as (one or more) new columns. A bicubic [Default],
+    bilinear, B-spline or nearest-neighbor (see -n) interpolation is used, requiring
+    boundary conditions at the limits of the region.
+
+    Parameters
+    ----------
+    table: pandas.DataFrame
+        Table with (x, y) or (lon, lat) values in the first two columns. More columns
+        may be present.
+
+    grid: xarray.DataArray
+        Gridded array from which to sample values from.
+
+    newcolname: str
+        Name for the new column in the table where the sampled values will be placed.
+        Defaults to "z_".
+
+    Returns
+    -------
+    ret: pandas.DataFrame
+        Table with (x, y, ..., z_) or (lon, lat, ..., z_) values.
+
+    """
+    with GMTTempFile(suffix=".csv") as tmpfile:
+        with Session() as lib:
+            # Store the pandas.DataFrame table in virtualfile
+            if data_kind(table) == "matrix":
+                table_context = lib.virtualfile_from_matrix(table.values)
+            else:
+                raise GMTInvalidInput(f"Unrecognized data type {type(table)}")
+
+            # Store the xarray.DataArray grid in virtualfile
+            if data_kind(grid) == "grid":
+                grid_context = lib.virtualfile_from_grid(grid)
+            else:
+                raise GMTInvalidInput(f"Unrecognized data type {type(grid)}")
+
+            # Run grdtrack on the temporary (csv) table and (netcdf) grid virtualfiles
+            with table_context as csvfile:
+                with grid_context as grdfile:
+                    kwargs = {"G": grdfile}
+                    arg_str = " ".join(
+                        [csvfile, build_arg_string(kwargs), "->" + tmpfile.name]
+                    )
+                    lib.call_module(module="grdtrack", args=arg_str)
+
+        # Read temporary csv output to a pandas table
+        column_names = table.columns.to_list() + [newcolname]
+        result = pd.read_csv(tmpfile.name, sep="\t", names=column_names)
+
+    return result
diff --git a/pygmt/tests/test_grdtrack.py b/pygmt/tests/test_grdtrack.py
new file mode 100644
index 00000000000..c18c7b88131
--- /dev/null
+++ b/pygmt/tests/test_grdtrack.py
@@ -0,0 +1,66 @@
+"""
+Tests for grdtrack
+"""
+
+import pandas as pd
+import pytest
+
+from .. import grdtrack
+from ..datasets import load_east_pacific_rise_grid, load_ocean_ridge_points
+from ..exceptions import GMTInvalidInput
+from ..helpers import data_kind
+
+
+def test_grdtrack_input_dataframe_and_dataarray():
+    """
+    Run grdtrack by passing in a pandas.DataFrame and xarray.DataArray as inputs
+    """
+    dataframe = load_ocean_ridge_points()
+    dataarray = load_east_pacific_rise_grid()
+
+    output = grdtrack(table=dataframe, grid=dataarray)
+    assert isinstance(output, pd.DataFrame)
+    assert output.columns.to_list() == ["longitude", "latitude", "z_"]
+    assert output.iloc[0].to_list() == [-110.9536, -42.2489, -2950.49576833]
+
+    return output
+
+
+def test_grdtrack_input_wrong_kind_of_table():
+    """
+    Run grdtrack using table input that is not a pandas.DataFrame (matrix)
+    """
+    dataframe = load_ocean_ridge_points()
+    invalid_table = dataframe.longitude.to_xarray()
+    dataarray = load_east_pacific_rise_grid()
+
+    assert data_kind(invalid_table) == "grid"
+    with pytest.raises(GMTInvalidInput):
+        grdtrack(table=invalid_table, grid=dataarray)
+
+
+def test_grdtrack_input_wrong_kind_of_grid():
+    """
+    Run grdtrack using grid input that is not an xarray.DataArray (grid)
+    """
+    dataframe = load_ocean_ridge_points()
+    dataarray = load_east_pacific_rise_grid()
+    invalid_grid = dataarray.to_dataset()
+
+    assert data_kind(invalid_grid) == "matrix"
+    with pytest.raises(GMTInvalidInput):
+        grdtrack(table=dataframe, grid=invalid_grid)
+
+
+def test_grdtrack_newcolname_setting():
+    """
+    Run grdtrack by passing in a non-default newcolname parameter setting
+    """
+    dataframe = load_ocean_ridge_points()
+    dataarray = load_east_pacific_rise_grid()
+
+    output = grdtrack(table=dataframe, grid=dataarray, newcolname="bathymetry")
+    assert output.columns.to_list() == ["longitude", "latitude", "bathymetry"]
+    assert output.iloc[0].to_list() == [-110.9536, -42.2489, -2950.49576833]
+
+    return output

From 97c2b4ac57068525956e16e7f2cf9f01c8e29fc5 Mon Sep 17 00:00:00 2001
From: Wei Ji <weiji.leong@vuw.ac.nz>
Date: Wed, 22 May 2019 17:28:37 +0200
Subject: [PATCH 02/11] Enable netcdf file input for grdtrack

Enable passing in netcdf file as input into grdtrack, instead of just xarray.DataArray. Also fix bug with kwargs being overwritten instead of appended to...
---
 pygmt/sampling.py            | 14 +++++++++++---
 pygmt/tests/test_grdtrack.py | 17 ++++++++++++++++-
 2 files changed, 27 insertions(+), 4 deletions(-)

diff --git a/pygmt/sampling.py b/pygmt/sampling.py
index 81d91f7ffa2..f96c5902e2c 100644
--- a/pygmt/sampling.py
+++ b/pygmt/sampling.py
@@ -5,7 +5,13 @@
 import xarray as xr
 
 from .clib import Session
-from .helpers import build_arg_string, fmt_docstring, GMTTempFile, data_kind
+from .helpers import (
+    build_arg_string,
+    fmt_docstring,
+    GMTTempFile,
+    data_kind,
+    dummy_context,
+)
 from .exceptions import GMTInvalidInput
 
 
@@ -27,7 +33,7 @@ def grdtrack(table: pd.DataFrame, grid: xr.DataArray, newcolname: str = "z_", **
         Table with (x, y) or (lon, lat) values in the first two columns. More columns
         may be present.
 
-    grid: xarray.DataArray
+    grid: xarray.DataArray or file (netcdf)
         Gridded array from which to sample values from.
 
     newcolname: str
@@ -51,13 +57,15 @@ def grdtrack(table: pd.DataFrame, grid: xr.DataArray, newcolname: str = "z_", **
             # Store the xarray.DataArray grid in virtualfile
             if data_kind(grid) == "grid":
                 grid_context = lib.virtualfile_from_grid(grid)
+            elif data_kind(grid) == "file":
+                grid_context = dummy_context(grid)
             else:
                 raise GMTInvalidInput(f"Unrecognized data type {type(grid)}")
 
             # Run grdtrack on the temporary (csv) table and (netcdf) grid virtualfiles
             with table_context as csvfile:
                 with grid_context as grdfile:
-                    kwargs = {"G": grdfile}
+                    kwargs.update({"G": grdfile})
                     arg_str = " ".join(
                         [csvfile, build_arg_string(kwargs), "->" + tmpfile.name]
                     )
diff --git a/pygmt/tests/test_grdtrack.py b/pygmt/tests/test_grdtrack.py
index c18c7b88131..4ceeb666d2d 100644
--- a/pygmt/tests/test_grdtrack.py
+++ b/pygmt/tests/test_grdtrack.py
@@ -6,6 +6,7 @@
 import pytest
 
 from .. import grdtrack
+from .. import which
 from ..datasets import load_east_pacific_rise_grid, load_ocean_ridge_points
 from ..exceptions import GMTInvalidInput
 from ..helpers import data_kind
@@ -26,6 +27,20 @@ def test_grdtrack_input_dataframe_and_dataarray():
     return output
 
 
+def test_grdtrack_input_dataframe_and_ncfile():
+    """
+    Run grdtrack by passing in a pandas.DataFrame and netcdf file as inputs
+    """
+    dataframe = load_ocean_ridge_points()
+    ncfile = which("@spac_33.nc", download="c")
+
+    output = grdtrack(table=dataframe, grid=ncfile)
+    assert isinstance(output, pd.DataFrame)
+    assert output.columns.to_list() == ["longitude", "latitude", "z_"]
+
+    return output
+
+
 def test_grdtrack_input_wrong_kind_of_table():
     """
     Run grdtrack using table input that is not a pandas.DataFrame (matrix)
@@ -41,7 +56,7 @@ def test_grdtrack_input_wrong_kind_of_table():
 
 def test_grdtrack_input_wrong_kind_of_grid():
     """
-    Run grdtrack using grid input that is not an xarray.DataArray (grid)
+    Run grdtrack using grid input that is not as xarray.DataArray (grid) or file
     """
     dataframe = load_ocean_ridge_points()
     dataarray = load_east_pacific_rise_grid()

From bcfab7f6154ecea7a65b67304f94045292382533 Mon Sep 17 00:00:00 2001
From: Wei Ji <weiji.leong@vuw.ac.nz>
Date: Thu, 23 May 2019 20:46:35 +0200
Subject: [PATCH 03/11] Modify grdtrack parameter names according to code
 review

Rename 'table' parameter to 'points', remove default 'newcolname' parameter (require it to be explicitly set) and change the output 'ret' to 'track' in docstring. Unit tests updated accordingly.
---
 pygmt/sampling.py            | 29 ++++++++++++++++++-----------
 pygmt/tests/test_grdtrack.py | 34 ++++++++++++++++------------------
 2 files changed, 34 insertions(+), 29 deletions(-)

diff --git a/pygmt/sampling.py b/pygmt/sampling.py
index f96c5902e2c..28d1439c7f7 100644
--- a/pygmt/sampling.py
+++ b/pygmt/sampling.py
@@ -16,7 +16,9 @@
 
 
 @fmt_docstring
-def grdtrack(table: pd.DataFrame, grid: xr.DataArray, newcolname: str = "z_", **kwargs):
+def grdtrack(
+    points: pd.DataFrame, grid: xr.DataArray, newcolname: str = None, **kwargs
+):
     """
     Sample grids at specified (x,y) locations.
 
@@ -29,7 +31,7 @@ def grdtrack(table: pd.DataFrame, grid: xr.DataArray, newcolname: str = "z_", **
 
     Parameters
     ----------
-    table: pandas.DataFrame
+    points: pandas.DataFrame
         Table with (x, y) or (lon, lat) values in the first two columns. More columns
         may be present.
 
@@ -38,21 +40,26 @@ def grdtrack(table: pd.DataFrame, grid: xr.DataArray, newcolname: str = "z_", **
 
     newcolname: str
         Name for the new column in the table where the sampled values will be placed.
-        Defaults to "z_".
 
     Returns
     -------
-    ret: pandas.DataFrame
-        Table with (x, y, ..., z_) or (lon, lat, ..., z_) values.
+    track: pandas.DataFrame
+        Table with (x, y, ..., newcolname) or (lon, lat, ..., newcolname) values.
 
     """
+
+    try:
+        assert isinstance(newcolname, str)
+    except AssertionError:
+        raise GMTInvalidInput("Please pass in a str to 'newcolname'")
+
     with GMTTempFile(suffix=".csv") as tmpfile:
         with Session() as lib:
-            # Store the pandas.DataFrame table in virtualfile
-            if data_kind(table) == "matrix":
-                table_context = lib.virtualfile_from_matrix(table.values)
+            # Store the pandas.DataFrame points table in virtualfile
+            if data_kind(points) == "matrix":
+                table_context = lib.virtualfile_from_matrix(points.values)
             else:
-                raise GMTInvalidInput(f"Unrecognized data type {type(table)}")
+                raise GMTInvalidInput(f"Unrecognized data type {type(points)}")
 
             # Store the xarray.DataArray grid in virtualfile
             if data_kind(grid) == "grid":
@@ -62,7 +69,7 @@ def grdtrack(table: pd.DataFrame, grid: xr.DataArray, newcolname: str = "z_", **
             else:
                 raise GMTInvalidInput(f"Unrecognized data type {type(grid)}")
 
-            # Run grdtrack on the temporary (csv) table and (netcdf) grid virtualfiles
+            # Run grdtrack on the temp (csv) points table and (netcdf) grid virtualfiles
             with table_context as csvfile:
                 with grid_context as grdfile:
                     kwargs.update({"G": grdfile})
@@ -72,7 +79,7 @@ def grdtrack(table: pd.DataFrame, grid: xr.DataArray, newcolname: str = "z_", **
                     lib.call_module(module="grdtrack", args=arg_str)
 
         # Read temporary csv output to a pandas table
-        column_names = table.columns.to_list() + [newcolname]
+        column_names = points.columns.to_list() + [newcolname]
         result = pd.read_csv(tmpfile.name, sep="\t", names=column_names)
 
     return result
diff --git a/pygmt/tests/test_grdtrack.py b/pygmt/tests/test_grdtrack.py
index 4ceeb666d2d..5525178f538 100644
--- a/pygmt/tests/test_grdtrack.py
+++ b/pygmt/tests/test_grdtrack.py
@@ -19,9 +19,9 @@ def test_grdtrack_input_dataframe_and_dataarray():
     dataframe = load_ocean_ridge_points()
     dataarray = load_east_pacific_rise_grid()
 
-    output = grdtrack(table=dataframe, grid=dataarray)
+    output = grdtrack(points=dataframe, grid=dataarray, newcolname="bathymetry")
     assert isinstance(output, pd.DataFrame)
-    assert output.columns.to_list() == ["longitude", "latitude", "z_"]
+    assert output.columns.to_list() == ["longitude", "latitude", "bathymetry"]
     assert output.iloc[0].to_list() == [-110.9536, -42.2489, -2950.49576833]
 
     return output
@@ -34,27 +34,28 @@ def test_grdtrack_input_dataframe_and_ncfile():
     dataframe = load_ocean_ridge_points()
     ncfile = which("@spac_33.nc", download="c")
 
-    output = grdtrack(table=dataframe, grid=ncfile)
+    output = grdtrack(points=dataframe, grid=ncfile, newcolname="bathymetry")
     assert isinstance(output, pd.DataFrame)
-    assert output.columns.to_list() == ["longitude", "latitude", "z_"]
+    assert output.columns.to_list() == ["longitude", "latitude", "bathymetry"]
+    assert output.iloc[0].to_list() == [-110.9536, -42.2489, -2950.49576833]
 
     return output
 
 
-def test_grdtrack_input_wrong_kind_of_table():
+def test_grdtrack_wrong_kind_of_points_input():
     """
-    Run grdtrack using table input that is not a pandas.DataFrame (matrix)
+    Run grdtrack using points input that is not a pandas.DataFrame (matrix)
     """
     dataframe = load_ocean_ridge_points()
-    invalid_table = dataframe.longitude.to_xarray()
+    invalid_points = dataframe.longitude.to_xarray()
     dataarray = load_east_pacific_rise_grid()
 
-    assert data_kind(invalid_table) == "grid"
+    assert data_kind(invalid_points) == "grid"
     with pytest.raises(GMTInvalidInput):
-        grdtrack(table=invalid_table, grid=dataarray)
+        grdtrack(points=invalid_points, grid=dataarray, newcolname="bathymetry")
 
 
-def test_grdtrack_input_wrong_kind_of_grid():
+def test_grdtrack_wrong_kind_of_grid_input():
     """
     Run grdtrack using grid input that is not as xarray.DataArray (grid) or file
     """
@@ -64,18 +65,15 @@ def test_grdtrack_input_wrong_kind_of_grid():
 
     assert data_kind(invalid_grid) == "matrix"
     with pytest.raises(GMTInvalidInput):
-        grdtrack(table=dataframe, grid=invalid_grid)
+        grdtrack(points=dataframe, grid=invalid_grid, newcolname="bathymetry")
 
 
-def test_grdtrack_newcolname_setting():
+def test_grdtrack_without_newcolname_setting():
     """
-    Run grdtrack by passing in a non-default newcolname parameter setting
+    Run grdtrack by not passing in newcolname parameter setting
     """
     dataframe = load_ocean_ridge_points()
     dataarray = load_east_pacific_rise_grid()
 
-    output = grdtrack(table=dataframe, grid=dataarray, newcolname="bathymetry")
-    assert output.columns.to_list() == ["longitude", "latitude", "bathymetry"]
-    assert output.iloc[0].to_list() == [-110.9536, -42.2489, -2950.49576833]
-
-    return output
+    with pytest.raises(GMTInvalidInput):
+        grdtrack(points=dataframe, grid=dataarray)

From 028b55a35fd6130661ec97e109d0c660551d6524 Mon Sep 17 00:00:00 2001
From: Wei Ji <weiji.leong@vuw.ac.nz>
Date: Thu, 23 May 2019 21:55:55 +0200
Subject: [PATCH 04/11] Remove east_pacific_rise_grid, add test for
 load_ocean_ridge_points

One less dataset to test by replacing the east_pacific_rise_grid with the earth_relief 60m grid (sliced to the same area). Also add tests for loading the ocean_ridge_points dataset.
---
 pygmt/datasets/__init__.py   |  1 -
 pygmt/datasets/tutorial.py   | 22 ----------------------
 pygmt/tests/test_datasets.py | 12 ++++++++++++
 pygmt/tests/test_grdtrack.py | 16 ++++++++--------
 4 files changed, 20 insertions(+), 31 deletions(-)

diff --git a/pygmt/datasets/__init__.py b/pygmt/datasets/__init__.py
index 17b4020b900..4454903ee11 100644
--- a/pygmt/datasets/__init__.py
+++ b/pygmt/datasets/__init__.py
@@ -3,7 +3,6 @@
 # Load sample data included with GMT (downloaded from the GMT cache server).
 
 from .tutorial import (
-    load_east_pacific_rise_grid,
     load_japan_quakes,
     load_ocean_ridge_points,
     load_sample_bathymetry,
diff --git a/pygmt/datasets/tutorial.py b/pygmt/datasets/tutorial.py
index 23e0aebe8f9..10e394a213a 100644
--- a/pygmt/datasets/tutorial.py
+++ b/pygmt/datasets/tutorial.py
@@ -7,28 +7,6 @@
 from .. import which
 
 
-def load_east_pacific_rise_grid():
-    """
-    Load a grid of bathymetry over part of the East Pacific Rise as a xarray.DataArray.
-
-    This is the ``@spac_33.nc`` 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 : xarray.DataArray
-        The data grid. Coordinates in longitude (lon) and latitude (lat).
-        Data attributes: bathymetry (z) in metres.
-    """
-    fname = which("@spac_33.nc", download="c")
-    with xr.open_dataarray(fname) as dataarray:
-        data = dataarray.load()
-    return data
-
-
 def load_japan_quakes():
     """
     Load a table of earthquakes around Japan as a pandas.Dataframe.
diff --git a/pygmt/tests/test_datasets.py b/pygmt/tests/test_datasets.py
index 23594de5031..24c1c721f91 100644
--- a/pygmt/tests/test_datasets.py
+++ b/pygmt/tests/test_datasets.py
@@ -8,6 +8,7 @@
 from ..datasets import (
     load_japan_quakes,
     load_earth_relief,
+    load_ocean_ridge_points,
     load_sample_bathymetry,
     load_usgs_quakes,
 )
@@ -27,6 +28,17 @@ def test_japan_quakes():
     assert summary.loc["max", "day"] == 31
 
 
+def test_ocean_ridge_points():
+    "Check that the @ridge.txt dataset loads without errors"
+    data = load_ocean_ridge_points()
+    assert data.shape == (4146, 2)
+    summary = data.describe()
+    assert summary.loc["min", "longitude"] == -179.9401
+    assert summary.loc["max", "longitude"] == 179.935
+    assert summary.loc["min", "latitude"] == -65.6182
+    assert summary.loc["max", "latitude"] == 86.8
+
+
 def test_sample_bathymetry():
     "Check that the @tut_ship.xyz dataset loads without errors"
     data = load_sample_bathymetry()
diff --git a/pygmt/tests/test_grdtrack.py b/pygmt/tests/test_grdtrack.py
index 5525178f538..fbf3b972dce 100644
--- a/pygmt/tests/test_grdtrack.py
+++ b/pygmt/tests/test_grdtrack.py
@@ -7,7 +7,7 @@
 
 from .. import grdtrack
 from .. import which
-from ..datasets import load_east_pacific_rise_grid, load_ocean_ridge_points
+from ..datasets import load_earth_relief, load_ocean_ridge_points
 from ..exceptions import GMTInvalidInput
 from ..helpers import data_kind
 
@@ -17,12 +17,12 @@ def test_grdtrack_input_dataframe_and_dataarray():
     Run grdtrack by passing in a pandas.DataFrame and xarray.DataArray as inputs
     """
     dataframe = load_ocean_ridge_points()
-    dataarray = load_east_pacific_rise_grid()
+    dataarray = load_earth_relief().sel(lat=slice(-49, -42), lon=slice(-118, -107))
 
     output = grdtrack(points=dataframe, grid=dataarray, newcolname="bathymetry")
     assert isinstance(output, pd.DataFrame)
     assert output.columns.to_list() == ["longitude", "latitude", "bathymetry"]
-    assert output.iloc[0].to_list() == [-110.9536, -42.2489, -2950.49576833]
+    assert output.iloc[0].to_list() == [-110.9536, -42.2489, -2823.96637605]
 
     return output
 
@@ -32,12 +32,12 @@ def test_grdtrack_input_dataframe_and_ncfile():
     Run grdtrack by passing in a pandas.DataFrame and netcdf file as inputs
     """
     dataframe = load_ocean_ridge_points()
-    ncfile = which("@spac_33.nc", download="c")
+    ncfile = which("@earth_relief_60m", download="c")
 
     output = grdtrack(points=dataframe, grid=ncfile, newcolname="bathymetry")
     assert isinstance(output, pd.DataFrame)
     assert output.columns.to_list() == ["longitude", "latitude", "bathymetry"]
-    assert output.iloc[0].to_list() == [-110.9536, -42.2489, -2950.49576833]
+    assert output.iloc[0].to_list() == [-32.2971, 37.4118, -1697.87197487]
 
     return output
 
@@ -48,7 +48,7 @@ def test_grdtrack_wrong_kind_of_points_input():
     """
     dataframe = load_ocean_ridge_points()
     invalid_points = dataframe.longitude.to_xarray()
-    dataarray = load_east_pacific_rise_grid()
+    dataarray = load_earth_relief().sel(lat=slice(-49, -42), lon=slice(-118, -107))
 
     assert data_kind(invalid_points) == "grid"
     with pytest.raises(GMTInvalidInput):
@@ -60,7 +60,7 @@ def test_grdtrack_wrong_kind_of_grid_input():
     Run grdtrack using grid input that is not as xarray.DataArray (grid) or file
     """
     dataframe = load_ocean_ridge_points()
-    dataarray = load_east_pacific_rise_grid()
+    dataarray = load_earth_relief().sel(lat=slice(-49, -42), lon=slice(-118, -107))
     invalid_grid = dataarray.to_dataset()
 
     assert data_kind(invalid_grid) == "matrix"
@@ -73,7 +73,7 @@ def test_grdtrack_without_newcolname_setting():
     Run grdtrack by not passing in newcolname parameter setting
     """
     dataframe = load_ocean_ridge_points()
-    dataarray = load_east_pacific_rise_grid()
+    dataarray = load_earth_relief().sel(lat=slice(-49, -42), lon=slice(-118, -107))
 
     with pytest.raises(GMTInvalidInput):
         grdtrack(points=dataframe, grid=dataarray)

From 3776214536a2433b348bb948f56434150dd44d44 Mon Sep 17 00:00:00 2001
From: Wei Ji <weiji.leong@vuw.ac.nz>
Date: Thu, 23 May 2019 22:12:53 +0200
Subject: [PATCH 05/11] Add grdtrack and load_ocean_ridge_points tutorial data
 to docs

PyGMT grdtrack added under "Operations on grids" (though it does require a tabular data input). Also added the ocean_ridge_points tutorial dataset used in grdtrack's unit tests and reordered the tutorial datasets alphabetically.
---
 doc/api/index.rst | 7 ++++---
 pygmt/sampling.py | 2 ++
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/doc/api/index.rst b/doc/api/index.rst
index 3b2a4aeab6d..adbdb1886b4 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -59,6 +59,7 @@ Operations on grids:
     :toctree: generated
 
     grdinfo
+    grdtrack
 
 
 Miscellaneous
@@ -87,10 +88,10 @@ and store them in the GMT cache folder.
     :toctree: generated
 
     datasets.load_earth_relief
-    datasets.load_usgs_quakes
-    datasets.load_sample_bathymetry
     datasets.load_japan_quakes
-
+    datasets.load_ocean_ridge_points
+    datasets.load_sample_bathymetry
+    datasets.load_usgs_quakes
 
 .. automodule:: pygmt.exceptions
 
diff --git a/pygmt/sampling.py b/pygmt/sampling.py
index 28d1439c7f7..9c988385089 100644
--- a/pygmt/sampling.py
+++ b/pygmt/sampling.py
@@ -29,6 +29,8 @@ def grdtrack(
     bilinear, B-spline or nearest-neighbor (see -n) interpolation is used, requiring
     boundary conditions at the limits of the region.
 
+    Full option list at :gmt-docs:`grdtrack.html`
+
     Parameters
     ----------
     points: pandas.DataFrame

From a37fb0aa3fef025eb1c723039060caa0489fcab1 Mon Sep 17 00:00:00 2001
From: Wei Ji <weiji.leong@vuw.ac.nz>
Date: Mon, 27 May 2019 21:45:03 +0200
Subject: [PATCH 06/11] Enable ascii file input for grdtrack

Enable passing in ascii file inputs (csv, txt, etc) into grdtrack's 'points' parameter instead of just pandas.DataFrame. This requires a new 'outfile' parameter to be set. The type of 'points' input determines the type of 'track' returned, i.e. pd.DataFrame in, pd.DataFrame out; filename in, filename out. Extra unit tests created to test the various new input combinations and associated outputs.
---
 pygmt/datasets/tutorial.py   |  1 -
 pygmt/sampling.py            | 45 ++++++++++++++++------------
 pygmt/tests/test_grdtrack.py | 57 +++++++++++++++++++++++++++++++++++-
 3 files changed, 83 insertions(+), 20 deletions(-)

diff --git a/pygmt/datasets/tutorial.py b/pygmt/datasets/tutorial.py
index 10e394a213a..3d9188e0dfa 100644
--- a/pygmt/datasets/tutorial.py
+++ b/pygmt/datasets/tutorial.py
@@ -2,7 +2,6 @@
 Functions to load sample data from the GMT tutorials.
 """
 import pandas as pd
-import xarray as xr
 
 from .. import which
 
diff --git a/pygmt/sampling.py b/pygmt/sampling.py
index 9c988385089..7df01eb2884 100644
--- a/pygmt/sampling.py
+++ b/pygmt/sampling.py
@@ -2,7 +2,6 @@
 GMT modules for Sampling of 1-D and 2-D Data
 """
 import pandas as pd
-import xarray as xr
 
 from .clib import Session
 from .helpers import (
@@ -16,9 +15,7 @@
 
 
 @fmt_docstring
-def grdtrack(
-    points: pd.DataFrame, grid: xr.DataArray, newcolname: str = None, **kwargs
-):
+def grdtrack(points, grid, newcolname=None, outfile=None, **kwargs):
     """
     Sample grids at specified (x,y) locations.
 
@@ -33,33 +30,40 @@ def grdtrack(
 
     Parameters
     ----------
-    points: pandas.DataFrame
-        Table with (x, y) or (lon, lat) values in the first two columns. More columns
-        may be present.
+    points: pandas.DataFrame or file (csv, txt, etc)
+        Either a table with (x, y) or (lon, lat) values in the first two columns,
+        or a data file name. More columns may be present.
 
     grid: xarray.DataArray or file (netcdf)
         Gridded array from which to sample values from.
 
     newcolname: str
-        Name for the new column in the table where the sampled values will be placed.
+        Required if 'points' is a pandas.DataFrame. The name for the new column in the
+        track pandas.DataFrame table where the sampled values will be placed.
+
+    outfile: str
+        Required if 'points' is a file. The file name for the output ASCII file.
 
     Returns
     -------
-    track: pandas.DataFrame
-        Table with (x, y, ..., newcolname) or (lon, lat, ..., newcolname) values.
+    track: pandas.DataFrame or None
+        Return type depends on whether the outfile parameter is set:
+        - pandas.DataFrame table with (x, y, ..., newcolname) if outfile is not set
+        - None if outfile is set (track output will be stored in outfile)
 
     """
 
-    try:
-        assert isinstance(newcolname, str)
-    except AssertionError:
-        raise GMTInvalidInput("Please pass in a str to 'newcolname'")
-
     with GMTTempFile(suffix=".csv") as tmpfile:
         with Session() as lib:
             # Store the pandas.DataFrame points table in virtualfile
             if data_kind(points) == "matrix":
+                if newcolname is None:
+                    raise GMTInvalidInput("Please pass in a str to 'newcolname'")
                 table_context = lib.virtualfile_from_matrix(points.values)
+            elif data_kind(points) == "file":
+                if outfile is None:
+                    raise GMTInvalidInput("Please pass in a str to 'outfile'")
+                table_context = dummy_context(points)
             else:
                 raise GMTInvalidInput(f"Unrecognized data type {type(points)}")
 
@@ -75,13 +79,18 @@ def grdtrack(
             with table_context as csvfile:
                 with grid_context as grdfile:
                     kwargs.update({"G": grdfile})
+                    if outfile is None:  # Output to tmpfile if outfile is not set
+                        outfile = tmpfile.name
                     arg_str = " ".join(
-                        [csvfile, build_arg_string(kwargs), "->" + tmpfile.name]
+                        [csvfile, build_arg_string(kwargs), "->" + outfile]
                     )
                     lib.call_module(module="grdtrack", args=arg_str)
 
         # Read temporary csv output to a pandas table
-        column_names = points.columns.to_list() + [newcolname]
-        result = pd.read_csv(tmpfile.name, sep="\t", names=column_names)
+        if outfile == tmpfile.name:  # if user did not set outfile, return pd.DataFrame
+            column_names = points.columns.to_list() + [newcolname]
+            result = pd.read_csv(tmpfile.name, sep="\t", names=column_names)
+        elif outfile != tmpfile.name:  # return None if outfile set, output in outfile
+            result = None
 
     return result
diff --git a/pygmt/tests/test_grdtrack.py b/pygmt/tests/test_grdtrack.py
index fbf3b972dce..4801ee0428f 100644
--- a/pygmt/tests/test_grdtrack.py
+++ b/pygmt/tests/test_grdtrack.py
@@ -1,6 +1,7 @@
 """
 Tests for grdtrack
 """
+import os
 
 import pandas as pd
 import pytest
@@ -11,6 +12,9 @@
 from ..exceptions import GMTInvalidInput
 from ..helpers import data_kind
 
+TEST_DATA_DIR = os.path.join(os.path.dirname(__file__), "data")
+TEMP_TRACK = os.path.join(TEST_DATA_DIR, "tmp_track.txt")
+
 
 def test_grdtrack_input_dataframe_and_dataarray():
     """
@@ -27,6 +31,26 @@ def test_grdtrack_input_dataframe_and_dataarray():
     return output
 
 
+def test_grdtrack_input_csvfile_and_dataarray():
+    """
+    Run grdtrack by passing in a csvfile and xarray.DataArray as inputs
+    """
+    csvfile = which("@ridge.txt", download="c")
+    dataarray = load_earth_relief().sel(lat=slice(-49, -42), lon=slice(-118, -107))
+
+    try:
+        output = grdtrack(points=csvfile, grid=dataarray, outfile=TEMP_TRACK)
+        assert output is None  # check that output is None since outfile is set
+        assert os.path.exists(path=TEMP_TRACK)  # check that outfile exists at path
+
+        track = pd.read_csv(TEMP_TRACK, sep="\t", header=None, comment=">")
+        assert track.iloc[0].to_list() == [-110.9536, -42.2489, -2823.96637605]
+    finally:
+        os.remove(path=TEMP_TRACK)
+
+    return output
+
+
 def test_grdtrack_input_dataframe_and_ncfile():
     """
     Run grdtrack by passing in a pandas.DataFrame and netcdf file as inputs
@@ -42,9 +66,29 @@ def test_grdtrack_input_dataframe_and_ncfile():
     return output
 
 
+def test_grdtrack_input_csvfile_and_ncfile():
+    """
+    Run grdtrack by passing in a csvfile and netcdf file as inputs
+    """
+    csvfile = which("@ridge.txt", download="c")
+    ncfile = which("@earth_relief_60m", download="c")
+
+    try:
+        output = grdtrack(points=csvfile, grid=ncfile, outfile=TEMP_TRACK)
+        assert output is None  # check that output is None since outfile is set
+        assert os.path.exists(path=TEMP_TRACK)  # check that outfile exists at path
+
+        track = pd.read_csv(TEMP_TRACK, sep="\t", header=None, comment=">")
+        assert track.iloc[0].to_list() == [-32.2971, 37.4118, -1697.87197487]
+    finally:
+        os.remove(path=TEMP_TRACK)
+
+    return output
+
+
 def test_grdtrack_wrong_kind_of_points_input():
     """
-    Run grdtrack using points input that is not a pandas.DataFrame (matrix)
+    Run grdtrack using points input that is not a pandas.DataFrame (matrix) or file
     """
     dataframe = load_ocean_ridge_points()
     invalid_points = dataframe.longitude.to_xarray()
@@ -77,3 +121,14 @@ def test_grdtrack_without_newcolname_setting():
 
     with pytest.raises(GMTInvalidInput):
         grdtrack(points=dataframe, grid=dataarray)
+
+
+def test_grdtrack_without_outfile_setting():
+    """
+    Run grdtrack by not passing in outfile parameter setting
+    """
+    csvfile = which("@ridge.txt", download="c")
+    dataarray = load_earth_relief().sel(lat=slice(-49, -42), lon=slice(-118, -107))
+
+    with pytest.raises(GMTInvalidInput):
+        grdtrack(points=csvfile, grid=dataarray)

From b1e7b75ab1b164830aa67e49d9e46147c2a7c01e Mon Sep 17 00:00:00 2001
From: Wei Ji <weiji.leong@vuw.ac.nz>
Date: Thu, 30 May 2019 21:48:23 +0200
Subject: [PATCH 07/11] Add gallery example 'sampling along tracks' for
 grdtrack

New raster grid gallery example for grdtrack! Plotting bathymetry sampled along ocean ridge points around the globe. As this 'grid' subsection is new to the gallery, I've added it to the ExplicitOrder part within the sphinx conf.py file.
---
 doc/conf.py                             |  6 +++-
 examples/gallery/grid/README.txt        |  2 ++
 examples/gallery/grid/track_sampling.py | 38 +++++++++++++++++++++++++
 3 files changed, 45 insertions(+), 1 deletion(-)
 create mode 100644 examples/gallery/grid/README.txt
 create mode 100644 examples/gallery/grid/track_sampling.py

diff --git a/doc/conf.py b/doc/conf.py
index 3d2876c03ed..787e4dfec17 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -56,7 +56,11 @@
     # path where to save gallery generated examples
     "gallery_dirs": ["gallery", "tutorials", "projections"],
     "subsection_order": ExplicitOrder(
-        ["../examples/gallery/coast", "../examples/gallery/plot"]
+        [
+            "../examples/gallery/coast",
+            "../examples/gallery/plot",
+            "../examples/gallery/grid",
+        ]
     ),
     # Patter to search for example files
     "filename_pattern": r"\.py",
diff --git a/examples/gallery/grid/README.txt b/examples/gallery/grid/README.txt
new file mode 100644
index 00000000000..912195c11b8
--- /dev/null
+++ b/examples/gallery/grid/README.txt
@@ -0,0 +1,2 @@
+Raster grids
+------------
diff --git a/examples/gallery/grid/track_sampling.py b/examples/gallery/grid/track_sampling.py
new file mode 100644
index 00000000000..3fc6603624a
--- /dev/null
+++ b/examples/gallery/grid/track_sampling.py
@@ -0,0 +1,38 @@
+"""
+Sampling along tracks
+---------------------
+
+The :meth:`pygmt.grdtrack` method samples a raster grid's value along specified points.
+We will need to input a 2D raster to ``grid`` which can be an ``xarray.DataArray``. The
+``points`` parameter can be a ``pandas.DataFrame`` table where the first two columns are
+x and y (or longitude and latitude). Note also that there is a ``newcolname`` parameter
+to set that will be used to name the new column of values we sampled from the grid.
+
+Alternatively, we can provide a NetCDF file path to ``grid``. An ASCII file path can
+also be accepted for ``points``, but an ``outfile`` parameter will then need to be set
+to name the resulting output ASCII file.
+"""
+
+import pygmt
+
+# Load sample grid and point datasets
+grid = pygmt.datasets.load_earth_relief()
+points = pygmt.datasets.load_ocean_ridge_points()
+# Sample the bathymetry along the world's ocean ridges at specified track points
+track = pygmt.grdtrack(points=points, grid=grid, newcolname="bathymetry")
+
+fig = pygmt.Figure()
+# Plot the earth relief grid on Cylindrical Stereographic projection, masking land areas
+fig.basemap(region="g", frame=True, projection="Cyl_stere/30/-20/8i")
+fig.grdimage(grid=grid, cmap="ocean")
+fig.coast(land="#666666")
+# Plot using circles (c) of 0.1cm, the sampled bathymetry points
+# Points are colored using elevation values (normalized for visual purposes)
+fig.plot(
+    x=track.longitude,
+    y=track.latitude,
+    style="c0.1c",
+    cmap="terra",
+    color=(track.bathymetry - track.bathymetry.mean()) / track.bathymetry.std(),
+)
+fig.show()

From 1f948965d8961df7b311376fb5ad55fd3ca0e6e5 Mon Sep 17 00:00:00 2001
From: Wei Ji <weiji.leong@vuw.ac.nz>
Date: Sun, 9 Jun 2019 12:38:02 +0200
Subject: [PATCH 08/11] Finalize grdtrack wrapper, use numpy assertion, change
 grid bg to gray

Using numpy.testing.assert_allclose in test_grdtrack.py unit tests and make some cosmetic changes to docs and gallery example. To make the gallery example plot clearer, grid background colormap was changed from 'ocean' to 'gray', and the points are made slightly bigger. Also changed map to be Pacific-centered!
---
 examples/gallery/grid/README.txt        |  2 +-
 examples/gallery/grid/track_sampling.py | 19 ++++++++++---------
 pygmt/tests/test_grdtrack.py            |  9 +++++----
 3 files changed, 16 insertions(+), 14 deletions(-)

diff --git a/examples/gallery/grid/README.txt b/examples/gallery/grid/README.txt
index 912195c11b8..23f2062152e 100644
--- a/examples/gallery/grid/README.txt
+++ b/examples/gallery/grid/README.txt
@@ -1,2 +1,2 @@
-Raster grids
+Grids
 ------------
diff --git a/examples/gallery/grid/track_sampling.py b/examples/gallery/grid/track_sampling.py
index 3fc6603624a..66bd9cb2316 100644
--- a/examples/gallery/grid/track_sampling.py
+++ b/examples/gallery/grid/track_sampling.py
@@ -2,11 +2,12 @@
 Sampling along tracks
 ---------------------
 
-The :meth:`pygmt.grdtrack` method samples a raster grid's value along specified points.
-We will need to input a 2D raster to ``grid`` which can be an ``xarray.DataArray``. The
-``points`` parameter can be a ``pandas.DataFrame`` table where the first two columns are
-x and y (or longitude and latitude). Note also that there is a ``newcolname`` parameter
-to set that will be used to name the new column of values we sampled from the grid.
+The :func:`pygmt.grdtrack` function samples a raster grid's value along specified
+points. We will need to input a 2D raster to ``grid`` which can be an
+``xarray.DataArray``. The ``points`` parameter can be a ``pandas.DataFrame`` table where
+the first two columns are x and y (or longitude and latitude). Note also that there is a
+``newcolname`` parameter that will be used to name the new column of values we sampled
+from the grid.
 
 Alternatively, we can provide a NetCDF file path to ``grid``. An ASCII file path can
 also be accepted for ``points``, but an ``outfile`` parameter will then need to be set
@@ -23,15 +24,15 @@
 
 fig = pygmt.Figure()
 # Plot the earth relief grid on Cylindrical Stereographic projection, masking land areas
-fig.basemap(region="g", frame=True, projection="Cyl_stere/30/-20/8i")
-fig.grdimage(grid=grid, cmap="ocean")
+fig.basemap(region="g", frame=True, projection="Cyl_stere/150/-20/8i")
+fig.grdimage(grid=grid, cmap="gray")
 fig.coast(land="#666666")
-# Plot using circles (c) of 0.1cm, the sampled bathymetry points
+# Plot using circles (c) of 0.15cm, the sampled bathymetry points
 # Points are colored using elevation values (normalized for visual purposes)
 fig.plot(
     x=track.longitude,
     y=track.latitude,
-    style="c0.1c",
+    style="c0.15c",
     cmap="terra",
     color=(track.bathymetry - track.bathymetry.mean()) / track.bathymetry.std(),
 )
diff --git a/pygmt/tests/test_grdtrack.py b/pygmt/tests/test_grdtrack.py
index 4801ee0428f..89359e8a507 100644
--- a/pygmt/tests/test_grdtrack.py
+++ b/pygmt/tests/test_grdtrack.py
@@ -3,6 +3,7 @@
 """
 import os
 
+import numpy.testing as npt
 import pandas as pd
 import pytest
 
@@ -26,7 +27,7 @@ def test_grdtrack_input_dataframe_and_dataarray():
     output = grdtrack(points=dataframe, grid=dataarray, newcolname="bathymetry")
     assert isinstance(output, pd.DataFrame)
     assert output.columns.to_list() == ["longitude", "latitude", "bathymetry"]
-    assert output.iloc[0].to_list() == [-110.9536, -42.2489, -2823.96637605]
+    npt.assert_allclose(output.iloc[0], [-110.9536, -42.2489, -2823.96637605])
 
     return output
 
@@ -44,7 +45,7 @@ def test_grdtrack_input_csvfile_and_dataarray():
         assert os.path.exists(path=TEMP_TRACK)  # check that outfile exists at path
 
         track = pd.read_csv(TEMP_TRACK, sep="\t", header=None, comment=">")
-        assert track.iloc[0].to_list() == [-110.9536, -42.2489, -2823.96637605]
+        npt.assert_allclose(track.iloc[0], [-110.9536, -42.2489, -2823.96637605])
     finally:
         os.remove(path=TEMP_TRACK)
 
@@ -61,7 +62,7 @@ def test_grdtrack_input_dataframe_and_ncfile():
     output = grdtrack(points=dataframe, grid=ncfile, newcolname="bathymetry")
     assert isinstance(output, pd.DataFrame)
     assert output.columns.to_list() == ["longitude", "latitude", "bathymetry"]
-    assert output.iloc[0].to_list() == [-32.2971, 37.4118, -1697.87197487]
+    npt.assert_allclose(output.iloc[0], [-32.2971, 37.4118, -1697.87197487])
 
     return output
 
@@ -79,7 +80,7 @@ def test_grdtrack_input_csvfile_and_ncfile():
         assert os.path.exists(path=TEMP_TRACK)  # check that outfile exists at path
 
         track = pd.read_csv(TEMP_TRACK, sep="\t", header=None, comment=">")
-        assert track.iloc[0].to_list() == [-32.2971, 37.4118, -1697.87197487]
+        npt.assert_allclose(track.iloc[0], [-32.2971, 37.4118, -1697.87197487])
     finally:
         os.remove(path=TEMP_TRACK)
 

From afb43828a669905658a7156207bd5c17e60376b9 Mon Sep 17 00:00:00 2001
From: Wei Ji <weiji.leong@vuw.ac.nz>
Date: Fri, 1 Nov 2019 17:45:42 +1300
Subject: [PATCH 09/11] Fix bathymetry values in test_grdtrack

Due to new SRTM15+V2 grids, see also #350.
---
 pygmt/tests/test_grdtrack.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/pygmt/tests/test_grdtrack.py b/pygmt/tests/test_grdtrack.py
index 89359e8a507..565343a5de3 100644
--- a/pygmt/tests/test_grdtrack.py
+++ b/pygmt/tests/test_grdtrack.py
@@ -27,7 +27,7 @@ def test_grdtrack_input_dataframe_and_dataarray():
     output = grdtrack(points=dataframe, grid=dataarray, newcolname="bathymetry")
     assert isinstance(output, pd.DataFrame)
     assert output.columns.to_list() == ["longitude", "latitude", "bathymetry"]
-    npt.assert_allclose(output.iloc[0], [-110.9536, -42.2489, -2823.96637605])
+    npt.assert_allclose(output.iloc[0], [-110.9536, -42.2489, -2799.085897])
 
     return output
 
@@ -45,7 +45,7 @@ def test_grdtrack_input_csvfile_and_dataarray():
         assert os.path.exists(path=TEMP_TRACK)  # check that outfile exists at path
 
         track = pd.read_csv(TEMP_TRACK, sep="\t", header=None, comment=">")
-        npt.assert_allclose(track.iloc[0], [-110.9536, -42.2489, -2823.96637605])
+        npt.assert_allclose(track.iloc[0], [-110.9536, -42.2489, -2799.085897])
     finally:
         os.remove(path=TEMP_TRACK)
 
@@ -62,7 +62,7 @@ def test_grdtrack_input_dataframe_and_ncfile():
     output = grdtrack(points=dataframe, grid=ncfile, newcolname="bathymetry")
     assert isinstance(output, pd.DataFrame)
     assert output.columns.to_list() == ["longitude", "latitude", "bathymetry"]
-    npt.assert_allclose(output.iloc[0], [-32.2971, 37.4118, -1697.87197487])
+    npt.assert_allclose(output.iloc[0], [-32.2971, 37.4118, -1685.745884])
 
     return output
 
@@ -80,7 +80,7 @@ def test_grdtrack_input_csvfile_and_ncfile():
         assert os.path.exists(path=TEMP_TRACK)  # check that outfile exists at path
 
         track = pd.read_csv(TEMP_TRACK, sep="\t", header=None, comment=">")
-        npt.assert_allclose(track.iloc[0], [-32.2971, 37.4118, -1697.87197487])
+        npt.assert_allclose(track.iloc[0], [-32.2971, 37.4118, -1685.745884])
     finally:
         os.remove(path=TEMP_TRACK)
 

From 0ed8eb4bd8e6541fc76a92da1b0f94f21c839acb Mon Sep 17 00:00:00 2001
From: Wei Ji <weiji.leong@vuw.ac.nz>
Date: Wed, 25 Mar 2020 10:40:26 +1300
Subject: [PATCH 10/11] Fix test values and truncate docstrings to 79 char

Make tests pass again, by fixing the test bathymetry values, and cutting our docstrings to 79 characters instead of 88.
---
 pygmt/datasets/tutorial.py   |  3 ++-
 pygmt/sampling.py            | 29 +++++++++++++++++------------
 pygmt/tests/test_grdtrack.py | 17 ++++++++++-------
 3 files changed, 29 insertions(+), 20 deletions(-)

diff --git a/pygmt/datasets/tutorial.py b/pygmt/datasets/tutorial.py
index 9263d22e0b5..b56e0ce0e1e 100644
--- a/pygmt/datasets/tutorial.py
+++ b/pygmt/datasets/tutorial.py
@@ -40,7 +40,8 @@ def load_japan_quakes():
 
 def load_ocean_ridge_points():
     """
-    Load a table of ocean ridge points for the entire world as a pandas.DataFrame.
+    Load a table of ocean ridge points for the entire world as a
+    pandas.DataFrame.
 
     This is the ``@ridge.txt`` dataset used in the GMT tutorials.
 
diff --git a/pygmt/sampling.py b/pygmt/sampling.py
index 7df01eb2884..4ed2a6d171b 100644
--- a/pygmt/sampling.py
+++ b/pygmt/sampling.py
@@ -20,35 +20,39 @@ def grdtrack(points, grid, newcolname=None, outfile=None, **kwargs):
     Sample grids at specified (x,y) locations.
 
     Grdtrack reads one or more grid files and a table with (x,y) [or (lon,lat)]
-    positions in the first two columns (more columns may be present). It interpolates
-    the grid(s) at the positions in the table and writes out the table with the
-    interpolated values added as (one or more) new columns. A bicubic [Default],
-    bilinear, B-spline or nearest-neighbor (see -n) interpolation is used, requiring
-    boundary conditions at the limits of the region.
+    positions in the first two columns (more columns may be present). It
+    interpolates the grid(s) at the positions in the table and writes out the
+    table with the interpolated values added as (one or more) new columns. A
+    bicubic [Default], bilinear, B-spline or nearest-neighbor (see -n)
+    interpolation is used, requiring boundary conditions at the limits of the
+    region.
 
     Full option list at :gmt-docs:`grdtrack.html`
 
     Parameters
     ----------
     points: pandas.DataFrame or file (csv, txt, etc)
-        Either a table with (x, y) or (lon, lat) values in the first two columns,
-        or a data file name. More columns may be present.
+        Either a table with (x, y) or (lon, lat) values in the first two
+        columns, or a data file name. More columns may be present.
 
     grid: xarray.DataArray or file (netcdf)
         Gridded array from which to sample values from.
 
     newcolname: str
-        Required if 'points' is a pandas.DataFrame. The name for the new column in the
-        track pandas.DataFrame table where the sampled values will be placed.
+        Required if 'points' is a pandas.DataFrame. The name for the new column
+        in the track pandas.DataFrame table where the sampled values will be
+        placed.
 
     outfile: str
-        Required if 'points' is a file. The file name for the output ASCII file.
+        Required if 'points' is a file. The file name for the output ASCII
+        file.
 
     Returns
     -------
     track: pandas.DataFrame or None
         Return type depends on whether the outfile parameter is set:
-        - pandas.DataFrame table with (x, y, ..., newcolname) if outfile is not set
+        - pandas.DataFrame table with (x, y, ..., newcolname) if outfile is not
+          set
         - None if outfile is set (track output will be stored in outfile)
 
     """
@@ -75,7 +79,8 @@ def grdtrack(points, grid, newcolname=None, outfile=None, **kwargs):
             else:
                 raise GMTInvalidInput(f"Unrecognized data type {type(grid)}")
 
-            # Run grdtrack on the temp (csv) points table and (netcdf) grid virtualfiles
+            # Run grdtrack on the temporary (csv) points table
+            # and (netcdf) grid virtualfile
             with table_context as csvfile:
                 with grid_context as grdfile:
                     kwargs.update({"G": grdfile})
diff --git a/pygmt/tests/test_grdtrack.py b/pygmt/tests/test_grdtrack.py
index 565343a5de3..acc1c5c00c7 100644
--- a/pygmt/tests/test_grdtrack.py
+++ b/pygmt/tests/test_grdtrack.py
@@ -19,7 +19,8 @@
 
 def test_grdtrack_input_dataframe_and_dataarray():
     """
-    Run grdtrack by passing in a pandas.DataFrame and xarray.DataArray as inputs
+    Run grdtrack by passing in a pandas.DataFrame and xarray.DataArray as
+    inputs
     """
     dataframe = load_ocean_ridge_points()
     dataarray = load_earth_relief().sel(lat=slice(-49, -42), lon=slice(-118, -107))
@@ -27,7 +28,7 @@ def test_grdtrack_input_dataframe_and_dataarray():
     output = grdtrack(points=dataframe, grid=dataarray, newcolname="bathymetry")
     assert isinstance(output, pd.DataFrame)
     assert output.columns.to_list() == ["longitude", "latitude", "bathymetry"]
-    npt.assert_allclose(output.iloc[0], [-110.9536, -42.2489, -2799.085897])
+    npt.assert_allclose(output.iloc[0], [-110.9536, -42.2489, -2797.482959])
 
     return output
 
@@ -45,7 +46,7 @@ def test_grdtrack_input_csvfile_and_dataarray():
         assert os.path.exists(path=TEMP_TRACK)  # check that outfile exists at path
 
         track = pd.read_csv(TEMP_TRACK, sep="\t", header=None, comment=">")
-        npt.assert_allclose(track.iloc[0], [-110.9536, -42.2489, -2799.085897])
+        npt.assert_allclose(track.iloc[0], [-110.9536, -42.2489, -2797.482959])
     finally:
         os.remove(path=TEMP_TRACK)
 
@@ -62,7 +63,7 @@ def test_grdtrack_input_dataframe_and_ncfile():
     output = grdtrack(points=dataframe, grid=ncfile, newcolname="bathymetry")
     assert isinstance(output, pd.DataFrame)
     assert output.columns.to_list() == ["longitude", "latitude", "bathymetry"]
-    npt.assert_allclose(output.iloc[0], [-32.2971, 37.4118, -1685.745884])
+    npt.assert_allclose(output.iloc[0], [-32.2971, 37.4118, -1686.692142])
 
     return output
 
@@ -80,7 +81,7 @@ def test_grdtrack_input_csvfile_and_ncfile():
         assert os.path.exists(path=TEMP_TRACK)  # check that outfile exists at path
 
         track = pd.read_csv(TEMP_TRACK, sep="\t", header=None, comment=">")
-        npt.assert_allclose(track.iloc[0], [-32.2971, 37.4118, -1685.745884])
+        npt.assert_allclose(track.iloc[0], [-32.2971, 37.4118, -1686.692142])
     finally:
         os.remove(path=TEMP_TRACK)
 
@@ -89,7 +90,8 @@ def test_grdtrack_input_csvfile_and_ncfile():
 
 def test_grdtrack_wrong_kind_of_points_input():
     """
-    Run grdtrack using points input that is not a pandas.DataFrame (matrix) or file
+    Run grdtrack using points input that is not a pandas.DataFrame (matrix) or
+    file
     """
     dataframe = load_ocean_ridge_points()
     invalid_points = dataframe.longitude.to_xarray()
@@ -102,7 +104,8 @@ def test_grdtrack_wrong_kind_of_points_input():
 
 def test_grdtrack_wrong_kind_of_grid_input():
     """
-    Run grdtrack using grid input that is not as xarray.DataArray (grid) or file
+    Run grdtrack using grid input that is not as xarray.DataArray (grid) or
+    file
     """
     dataframe = load_ocean_ridge_points()
     dataarray = load_earth_relief().sel(lat=slice(-49, -42), lon=slice(-118, -107))

From f0b6de3b7f725e5c2fd34a300703b3468e99b54f Mon Sep 17 00:00:00 2001
From: Wei Ji <weiji.leong@vuw.ac.nz>
Date: Thu, 26 Mar 2020 10:50:08 +1300
Subject: [PATCH 11/11] Minor formatting fix

---
 examples/gallery/grid/README.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/gallery/grid/README.txt b/examples/gallery/grid/README.txt
index 23f2062152e..6fd2c9aea1a 100644
--- a/examples/gallery/grid/README.txt
+++ b/examples/gallery/grid/README.txt
@@ -1,2 +1,2 @@
 Grids
-------------
+-----