Add function to import seafloor crustal age dataset (#1471)

New function to download the Earth seafloor age datasets
at various resolutions. This is the `@earth_age` grids on the
GMT remote data server.

* add earth_age.py
* Add test_datasets_earth_age.py
* add seafloor age datasets to helpers/testing.py
* add load_earth_age to index.rst
* add 05m grid tile for caching

Co-authored-by: Wei Ji <[email protected]>

doc/api/index.rst

@@ -170,6 +170,7 @@ and store them in the GMT cache folder.
 .. autosummary::
  :toctree: generated
 
+ datasets.load_earth_age
  datasets.load_earth_relief
  datasets.load_fractures_compilation
  datasets.load_hotspots

pygmt/datasets/__init__.py

@@ -2,6 +2,7 @@
 #
 # Load sample data included with GMT (downloaded from the GMT cache server).
 
+from pygmt.datasets.earth_age import load_earth_age
 from pygmt.datasets.earth_relief import load_earth_relief
 from pygmt.datasets.samples import (
  load_fractures_compilation,

pygmt/datasets/earth_age.py

@@ -0,0 +1,113 @@
+"""
+Function to download the Earth seafloor age datasets from the GMT data server,
+and load as :class:`xarray.DataArray`.
+
+The grids are available in various resolutions.
+"""
+from pygmt.exceptions import GMTInvalidInput
+from pygmt.helpers import kwargs_to_strings
+from pygmt.io import load_dataarray
+from pygmt.src import grdcut, which
+
+
+@kwargs_to_strings(region="sequence")
+def load_earth_age(resolution="01d", region=None, registration=None):
+ r"""
+ Load Earth seafloor crustal ages in various resolutions.
+
+ The grids are downloaded to a user data directory
+ (usually ``~/.gmt/server/earth/earth_age/``) the first time you invoke
+ this function. Afterwards, it will load the grid from the data directory.
+ So you'll need an internet connection the first time around.
+
+ These grids can also be accessed by passing in the file name
+ **@earth_age**\_\ *res*\[_\ *reg*] to any grid plotting/processing
+ function. *res* is the grid resolution (see below), and *reg* is grid
+ registration type (**p** for pixel registration or **g** for gridline
+ registration).
+
+ Refer to
+ :gmt-docs:`datasets/remote-data.html#global-earth-seafloor-crustal-age-grids`
+ for more details.
+
+ Parameters
+ ----------
+ resolution : str
+ The grid resolution. The suffix ``d`` and ``m`` stand for
+ arc-degree, arc-minute and arc-second. It can be ``'01d'``, ``'30m'``,
+ ``'20m'``, ``'15m'``, ``'10m'``, ``'06m'``, ``'05m'``, ``'04m'``,
+ ``'03m'``, ``'02m'``, or ``'01m'``.
+
+ region : str or list
+ The subregion of the grid to load, in the forms of a list
+ [*xmin*, *xmax*, *ymin*, *ymax*] or a string *xmin/xmax/ymin/ymax*.
+ Required for grids with resolutions higher than 5
+ arc-minute (i.e., ``05m``).
+
+ registration : str
+ Grid registration type. Either ``pixel`` for pixel registration or
+ ``gridline`` for gridline registration. Default is ``None``, where
+ a pixel-registered grid is returned unless only the
+ gridline-registered grid is available.
+
+ Returns
+ -------
+ grid : :class:`xarray.DataArray`
+ The Earth seafloor crustal age grid. Coordinates are latitude and
+ longitude in degrees. Age is in millions of years (Myr).
+
+ Notes
+ -----
+ The :class:`xarray.DataArray` grid doesn't support slice operation, for
+ Earth seafloor crustal age with resolutions of 5 arc-minutes or higher,
+ which are stored as smaller tiles.
+ """
+
+ # earth seafloor crust age data stored as single grids for low resolutions
+ non_tiled_resolutions = ["01d", "30m", "20m", "15m", "10m", "06m"]
+ # earth seafloor crust age data stored as tiles for high resolutions
+ tiled_resolutions = ["05m", "04m", "03m", "02m", "01m"]
+
+ if registration in ("pixel", "gridline", None):
+ # If None, let GMT decide on Pixel/Gridline type
+ reg = f"_{registration[0]}" if registration else ""
+ else:
+ raise GMTInvalidInput(
+ f"Invalid grid registration: '{registration}', should be either "
+ "'pixel', 'gridline' or None. Default is None, where a "
+ "pixel-registered grid is returned unless only the "
+ "gridline-registered grid is available."
+ )
+
+ if resolution not in non_tiled_resolutions + tiled_resolutions:
+ raise GMTInvalidInput(f"Invalid Earth relief resolution '{resolution}'.")
+
+ # Choose earth relief data prefix
+ earth_age_prefix = "earth_age_"
+
+ # different ways to load tiled and non-tiled earth relief data
+ # Known issue: tiled grids don't support slice operation
+ # See https://github.com/GenericMappingTools/pygmt/issues/524
+ if region is None:
+ if resolution not in non_tiled_resolutions:
+ raise GMTInvalidInput(
+ f"'region' is required for Earth age resolution '{resolution}'."
+ )
+ fname = which(f"@earth_age_{resolution}{reg}", download="a")
+ grid = load_dataarray(fname, engine="netcdf4")
+ else:
+ grid = grdcut(f"@{earth_age_prefix}{resolution}{reg}", region=region)
+
+ # Add some metadata to the grid
+ grid.name = "seafloor_age"
+ grid.attrs["long_name"] = "age of seafloor crust"
+ grid.attrs["units"] = "Myr"
+ grid.attrs["vertical_datum"] = "EMG96"
+ grid.attrs["horizontal_datum"] = "WGS84"
+ # Remove the actual range because it gets outdated when indexing the grid,
+ # which causes problems when exporting it to netCDF for usage on the
+ # command-line.
+ grid.attrs.pop("actual_range")
+ for coord in grid.coords:
+ grid[coord].attrs.pop("actual_range")
+ return grid

pygmt/helpers/testing.py

@@ -165,6 +165,9 @@ def download_test_data():
  "@N35E135.earth_relief_03s_g.nc",
  "@N37W120.earth_relief_03s_g.nc",
  "@N00W090.earth_relief_03m_p.nc",
+ # Earth seafloor age grids
+ "@earth_age_01d_g",
+ "S90W180.earth_age_05m_g.nc" # Specific grid for 05m test
  # Other cache files
  "@EGM96_to_36.txt",
  "@fractures_06.txt",

pygmt/tests/test_datasets_earth_age.py

@@ -0,0 +1,78 @@
+"""
+Test basic functionality for loading Earth seafloor crust age datasets.
+"""
+import numpy as np
+import numpy.testing as npt
+import pytest
+from pygmt.datasets import load_earth_age
+from pygmt.exceptions import GMTInvalidInput
+
+
+def test_earth_age_fails():
+ """
+ Make sure earth_age fails for invalid resolutions.
+ """
+ resolutions = "1m 1d bla 60d 001m 03".split()
+ resolutions.append(60)
+ for resolution in resolutions:
+ with pytest.raises(GMTInvalidInput):
+ load_earth_age(resolution=resolution)
+
+
+def test_earth_age_incorrect_registration():
+ """
+ Test loading earth_age with incorrect registration type.
+ """
+ with pytest.raises(GMTInvalidInput):
+ load_earth_age(registration="improper_type")
+
+
+def test_earth_age_01d():
+ """
+ Test some properties of the earth age 01d data.
+ """
+ data = load_earth_age(resolution="01d", registration="gridline")
+ assert data.shape == (181, 361)
+ npt.assert_allclose(data.lat, np.arange(-90, 91, 1))
+ npt.assert_allclose(data.lon, np.arange(-180, 181, 1))
+ npt.assert_allclose(data.min(), 0.167381, rtol=1e-5)
+ npt.assert_allclose(data.max(), 338.0274, rtol=1e-5)
+
+
+def test_earth_age_01d_with_region():
+ """
+ Test loading low-resolution earth age with 'region'.
+ """
+ data = load_earth_age(
+ resolution="01d", region=[-10, 10, -5, 5], registration="gridline"
+ )
+ assert data.shape == (11, 21)
+ npt.assert_allclose(data.lat, np.arange(-5, 6, 1))
+ npt.assert_allclose(data.lon, np.arange(-10, 11, 1))
+ npt.assert_allclose(data.min(), 11.293945)
+ npt.assert_allclose(data.max(), 125.1189)
+
+
+def test_earth_age_05m_with_region():
+ """
+ Test loading a subregion of high-resolution earth age.
+ """
+ data = load_earth_age(
+ resolution="05m", region=[-50, -40, 20, 30], registration="gridline"
+ )
+ assert data.coords["lat"].data.min() == 20.0
+ assert data.coords["lat"].data.max() == 30.0
+ assert data.coords["lon"].data.min() == -50.0
+ assert data.coords["lon"].data.max() == -40.0
+ npt.assert_allclose(data.data.min(), 0.040000916)
+ npt.assert_allclose(data.data.max(), 46.530003)
+ assert data.sizes["lat"] == 121
+ assert data.sizes["lon"] == 121
+
+
+def test_earth_age_05m_without_region():
+ """
+ Test loading high-resolution earth age without passing 'region'.
+ """
+ with pytest.raises(GMTInvalidInput):
+ load_earth_age("05m")