From 923111a8ab05b36e354b1b01be95e5621ff0fc20 Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Mon, 30 Sep 2024 17:50:33 +1300 Subject: [PATCH] Add function to load Black Marble dataset (#3469) Function to load the `@earth_night_` NASA Black Marble mosaic RGB images stored as GeoTIFF files. * Add entry for datasets.load_black_marble in doc/api/index.rst * Add unit tests for load_black_marble * Import rioxarray to register rio accessor * Add @earth_day_01d to pygmt/helpers/caching.py --- doc/api/index.rst | 1 + pygmt/datasets/__init__.py | 1 + pygmt/datasets/earth_day.py | 5 ++ pygmt/datasets/earth_night.py | 102 +++++++++++++++++++++++ pygmt/datasets/load_remote_dataset.py | 21 ++++- pygmt/helpers/caching.py | 1 + pygmt/tests/test_datasets_earth_night.py | 41 +++++++++ 7 files changed, 171 insertions(+), 1 deletion(-) create mode 100644 pygmt/datasets/earth_night.py create mode 100644 pygmt/tests/test_datasets_earth_night.py diff --git a/doc/api/index.rst b/doc/api/index.rst index 0ea40c565e9..01fac7d7a89 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -218,6 +218,7 @@ and store them in GMT's user data directory. :toctree: generated datasets.list_sample_data + datasets.load_black_marble datasets.load_blue_marble datasets.load_earth_age datasets.load_earth_free_air_anomaly diff --git a/pygmt/datasets/__init__.py b/pygmt/datasets/__init__.py index 3aeb56d44fd..d70eec5a1de 100644 --- a/pygmt/datasets/__init__.py +++ b/pygmt/datasets/__init__.py @@ -10,6 +10,7 @@ from pygmt.datasets.earth_geoid import load_earth_geoid from pygmt.datasets.earth_magnetic_anomaly import load_earth_magnetic_anomaly from pygmt.datasets.earth_mask import load_earth_mask +from pygmt.datasets.earth_night import load_black_marble from pygmt.datasets.earth_relief import load_earth_relief from pygmt.datasets.earth_vertical_gravity_gradient import ( load_earth_vertical_gravity_gradient, diff --git a/pygmt/datasets/earth_day.py b/pygmt/datasets/earth_day.py index 2e62d2a1b25..f95a26f7071 100644 --- a/pygmt/datasets/earth_day.py +++ b/pygmt/datasets/earth_day.py @@ -5,12 +5,17 @@ The images are available in various resolutions. """ +import contextlib from collections.abc import Sequence from typing import Literal import xarray as xr from pygmt.datasets.load_remote_dataset import _load_remote_dataset +with contextlib.suppress(ImportError): + # rioxarray is needed to register the rio accessor + import rioxarray # noqa: F401 + __doctest_skip__ = ["load_blue_marble"] diff --git a/pygmt/datasets/earth_night.py b/pygmt/datasets/earth_night.py new file mode 100644 index 00000000000..39742c8124e --- /dev/null +++ b/pygmt/datasets/earth_night.py @@ -0,0 +1,102 @@ +""" +Function to download the NASA Black Marble image datasets from the GMT data server, and +load as :class:`xarray.DataArray`. + +The images are available in various resolutions. +""" + +import contextlib +from collections.abc import Sequence +from typing import Literal + +import xarray as xr +from pygmt.datasets.load_remote_dataset import _load_remote_dataset + +with contextlib.suppress(ImportError): + # rioxarray is needed to register the rio accessor + import rioxarray # noqa: F401 + +__doctest_skip__ = ["load_black_marble"] + + +def load_black_marble( + resolution: Literal[ + "01d", + "30m", + "20m", + "15m", + "10m", + "06m", + "05m", + "04m", + "03m", + "02m", + "01m", + "30s", + ] = "01d", + region: Sequence[float] | str | None = None, +) -> xr.DataArray: + r""" + Load NASA Black Marble images in various resolutions. + + .. figure:: https://www.generic-mapping-tools.org/remote-datasets/_images/GMT_earth_daynight.jpg + :width: 80% + :align: center + + Earth day/night dataset. + + The images are downloaded to a user data directory (usually + ``~/.gmt/server/earth/earth_night/``) the first time you invoke this function. + Afterwards, it will load the image from the data directory. So you'll need an + internet connection the first time around. + + These images can also be accessed by passing in the file name + **@earth_night**\_\ *res* to any image processing function or plotting method. *res* + is the image resolution (see below). + + Refer to :gmt-datasets:`earth-daynight.html` for more details about available + datasets, including version information and references. + + Parameters + ---------- + resolution + The image resolution. The suffix ``d``, ``m``, and ``s`` stand for arc-degree, + arc-minute, and arc-second. + + region + The subregion of the image to load, in the form of a sequence [*xmin*, *xmax*, + *ymin*, *ymax*]. + + Returns + ------- + image + The NASA Black Marble image. Coordinates are latitude and longitude in degrees. + + Note + ---- + The registration and coordinate system type of the returned + :class:`xarray.DataArray` image can be accessed via the GMT accessors (i.e., + ``image.gmt.registration`` and ``image.gmt.gtype`` respectively). However, these + properties may be lost after specific image operations (such as slicing) and will + need to be manually set before passing the image to any PyGMT data processing or + plotting functions. Refer to :class:`pygmt.GMTDataArrayAccessor` for detailed + explanations and workarounds. + + Examples + -------- + + >>> from pygmt.datasets import load_black_marble + >>> # load the default image (pixel-registered 1 arc-degree image) + >>> image = load_black_marble() + """ + image = _load_remote_dataset( + name="earth_night", + prefix="earth_night", + resolution=resolution, + region=region, + registration="pixel", + ) + # If rioxarray is installed, set the coordinate reference system + if hasattr(image, "rio"): + image = image.rio.write_crs(input_crs="OGC:CRS84") + return image diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py index e1d12596f86..d3776f1b3fa 100644 --- a/pygmt/datasets/load_remote_dataset.py +++ b/pygmt/datasets/load_remote_dataset.py @@ -214,6 +214,25 @@ class GMTRemoteDataset(NamedTuple): "15s": Resolution("15s"), }, ), + "earth_night": GMTRemoteDataset( + description="NASA Night Images", + units=None, + extra_attributes={"long_name": "black_marble", "horizontal_datum": "WGS84"}, + resolutions={ + "01d": Resolution("01d", registrations=["pixel"]), + "30m": Resolution("30m", registrations=["pixel"]), + "20m": Resolution("20m", registrations=["pixel"]), + "15m": Resolution("15m", registrations=["pixel"]), + "10m": Resolution("10m", registrations=["pixel"]), + "06m": Resolution("06m", registrations=["pixel"]), + "05m": Resolution("05m", registrations=["pixel"]), + "04m": Resolution("04m", registrations=["pixel"]), + "03m": Resolution("03m", registrations=["pixel"]), + "02m": Resolution("02m", registrations=["pixel"]), + "01m": Resolution("01m", registrations=["pixel"]), + "30s": Resolution("30s", registrations=["pixel"]), + }, + ), "earth_vgg": GMTRemoteDataset( description="IGPP Earth vertical gravity gradient", units="Eotvos", @@ -428,7 +447,7 @@ def _load_remote_dataset( f"'region' is required for {dataset.description} resolution '{resolution}'." ) - kind = "image" if name in {"earth_day"} else "grid" + kind = "image" if name in {"earth_day", "earth_night"} else "grid" kwdict = { "R": region, # region can be None "T": "i" if kind == "image" else "g", diff --git a/pygmt/helpers/caching.py b/pygmt/helpers/caching.py index 714f12d890e..cfff7a28393 100644 --- a/pygmt/helpers/caching.py +++ b/pygmt/helpers/caching.py @@ -22,6 +22,7 @@ def cache_data(): "@earth_mag_01d_g", "@earth_mag4km_01d_g", "@earth_mask_01d_g", + "@earth_night_01d", "@earth_relief_01d_g", "@earth_relief_01d_p", "@earth_relief_10m_g", diff --git a/pygmt/tests/test_datasets_earth_night.py b/pygmt/tests/test_datasets_earth_night.py new file mode 100644 index 00000000000..dd6b3374033 --- /dev/null +++ b/pygmt/tests/test_datasets_earth_night.py @@ -0,0 +1,41 @@ +""" +Test basic functionality for loading Black Marble datasets. +""" + +import numpy as np +import numpy.testing as npt +from pygmt.datasets import load_black_marble + + +def test_black_marble_01d(): + """ + Test some properties of the Black Marble 01d data. + """ + data = load_black_marble(resolution="01d") + assert data.name == "z" + assert data.long_name == "black_marble" + assert data.attrs["horizontal_datum"] == "WGS84" + assert data.attrs["description"] == "NASA Night Images" + assert data.shape == (3, 180, 360) + assert data.dtype == "uint8" + assert data.gmt.registration == 1 + assert data.gmt.gtype == 1 + npt.assert_allclose(data.y, np.arange(89.5, -90.5, -1)) + npt.assert_allclose(data.x, np.arange(-179.5, 180.5, 1)) + npt.assert_allclose(data.min(), 3, atol=1) + npt.assert_allclose(data.max(), 174, atol=1) + + +def test_black_marble_01d_with_region(): + """ + Test loading low-resolution Black Marble with 'region'. + """ + data = load_black_marble(resolution="01d", region=[-10, 10, -5, 5]) + assert data.shape == (3, 10, 20) + assert data.dtype == "uint8" + assert data.gmt.registration == 1 + assert data.gmt.gtype == 1 + npt.assert_allclose(data.y, np.arange(4.5, -5.5, -1)) + npt.assert_allclose(data.x, np.arange(-9.5, 10.5, 1)) + npt.assert_allclose(data.min(), 3, atol=1) + npt.assert_allclose(data.max(), 40, atol=1)