From 6083975ccd7319ec2f42fcb9a424eb6dff2ab913 Mon Sep 17 00:00:00 2001 From: Will Schlitzer Date: Thu, 22 Dec 2022 22:17:44 -0500 Subject: [PATCH] Add load_earth_vertical_gravity_gradient function for Earth vertical gravity gradient dataset (#2240) Co-authored-by: Dongdong Tian Co-authored-by: Michael Grund <23025878+michaelgrund@users.noreply.github.com> --- doc/api/index.rst | 1 + pygmt/datasets/__init__.py | 3 + .../earth_vertical_gravity_gradient.py | 71 ++++++++++++++ pygmt/datasets/load_remote_dataset.py | 20 ++++ pygmt/helpers/testing.py | 3 + ...atasets_earth_vertical_gravity_gradient.py | 98 +++++++++++++++++++ 6 files changed, 196 insertions(+) create mode 100644 pygmt/datasets/earth_vertical_gravity_gradient.py create mode 100644 pygmt/tests/test_datasets_earth_vertical_gravity_gradient.py diff --git a/doc/api/index.rst b/doc/api/index.rst index 21fec1f16bd..618e95b6786 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -226,6 +226,7 @@ and store them in GMT's user data directory. datasets.load_earth_geoid datasets.load_earth_magnetic_anomaly datasets.load_earth_relief + datasets.load_earth_vertical_gravity_gradient datasets.load_sample_data The following functions are deprecated since v0.6.0 and will be removed in v0.9.0. diff --git a/pygmt/datasets/__init__.py b/pygmt/datasets/__init__.py index 735be59d83b..87595e0b6bc 100644 --- a/pygmt/datasets/__init__.py +++ b/pygmt/datasets/__init__.py @@ -7,6 +7,9 @@ from pygmt.datasets.earth_geoid import load_earth_geoid from pygmt.datasets.earth_magnetic_anomaly import load_earth_magnetic_anomaly from pygmt.datasets.earth_relief import load_earth_relief +from pygmt.datasets.earth_vertical_gravity_gradient import ( + load_earth_vertical_gravity_gradient, +) from pygmt.datasets.samples import ( list_sample_data, load_fractures_compilation, diff --git a/pygmt/datasets/earth_vertical_gravity_gradient.py b/pygmt/datasets/earth_vertical_gravity_gradient.py new file mode 100644 index 00000000000..b47daad94c7 --- /dev/null +++ b/pygmt/datasets/earth_vertical_gravity_gradient.py @@ -0,0 +1,71 @@ +""" +Function to download the IGPP Global Earth Vertical Gravity Gradient from the +GMT data server, and load as :class:`xarray.DataArray`. + +The grids are available in various resolutions. +""" +from pygmt.datasets.load_remote_dataset import _load_remote_dataset +from pygmt.helpers import kwargs_to_strings + + +@kwargs_to_strings(region="sequence") +def load_earth_vertical_gravity_gradient( + resolution="01d", region=None, registration=None +): + r""" + Load the IGPP Global Earth Vertical Gravity Gradient in various + resolutions. + + The grids are downloaded to a user data directory + (usually ``~/.gmt/server/earth/earth_vgg/``) 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_vgg**\_\ *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-datasets:`earth-vgg.html` for more details. + + Parameters + ---------- + resolution : str + The grid resolution. The suffix ``d`` and ``m`` stand for + arc-degrees and arc-minutes. 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 form of a list + [*xmin*, *xmax*, *ymin*, *ymax*] or a string *xmin/xmax/ymin/ymax*. + Required for grids with resolutions higher than 5 + arc-minutes (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 vertical gravity gradient grid. Coordinates are latitude and + longitude in degrees. Units are in Eotvos. + + Note + ---- + The :class:`xarray.DataArray` grid doesn't support slice operation, for + Earth vertical gravity gradient grids with resolutions of 5 arc-minutes or + higher, which are stored as smaller tiles. + """ + grid = _load_remote_dataset( + dataset_name="earth_vgg", + dataset_prefix="earth_vgg_", + resolution=resolution, + region=region, + registration=registration, + ) + return grid diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py index 81614b40fe2..915f31b15bd 100644 --- a/pygmt/datasets/load_remote_dataset.py +++ b/pygmt/datasets/load_remote_dataset.py @@ -167,6 +167,26 @@ class GMTRemoteDataset(NamedTuple): "01m": Resolution(["pixel"], True), }, ), + "earth_vgg": GMTRemoteDataset( + title="Earth vertical gravity gradient", + name="earth_vgg", + long_name="IGPP Global Earth Vertical Gravity Gradient", + units="Eotvos", + extra_attributes={"horizontal_datum": "WGS84"}, + resolutions={ + "01d": Resolution(["pixel", "gridline"], False), + "30m": Resolution(["pixel", "gridline"], False), + "20m": Resolution(["pixel", "gridline"], False), + "15m": Resolution(["pixel", "gridline"], False), + "10m": Resolution(["pixel", "gridline"], False), + "06m": Resolution(["pixel", "gridline"], False), + "05m": Resolution(["pixel", "gridline"], True), + "04m": Resolution(["pixel", "gridline"], True), + "03m": Resolution(["pixel", "gridline"], True), + "02m": Resolution(["pixel", "gridline"], True), + "01m": Resolution(["pixel"], True), + }, + ), } diff --git a/pygmt/helpers/testing.py b/pygmt/helpers/testing.py index 86f974571ea..06eb3af0482 100644 --- a/pygmt/helpers/testing.py +++ b/pygmt/helpers/testing.py @@ -191,6 +191,9 @@ def download_test_data(): # Earth free-air anomaly grids "@earth_faa_01d_g", "@S90W180.earth_faa_05m_g.nc", # Specific grid for 05m test + # Earth vertical gravity gradient grids + "@earth_vgg_01d_g", + "@S90W180.earth_vgg_05m_g.nc", # Specific grid for 05m test # Other cache files "@capitals.gmt", "@earth_relief_20m_holes.grd", diff --git a/pygmt/tests/test_datasets_earth_vertical_gravity_gradient.py b/pygmt/tests/test_datasets_earth_vertical_gravity_gradient.py new file mode 100644 index 00000000000..2bbb85764d5 --- /dev/null +++ b/pygmt/tests/test_datasets_earth_vertical_gravity_gradient.py @@ -0,0 +1,98 @@ +""" +Test basic functionality for loading Earth vertical gravity gradient datasets. +""" +import numpy as np +import numpy.testing as npt +import pytest +from pygmt.datasets import load_earth_vertical_gravity_gradient +from pygmt.exceptions import GMTInvalidInput + + +def test_earth_vertical_gravity_gradient_fails(): + """ + Make sure load_earth_vertical_gravity_gradient 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_vertical_gravity_gradient(resolution=resolution) + + +def test_earth_vertical_gravity_gradient_incorrect_registration(): + """ + Test loading load_earth_vertical_gravity_gradient with incorrect + registration type. + """ + with pytest.raises(GMTInvalidInput): + load_earth_vertical_gravity_gradient(registration="improper_type") + + +def test_earth_vertical_gravity_gradient_01d(): + """ + Test some properties of the earth vgg 01d data. + """ + data = load_earth_vertical_gravity_gradient( + resolution="01d", registration="gridline" + ) + assert data.name == "earth_vgg" + assert data.attrs["units"] == "Eotvos" + assert data.attrs["long_name"] == "IGPP Global Earth Vertical Gravity Gradient" + assert data.attrs["horizontal_datum"] == "WGS84" + 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(), -136.34375) + npt.assert_allclose(data.max(), 104.59375) + assert data[1, 1].isnull() + + +def test_earth_vertical_gravity_gradient_01d_with_region(): + """ + Test loading low-resolution earth vgg with 'region'. + """ + data = load_earth_vertical_gravity_gradient( + 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(), -16.34375) + npt.assert_allclose(data.max(), 19.78125) + + +def test_earth_vertical_gravity_gradient_05m_with_region(): + """ + Test loading a subregion of high-resolution earth vgg. + """ + data = load_earth_vertical_gravity_gradient( + resolution="05m", region=[-50, -40, 20, 26], registration="gridline" + ) + assert data.coords["lat"].data.min() == 20.0 + assert data.coords["lat"].data.max() == 26.0 + assert data.coords["lon"].data.min() == -50.0 + assert data.coords["lon"].data.max() == -40.0 + npt.assert_allclose(data.min(), -107.625) + npt.assert_allclose(data.max(), 159.75) + assert data.sizes["lat"] == 73 + assert data.sizes["lon"] == 121 + + +def test_earth_vertical_gravity_gradient_05m_without_region(): + """ + Test loading high-resolution earth vgg without passing 'region'. + """ + with pytest.raises(GMTInvalidInput): + load_earth_vertical_gravity_gradient("05m") + + +def test_earth_vertical_gravity_gradient_incorrect_resolution_registration(): + """ + Test that an error is raised when trying to load a grid registration with + an unavailable resolution. + """ + with pytest.raises(GMTInvalidInput): + load_earth_vertical_gravity_gradient( + resolution="01m", region=[0, 1, 3, 5], registration="gridline" + )