From 66f258aa8bb98d55b6768067ad00b60cd06c89b4 Mon Sep 17 00:00:00 2001 From: Will Schlitzer Date: Sun, 29 Sep 2024 22:25:40 -0400 Subject: [PATCH] Add function to load Blue Marble dataset (#2235) Function to load the `@earth_day_` NASA Blue Marble mosaic RGB images stored as GeoTIFF files. --------- Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com> Co-authored-by: Dongdong Tian --- doc/api/index.rst | 1 + pygmt/datasets/__init__.py | 1 + pygmt/datasets/earth_day.py | 97 ++++++++++++++++++++++++++ pygmt/datasets/load_remote_dataset.py | 30 ++++++-- pygmt/tests/test_datasets_earth_day.py | 41 +++++++++++ pygmt/tests/test_grdimage_image.py | 12 ++-- 6 files changed, 171 insertions(+), 11 deletions(-) create mode 100644 pygmt/datasets/earth_day.py create mode 100644 pygmt/tests/test_datasets_earth_day.py diff --git a/doc/api/index.rst b/doc/api/index.rst index 1c80d8d163f..0ea40c565e9 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_blue_marble datasets.load_earth_age datasets.load_earth_free_air_anomaly datasets.load_earth_geoid diff --git a/pygmt/datasets/__init__.py b/pygmt/datasets/__init__.py index da8fe741c12..3aeb56d44fd 100644 --- a/pygmt/datasets/__init__.py +++ b/pygmt/datasets/__init__.py @@ -5,6 +5,7 @@ """ from pygmt.datasets.earth_age import load_earth_age +from pygmt.datasets.earth_day import load_blue_marble from pygmt.datasets.earth_free_air_anomaly import load_earth_free_air_anomaly from pygmt.datasets.earth_geoid import load_earth_geoid from pygmt.datasets.earth_magnetic_anomaly import load_earth_magnetic_anomaly diff --git a/pygmt/datasets/earth_day.py b/pygmt/datasets/earth_day.py new file mode 100644 index 00000000000..2e62d2a1b25 --- /dev/null +++ b/pygmt/datasets/earth_day.py @@ -0,0 +1,97 @@ +""" +Function to download the NASA Blue Marble image datasets from the GMT data server, and +load as :class:`xarray.DataArray`. + +The images are available in various resolutions. +""" + +from collections.abc import Sequence +from typing import Literal + +import xarray as xr +from pygmt.datasets.load_remote_dataset import _load_remote_dataset + +__doctest_skip__ = ["load_blue_marble"] + + +def load_blue_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 Blue 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_day/``) 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_day**\_\ *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 Blue 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_blue_marble + >>> # load the default image (pixel-registered 1 arc-degree image) + >>> image = load_blue_marble() + """ + image = _load_remote_dataset( + name="earth_day", + prefix="earth_day", + 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 b03908414a5..e1d12596f86 100644 --- a/pygmt/datasets/load_remote_dataset.py +++ b/pygmt/datasets/load_remote_dataset.py @@ -78,6 +78,25 @@ class GMTRemoteDataset(NamedTuple): "01m": Resolution("01m", registrations=["gridline"], tiled=True), }, ), + "earth_day": GMTRemoteDataset( + description="NASA Day Images", + units=None, + extra_attributes={"long_name": "blue_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_faa": GMTRemoteDataset( description="IGPP Earth free-air anomaly", units="mGal", @@ -409,15 +428,18 @@ def _load_remote_dataset( f"'region' is required for {dataset.description} resolution '{resolution}'." ) - # Currently, only grids are supported. Will support images in the future. - kwdict = {"T": "g", "R": region} # region can be None + kind = "image" if name in {"earth_day"} else "grid" + kwdict = { + "R": region, # region can be None + "T": "i" if kind == "image" else "g", + } with Session() as lib: - with lib.virtualfile_out(kind="grid") as voutgrd: + with lib.virtualfile_out(kind=kind) as voutgrd: lib.call_module( module="read", args=[fname, voutgrd, *build_arg_list(kwdict)], ) - grid = lib.virtualfile_to_raster(outgrid=None, vfname=voutgrd) + grid = lib.virtualfile_to_raster(kind=kind, outgrid=None, vfname=voutgrd) # Full path to the grid if not tiled grids. source = which(fname, download="a") if not resinfo.tiled else None diff --git a/pygmt/tests/test_datasets_earth_day.py b/pygmt/tests/test_datasets_earth_day.py new file mode 100644 index 00000000000..f04d7f7f970 --- /dev/null +++ b/pygmt/tests/test_datasets_earth_day.py @@ -0,0 +1,41 @@ +""" +Test basic functionality for loading Blue Marble datasets. +""" + +import numpy as np +import numpy.testing as npt +from pygmt.datasets import load_blue_marble + + +def test_blue_marble_01d(): + """ + Test some properties of the Blue Marble 01d data. + """ + data = load_blue_marble(resolution="01d") + assert data.name == "z" + assert data.long_name == "blue_marble" + assert data.attrs["horizontal_datum"] == "WGS84" + assert data.attrs["description"] == "NASA Day 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(), 10, atol=1) + npt.assert_allclose(data.max(), 255, atol=1) + + +def test_blue_marble_01d_with_region(): + """ + Test loading low-resolution Blue Marble with 'region'. + """ + data = load_blue_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(), 10, atol=1) + npt.assert_allclose(data.max(), 77, atol=1) diff --git a/pygmt/tests/test_grdimage_image.py b/pygmt/tests/test_grdimage_image.py index c7d74b206ed..e2ec8980d18 100644 --- a/pygmt/tests/test_grdimage_image.py +++ b/pygmt/tests/test_grdimage_image.py @@ -3,7 +3,8 @@ """ import pytest -from pygmt import Figure, which +from pygmt import Figure +from pygmt.datasets import load_blue_marble rioxarray = pytest.importorskip("rioxarray") @@ -14,12 +15,9 @@ def fixture_xr_image(): Load the image data from Blue Marble as an xarray.DataArray with shape {"band": 3, "y": 180, "x": 360}. """ - geotiff = which(fname="@earth_day_01d_p", download="c") - with rioxarray.open_rasterio(filename=geotiff) as rda: - if len(rda.band) == 3: - xr_image = rda.load() - assert xr_image.sizes == {"band": 3, "y": 180, "x": 360} - return xr_image + xr_image = load_blue_marble(resolution="01d") + assert xr_image.sizes == {"band": 3, "y": 180, "x": 360} + return xr_image @pytest.mark.mpl_image_compare