From 9bde7f7597fcad257a5e001c93d63de431d2e27c Mon Sep 17 00:00:00 2001
From: Will Schlitzer <schlitzer90@gmail.com>
Date: Wed, 7 Dec 2022 07:59:03 -0500
Subject: [PATCH 01/31] initial commit for blue_marble dataset

---
 pygmt/datasets/__init__.py            |  1 +
 pygmt/datasets/blue_marble.py         | 55 +++++++++++++++++++++++++++
 pygmt/datasets/load_remote_dataset.py | 21 ++++++++++
 3 files changed, 77 insertions(+)
 create mode 100644 pygmt/datasets/blue_marble.py

diff --git a/pygmt/datasets/__init__.py b/pygmt/datasets/__init__.py
index 93b7210a73e..e5f51998f1c 100644
--- a/pygmt/datasets/__init__.py
+++ b/pygmt/datasets/__init__.py
@@ -2,6 +2,7 @@
 #
 # Load sample data included with GMT (downloaded from the GMT cache server).
 
+from pygmt.datasets.blue_marble import load_blue_marble
 from pygmt.datasets.earth_age import load_earth_age
 from pygmt.datasets.earth_relief import load_earth_relief
 from pygmt.datasets.samples import (
diff --git a/pygmt/datasets/blue_marble.py b/pygmt/datasets/blue_marble.py
new file mode 100644
index 00000000000..590280a75b1
--- /dev/null
+++ b/pygmt/datasets/blue_marble.py
@@ -0,0 +1,55 @@
+"""
+Function to download the NASA Blue Marble image datasets 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
+
+__doctest_skip__ = ["load_earth_age"]
+
+
+@kwargs_to_strings(region="sequence")
+def load_blue_marble(resolution="01d"):
+    r"""
+    Load NASA Blue Marble images 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-datasets:`earth-daynight.html` for more details.
+
+    Parameters
+    ----------
+    resolution : str
+        The grid resolution. The suffix ``d``, ``m``, and ``s`` stand for
+        arc-degree, arc-minute, and arc-second. It can be ``"01d"``, ``"30m"``,
+        ``"20m"``, ``"15m"``, ``"10m"``, ``"06m"``, ``"05m"``, ``"04m"``,
+        ``"03m"``, ``"02m"``, ``"01m"``, or ``"30s"``..
+
+
+    Returns
+    -------
+    grid : :class:`xarray.DataArray`
+        The NASA Blue Marble grid. Coordinates are latitude and
+        longitude in degrees.
+    """
+    dataset_prefix = "earth_day_"
+    dataset_name = "blue_marble"
+    grid = _load_remote_dataset(
+        dataset_name=dataset_name,
+        dataset_prefix=dataset_prefix,
+        resolution=resolution,
+        region=None,
+        registration=None,
+    )
+    return grid
diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py
index a60ccb85423..2abefb8f786 100644
--- a/pygmt/datasets/load_remote_dataset.py
+++ b/pygmt/datasets/load_remote_dataset.py
@@ -108,6 +108,27 @@ class GMTRemoteDataset(NamedTuple):
             "01m": Resolution(["gridline"], True),
         },
     ),
+    "blue_marble": GMTRemoteDataset(
+        title="NASA Day/Night Images",
+        name="blue_marble",
+        long_name="NASA Day/Night Images",
+        units="N/A",
+        extra_attributes={"horizontal_datum": "WGS84"},
+        resolutions={
+            "01d": Resolution([None], False),
+            "30m": Resolution([None], False),
+            "20m": Resolution([None], False),
+            "15m": Resolution([None], False),
+            "10m": Resolution([None], False),
+            "06m": Resolution([None], False),
+            "05m": Resolution([None], False),
+            "04m": Resolution([None], False),
+            "03m": Resolution([None], False),
+            "02m": Resolution([None], False),
+            "01m": Resolution([None], False),
+            "30s": Resolution([None], False),
+        },
+    ),
 }
 
 

From d8c8759efc37a879737030e52b734f2bf698a70c Mon Sep 17 00:00:00 2001
From: Will Schlitzer <schlitzer90@gmail.com>
Date: Wed, 7 Dec 2022 08:02:24 -0500
Subject: [PATCH 02/31] update docstring

---
 pygmt/datasets/blue_marble.py | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/pygmt/datasets/blue_marble.py b/pygmt/datasets/blue_marble.py
index 590280a75b1..95302c2dbf3 100644
--- a/pygmt/datasets/blue_marble.py
+++ b/pygmt/datasets/blue_marble.py
@@ -16,15 +16,13 @@ def load_blue_marble(resolution="01d"):
     Load NASA Blue Marble images in various resolutions.
 
     The grids are downloaded to a user data directory
-    (usually ``~/.gmt/server/earth/earth_age/``) the first time you invoke
+    (usually ``~/.gmt/server/earth/earth_day/``) 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).
+    **@earth_day**\_\ *res*\ to any grid plotting/processing
+    function. *res* is the grid resolution (see below).
 
     Refer to :gmt-datasets:`earth-daynight.html` for more details.
 

From 093fc26b34e8cccbb481edf2df7748bcfc07fe5b Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Sat, 17 Feb 2024 15:46:20 +1300
Subject: [PATCH 03/31] Rename load_blue_marble.py to earth_daynight.py and
 update docs

---
 pygmt/datasets/__init__.py                    |  2 +-
 .../{blue_marble.py => earth_daynight.py}     | 33 ++++++++++---------
 2 files changed, 19 insertions(+), 16 deletions(-)
 rename pygmt/datasets/{blue_marble.py => earth_daynight.py} (56%)

diff --git a/pygmt/datasets/__init__.py b/pygmt/datasets/__init__.py
index 37333578bed..fe0636245fd 100644
--- a/pygmt/datasets/__init__.py
+++ b/pygmt/datasets/__init__.py
@@ -4,8 +4,8 @@
 Data are downloaded from the GMT data server.
 """
 
-from pygmt.datasets.blue_marble import load_blue_marble
 from pygmt.datasets.earth_age import load_earth_age
+from pygmt.datasets.earth_daynight 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/blue_marble.py b/pygmt/datasets/earth_daynight.py
similarity index 56%
rename from pygmt/datasets/blue_marble.py
rename to pygmt/datasets/earth_daynight.py
index 95302c2dbf3..e0ad9a39b69 100644
--- a/pygmt/datasets/blue_marble.py
+++ b/pygmt/datasets/earth_daynight.py
@@ -7,7 +7,7 @@
 from pygmt.datasets.load_remote_dataset import _load_remote_dataset
 from pygmt.helpers import kwargs_to_strings
 
-__doctest_skip__ = ["load_earth_age"]
+__doctest_skip__ = ["load_blue_marble"]
 
 
 @kwargs_to_strings(region="sequence")
@@ -15,16 +15,23 @@ def load_blue_marble(resolution="01d"):
     r"""
     Load NASA Blue Marble images in various resolutions.
 
-    The grids 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 grid from the data directory.
-    So you'll need an internet connection the first time around.
+    .. figure:: https://www.generic-mapping-tools.org/remote-datasets/_images/GMT_earth_daynight.jpg
+       :width: 80%
+       :align: center
+
+       Earth day/night dataset.
+
+    The grids 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 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_day**\_\ *res*\ to any grid plotting/processing
-    function. *res* is the grid resolution (see below).
+    **@earth_day**\_\ *res* to any grid processing function or plotting method. *res*
+    is the grid resolution (see below).
 
-    Refer to :gmt-datasets:`earth-daynight.html` for more details.
+    Refer to :gmt-datasets:`earth-daynight.html` for more details about available
+    datasets, including version information and references.
 
     Parameters
     ----------
@@ -34,18 +41,14 @@ def load_blue_marble(resolution="01d"):
         ``"20m"``, ``"15m"``, ``"10m"``, ``"06m"``, ``"05m"``, ``"04m"``,
         ``"03m"``, ``"02m"``, ``"01m"``, or ``"30s"``..
 
-
     Returns
     -------
     grid : :class:`xarray.DataArray`
-        The NASA Blue Marble grid. Coordinates are latitude and
-        longitude in degrees.
+        The NASA Blue Marble grid. Coordinates are latitude and longitude in degrees.
     """
-    dataset_prefix = "earth_day_"
-    dataset_name = "blue_marble"
     grid = _load_remote_dataset(
-        dataset_name=dataset_name,
-        dataset_prefix=dataset_prefix,
+        dataset_name="blue_marble",
+        dataset_prefix="earth_day_",
         resolution=resolution,
         region=None,
         registration=None,

From 2f294581e2aee18e63421534806ff415d14bcae6 Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Sat, 17 Feb 2024 15:48:12 +1300
Subject: [PATCH 04/31] Add load_blue_marble funtion to doc/api/index.rst

---
 doc/api/index.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/doc/api/index.rst b/doc/api/index.rst
index 76b1a19fbe0..dffe3f661cd 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -219,6 +219,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

From 1c99f21cf265ef36380e0dc70d7183876f0c535c Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Sat, 17 Feb 2024 15:58:07 +1300
Subject: [PATCH 05/31] Fix mypy errors in load_remote_dataset.py

Need modify arguments into the Resolution class.
---
 pygmt/datasets/load_remote_dataset.py | 42 +++++++++++++--------------
 1 file changed, 21 insertions(+), 21 deletions(-)

diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py
index 1c96a91d843..ef0737e7f98 100644
--- a/pygmt/datasets/load_remote_dataset.py
+++ b/pygmt/datasets/load_remote_dataset.py
@@ -83,6 +83,27 @@ class GMTRemoteDataset(NamedTuple):
             "01m": Resolution("01m", registrations=["gridline"], tiled=True),
         },
     ),
+    "earth_day": GMTRemoteDataset(
+        title="NASA Day Images",
+        name="blue_marble",
+        long_name="NASA Day Images",
+        units="N/A",
+        extra_attributes={"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_free_air_anomaly": GMTRemoteDataset(
         title="free air anomaly",
         name="free_air_anomaly",
@@ -334,27 +355,6 @@ class GMTRemoteDataset(NamedTuple):
             "01m": Resolution("01m", registrations=["gridline"], tiled=True),
         },
     ),
-    "blue_marble": GMTRemoteDataset(
-        title="NASA Day/Night Images",
-        name="blue_marble",
-        long_name="NASA Day/Night Images",
-        units="N/A",
-        extra_attributes={"horizontal_datum": "WGS84"},
-        resolutions={
-            "01d": Resolution([None], False),
-            "30m": Resolution([None], False),
-            "20m": Resolution([None], False),
-            "15m": Resolution([None], False),
-            "10m": Resolution([None], False),
-            "06m": Resolution([None], False),
-            "05m": Resolution([None], False),
-            "04m": Resolution([None], False),
-            "03m": Resolution([None], False),
-            "02m": Resolution([None], False),
-            "01m": Resolution([None], False),
-            "30s": Resolution([None], False),
-        },
-    ),
 }
 
 

From e839ef79fb4d0176729e09a338fc69a3d1ed5b8b Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Sat, 17 Feb 2024 16:08:14 +1300
Subject: [PATCH 06/31] Add doctest and rename returned output from grid to
 image

A short doctest to test loading the Blue Marble dataset, and rename grid to image since a 3-band RGB image is returned.
---
 pygmt/datasets/earth_daynight.py | 18 +++++++++++++-----
 1 file changed, 13 insertions(+), 5 deletions(-)

diff --git a/pygmt/datasets/earth_daynight.py b/pygmt/datasets/earth_daynight.py
index e0ad9a39b69..5ead9d4932d 100644
--- a/pygmt/datasets/earth_daynight.py
+++ b/pygmt/datasets/earth_daynight.py
@@ -5,12 +5,13 @@
 The grids are available in various resolutions.
 """
 from pygmt.datasets.load_remote_dataset import _load_remote_dataset
-from pygmt.helpers import kwargs_to_strings
+
+# from pygmt.helpers import kwargs_to_strings
 
 __doctest_skip__ = ["load_blue_marble"]
 
 
-@kwargs_to_strings(region="sequence")
+# @kwargs_to_strings(region="sequence")
 def load_blue_marble(resolution="01d"):
     r"""
     Load NASA Blue Marble images in various resolutions.
@@ -43,11 +44,18 @@ def load_blue_marble(resolution="01d"):
 
     Returns
     -------
-    grid : :class:`xarray.DataArray`
-        The NASA Blue Marble grid. Coordinates are latitude and longitude in degrees.
+    image : :class:`xarray.DataArray`
+        The NASA Blue Marble image. Coordinates are latitude and longitude in degrees.
+
+    Examples
+    --------
+
+    >>> from pygmt.datasets import load_blue_marble
+    >>> # load the default image (pixel-registered 1 arc-degree image)
+    >>> image = load_blue_marble()
     """
     grid = _load_remote_dataset(
-        dataset_name="blue_marble",
+        dataset_name="earth_day",
         dataset_prefix="earth_day_",
         resolution=resolution,
         region=None,

From aa8a47477a13f22889491e00093da9488f091557 Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Sat, 17 Feb 2024 16:34:44 +1300
Subject: [PATCH 07/31] Use rioxarray.open_rasterio to open earth_day images

---
 pygmt/datasets/load_remote_dataset.py |  3 ++-
 pygmt/io.py                           | 10 +++++++++-
 2 files changed, 11 insertions(+), 2 deletions(-)

diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py
index ef0737e7f98..190e8f38c6d 100644
--- a/pygmt/datasets/load_remote_dataset.py
+++ b/pygmt/datasets/load_remote_dataset.py
@@ -435,7 +435,8 @@ def _load_remote_dataset(
                 f"'region' is required for {dataset.title} resolution '{resolution}'."
             )
         fname = which(f"@{dataset_prefix}{resolution}{reg}", download="a")
-        grid = load_dataarray(fname, engine="netcdf4")
+        engine = "rasterio" if dataset_prefix == "earth_day_" else "netcdf4"
+        grid = load_dataarray(fname, engine=engine)
     else:
         grid = grdcut(f"@{dataset_prefix}{resolution}{reg}", region=region)
 
diff --git a/pygmt/io.py b/pygmt/io.py
index 9c27e7b08b5..e612da5a692 100644
--- a/pygmt/io.py
+++ b/pygmt/io.py
@@ -39,7 +39,15 @@ def load_dataarray(filename_or_obj, **kwargs):
     if "cache" in kwargs:
         raise TypeError("cache has no effect in this context")
 
-    with xr.open_dataarray(filename_or_obj, **kwargs) as dataarray:
+    if kwargs.get("engine") == "rasterio":
+        import rioxarray
+
+        open_dataarray_func = rioxarray.open_rasterio
+        kwargs.pop("engine")
+    else:
+        open_dataarray_func = xr.open_dataarray
+
+    with open_dataarray_func(filename_or_obj, **kwargs) as dataarray:
         result = dataarray.load()
         _ = result.gmt  # load GMTDataArray accessor information
 

From cd4cab6fe6a07bc0b3fd9b2dfb7711d19dea7b2a Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Sat, 17 Feb 2024 16:35:14 +1300
Subject: [PATCH 08/31] Refactor fixture_xr_image to use load_blue_marble
 function

---
 pygmt/tests/test_grdimage_image.py | 12 +++++-------
 1 file changed, 5 insertions(+), 7 deletions(-)

diff --git a/pygmt/tests/test_grdimage_image.py b/pygmt/tests/test_grdimage_image.py
index c713957dc91..0ea47cd238b 100644
--- a/pygmt/tests/test_grdimage_image.py
+++ b/pygmt/tests/test_grdimage_image.py
@@ -2,7 +2,8 @@
 Test Figure.grdimage on 3-band RGB images.
 """
 import pytest
-from pygmt import Figure, which
+from pygmt import Figure
+from pygmt.datasets import load_blue_marble
 
 rioxarray = pytest.importorskip("rioxarray")
 
@@ -13,12 +14,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

From fc060eaca825d2dba5618b2cf29bfe63e6e4d529 Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Sun, 18 Feb 2024 00:49:24 -0500
Subject: [PATCH 09/31] Only load GMTDataArray accessor info when
 engine!=rasterio

Prevent segfault when loading accessor info from an xarray.Datarray loaded using rioxarray.open_rasterio.
---
 pygmt/io.py | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/pygmt/io.py b/pygmt/io.py
index e612da5a692..2efbed24d30 100644
--- a/pygmt/io.py
+++ b/pygmt/io.py
@@ -42,13 +42,12 @@ def load_dataarray(filename_or_obj, **kwargs):
     if kwargs.get("engine") == "rasterio":
         import rioxarray
 
-        open_dataarray_func = rioxarray.open_rasterio
         kwargs.pop("engine")
+        with rioxarray.open_rasterio(filename_or_obj, **kwargs) as dataarray:
+            result = dataarray.load()
     else:
-        open_dataarray_func = xr.open_dataarray
-
-    with open_dataarray_func(filename_or_obj, **kwargs) as dataarray:
-        result = dataarray.load()
-        _ = result.gmt  # load GMTDataArray accessor information
+        with xr.open_dataarray(filename_or_obj, **kwargs) as dataarray:
+            result = dataarray.load()
+            _ = result.gmt  # load GMTDataArray accessor information
 
     return result

From c3017d54bb0129f3478c174f3795b09ea45b6bb7 Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Sun, 18 Feb 2024 10:40:46 -0500
Subject: [PATCH 10/31] Set tiled=True for spatial resolutions 05m and higher

Also set units=None, and renamed 'grid' to 'image' in the docstrings.
---
 pygmt/datasets/earth_daynight.py      | 24 ++++++++++++------------
 pygmt/datasets/load_remote_dataset.py | 14 +++++++-------
 2 files changed, 19 insertions(+), 19 deletions(-)

diff --git a/pygmt/datasets/earth_daynight.py b/pygmt/datasets/earth_daynight.py
index 5ead9d4932d..dd940b9b7f5 100644
--- a/pygmt/datasets/earth_daynight.py
+++ b/pygmt/datasets/earth_daynight.py
@@ -22,14 +22,14 @@ def load_blue_marble(resolution="01d"):
 
        Earth day/night dataset.
 
-    The grids are downloaded to a user data directory (usually
+    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 grid from the data directory. So you'll need an
+    Afterwards, it will load the image 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_day**\_\ *res* to any grid processing function or plotting method. *res*
-    is the grid resolution (see below).
+    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.
@@ -37,10 +37,10 @@ def load_blue_marble(resolution="01d"):
     Parameters
     ----------
     resolution : str
-        The grid resolution. The suffix ``d``, ``m``, and ``s`` stand for
-        arc-degree, arc-minute, and arc-second. It can be ``"01d"``, ``"30m"``,
-        ``"20m"``, ``"15m"``, ``"10m"``, ``"06m"``, ``"05m"``, ``"04m"``,
-        ``"03m"``, ``"02m"``, ``"01m"``, or ``"30s"``..
+        The image resolution. The suffix ``d``, ``m``, and ``s`` stand for arc-degree,
+        arc-minute, and arc-second. It can be ``"01d"``, ``"30m"``, ``"20m"``,
+        ``"15m"``, ``"10m"``, ``"06m"``, ``"05m"``, ``"04m"``, ``"03m"``, ``"02m"``,
+        ``"01m"``, or ``"30s"``.
 
     Returns
     -------
@@ -54,11 +54,11 @@ def load_blue_marble(resolution="01d"):
     >>> # load the default image (pixel-registered 1 arc-degree image)
     >>> image = load_blue_marble()
     """
-    grid = _load_remote_dataset(
+    image = _load_remote_dataset(
         dataset_name="earth_day",
         dataset_prefix="earth_day_",
         resolution=resolution,
         region=None,
-        registration=None,
+        registration="pixel",
     )
-    return grid
+    return image
diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py
index 190e8f38c6d..40668ed8c60 100644
--- a/pygmt/datasets/load_remote_dataset.py
+++ b/pygmt/datasets/load_remote_dataset.py
@@ -87,7 +87,7 @@ class GMTRemoteDataset(NamedTuple):
         title="NASA Day Images",
         name="blue_marble",
         long_name="NASA Day Images",
-        units="N/A",
+        units=None,
         extra_attributes={"horizontal_datum": "WGS84"},
         resolutions={
             "01d": Resolution("01d", registrations=["pixel"]),
@@ -96,12 +96,12 @@ class GMTRemoteDataset(NamedTuple):
             "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"]),
+            "05m": Resolution("05m", registrations=["pixel"], tiled=True),
+            "04m": Resolution("04m", registrations=["pixel"], tiled=True),
+            "03m": Resolution("03m", registrations=["pixel"], tiled=True),
+            "02m": Resolution("02m", registrations=["pixel"], tiled=True),
+            "01m": Resolution("01m", registrations=["pixel"], tiled=True),
+            "30s": Resolution("30s", registrations=["pixel"], tiled=True),
         },
     ),
     "earth_free_air_anomaly": GMTRemoteDataset(

From c3e2e66584f567a1be246e5c04a2b2081df8dc22 Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Sun, 18 Feb 2024 10:41:48 -0500
Subject: [PATCH 11/31] Add unit test to load_blue_marble 01d

---
 pygmt/tests/test_datasets_earth_daynight.py | 23 +++++++++++++++++++++
 1 file changed, 23 insertions(+)
 create mode 100644 pygmt/tests/test_datasets_earth_daynight.py

diff --git a/pygmt/tests/test_datasets_earth_daynight.py b/pygmt/tests/test_datasets_earth_daynight.py
new file mode 100644
index 00000000000..307663bc43e
--- /dev/null
+++ b/pygmt/tests/test_datasets_earth_daynight.py
@@ -0,0 +1,23 @@
+"""
+Test basic functionality for loading Blue and Black 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 == "blue_marble"
+    assert data.attrs["horizontal_datum"] == "WGS84"
+    assert data.attrs["long_name"] == "NASA Day Images"
+    assert data.shape == (3, 180, 360)
+    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)

From a2a18a906c9104a17420b6fccc72869f683de2d1 Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Sun, 18 Feb 2024 10:54:32 -0500
Subject: [PATCH 12/31] Skip test_datasets_load_earth_daynight if rioxarray not
 installed

---
 pygmt/tests/test_datasets_earth_daynight.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/pygmt/tests/test_datasets_earth_daynight.py b/pygmt/tests/test_datasets_earth_daynight.py
index 307663bc43e..0c871e2bece 100644
--- a/pygmt/tests/test_datasets_earth_daynight.py
+++ b/pygmt/tests/test_datasets_earth_daynight.py
@@ -3,8 +3,11 @@
 """
 import numpy as np
 import numpy.testing as npt
+import pytest
 from pygmt.datasets import load_blue_marble
 
+rioxarray = pytest.importorskip("rioxarray")
+
 
 def test_blue_marble_01d():
     """

From 7c3875482ac53a7d1fe8f7fa6775ab42a5be1722 Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Sun, 18 Feb 2024 11:06:16 -0500
Subject: [PATCH 13/31] Add region parameter for load_blue_marble

Also added some failing unit tests, since grdcut doesn't allow cropping GeoTIFFs yet.
---
 pygmt/datasets/earth_daynight.py            | 14 ++++++----
 pygmt/tests/test_datasets_earth_daynight.py | 31 +++++++++++++++++++++
 2 files changed, 40 insertions(+), 5 deletions(-)

diff --git a/pygmt/datasets/earth_daynight.py b/pygmt/datasets/earth_daynight.py
index dd940b9b7f5..566472cf869 100644
--- a/pygmt/datasets/earth_daynight.py
+++ b/pygmt/datasets/earth_daynight.py
@@ -5,14 +5,13 @@
 The grids are available in various resolutions.
 """
 from pygmt.datasets.load_remote_dataset import _load_remote_dataset
-
-# from pygmt.helpers import kwargs_to_strings
+from pygmt.helpers import kwargs_to_strings
 
 __doctest_skip__ = ["load_blue_marble"]
 
 
-# @kwargs_to_strings(region="sequence")
-def load_blue_marble(resolution="01d"):
+@kwargs_to_strings(region="sequence")
+def load_blue_marble(resolution="01d", region=None):
     r"""
     Load NASA Blue Marble images in various resolutions.
 
@@ -42,6 +41,11 @@ def load_blue_marble(resolution="01d"):
         ``"15m"``, ``"10m"``, ``"06m"``, ``"05m"``, ``"04m"``, ``"03m"``, ``"02m"``,
         ``"01m"``, or ``"30s"``.
 
+    region : str or list
+        The subregion of the image to load, in the form of a list [*xmin*, *xmax*,
+        *ymin*, *ymax*] or a string *xmin/xmax/ymin/ymax*. Required for images with
+        resolutions higher than 5 arc-minutes (i.e., ``"05m"``).
+
     Returns
     -------
     image : :class:`xarray.DataArray`
@@ -58,7 +62,7 @@ def load_blue_marble(resolution="01d"):
         dataset_name="earth_day",
         dataset_prefix="earth_day_",
         resolution=resolution,
-        region=None,
+        region=region,
         registration="pixel",
     )
     return image
diff --git a/pygmt/tests/test_datasets_earth_daynight.py b/pygmt/tests/test_datasets_earth_daynight.py
index 0c871e2bece..c7477fc6cfa 100644
--- a/pygmt/tests/test_datasets_earth_daynight.py
+++ b/pygmt/tests/test_datasets_earth_daynight.py
@@ -24,3 +24,34 @@ def test_blue_marble_01d():
     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, 20, 10)
+    assert data.gmt.registration == 1
+    assert data.gmt.gtype == 1
+    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(), 0, atol=1)
+    npt.assert_allclose(data.max(), 255, atol=1)
+
+
+def test_blue_marble_01m_default_registration():
+    """
+    Test that the grid returned by default for the 1 arc-minute resolution has
+    a "pixel" registration.
+    """
+    data = load_blue_marble(resolution="01m", region=[-10, -9, 3, 5])
+    assert data.shape == (3, 12, 3)
+    assert data.gmt.registration == 1
+    assert data.gmt.gtype == 1
+    assert data.coords["lat"].data.min() == 3.0
+    assert data.coords["lat"].data.max() == 5.0
+    assert data.coords["lon"].data.min() == -10.0
+    assert data.coords["lon"].data.max() == -9.0
+    npt.assert_allclose(data.min(), 0, atol=1)
+    npt.assert_allclose(data.max(), 255, atol=1)

From 767d7c7e4557799df39471d690831eff9d0dc16b Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Sat, 16 Mar 2024 18:28:04 +1300
Subject: [PATCH 14/31] Allow grdcut to slice GeoTIFFs so load_blue_marble's
 region param works

New engine parameter in `grdcut` that allows passing in the `engine='rasterio'` to enable slicing the Blue Marble GeoTIFF files. Updated some numbers in test_datasets_earth_daynight.py to get unit tests to pass.
---
 pygmt/datasets/load_remote_dataset.py       | 10 ++++++--
 pygmt/src/grdcut.py                         | 13 +++++++---
 pygmt/tests/test_datasets_earth_daynight.py | 27 ++++++++++++---------
 3 files changed, 33 insertions(+), 17 deletions(-)

diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py
index 2d5df5b3d0d..7984324d397 100644
--- a/pygmt/datasets/load_remote_dataset.py
+++ b/pygmt/datasets/load_remote_dataset.py
@@ -430,16 +430,22 @@ def _load_remote_dataset(
     # different ways to load tiled and non-tiled grids.
     # Known issue: tiled grids don't support slice operation
     # See https://github.com/GenericMappingTools/pygmt/issues/524
+    engine = "rasterio" if dataset_prefix == "earth_day_" else "netcdf4"
     if region is None:
         if dataset.resolutions[resolution].tiled:
             raise GMTInvalidInput(
                 f"'region' is required for {dataset.title} resolution '{resolution}'."
             )
         fname = which(f"@{dataset_prefix}{resolution}{reg}", download="a")
-        engine = "rasterio" if dataset_prefix == "earth_day_" else "netcdf4"
         grid = load_dataarray(fname, engine=engine)
     else:
-        grid = grdcut(f"@{dataset_prefix}{resolution}{reg}", region=region)
+        grid = grdcut(
+            f"@{dataset_prefix}{resolution}{reg}", region=region, engine=engine
+        )
+
+    if registration == "pixel":
+        grid.gmt.registration = 1
+    grid.gmt.gtype = 1  # GMT remote datasets are always geographic?
 
     # Add some metadata to the grid
     grid.name = dataset.name
diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py
index 9ffcda213d6..26f29782878 100644
--- a/pygmt/src/grdcut.py
+++ b/pygmt/src/grdcut.py
@@ -27,7 +27,7 @@
     f="coltypes",
 )
 @kwargs_to_strings(R="sequence")
-def grdcut(grid, **kwargs):
+def grdcut(grid, engine: str = None, **kwargs):
     r"""
     Extract subregion from a grid.
 
@@ -74,6 +74,10 @@ def grdcut(grid, **kwargs):
         NaNs, append **+N** to strip off such columns before (optionally)
         considering the range of the core subset for further reduction of the
         area.
+    engine : str
+        Optional. Installed backend or subclass of xarray.backends.BackendEntrypoint.
+        Engine to use when reading files. If not provided, the default engine is chosen
+        based on available dependencies, with a preference for "netcdf4".
 
     {verbose}
     {coltypes}
@@ -99,7 +103,8 @@ def grdcut(grid, **kwargs):
     >>> # 12° E to 15° E and a latitude range of 21° N to 24° N
     >>> new_grid = pygmt.grdcut(grid=grid, region=[12, 15, 21, 24])
     """
-    with GMTTempFile(suffix=".nc") as tmpfile:
+    suffix = ".tif" if engine == "rasterio" else ".nc"
+    with GMTTempFile(suffix=suffix) as tmpfile:
         with Session() as lib:
             with lib.virtualfile_in(check_kind="raster", data=grid) as vingrd:
                 if (outgrid := kwargs.get("G")) is None:
@@ -108,4 +113,6 @@ def grdcut(grid, **kwargs):
                     module="grdcut", args=build_arg_string(kwargs, infile=vingrd)
                 )
 
-        return load_dataarray(outgrid) if outgrid == tmpfile.name else None
+        return (
+            load_dataarray(outgrid, engine=engine) if outgrid == tmpfile.name else None
+        )
diff --git a/pygmt/tests/test_datasets_earth_daynight.py b/pygmt/tests/test_datasets_earth_daynight.py
index c7477fc6cfa..ff315837b03 100644
--- a/pygmt/tests/test_datasets_earth_daynight.py
+++ b/pygmt/tests/test_datasets_earth_daynight.py
@@ -18,6 +18,7 @@ def test_blue_marble_01d():
     assert data.attrs["horizontal_datum"] == "WGS84"
     assert data.attrs["long_name"] == "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))
@@ -31,13 +32,14 @@ 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, 20, 10)
+    assert data.shape == (3, 10, 20)
+    assert data.dtype == "uint8"
     assert data.gmt.registration == 1
     assert data.gmt.gtype == 1
-    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(), 0, atol=1)
-    npt.assert_allclose(data.max(), 255, atol=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)
 
 
 def test_blue_marble_01m_default_registration():
@@ -46,12 +48,13 @@ def test_blue_marble_01m_default_registration():
     a "pixel" registration.
     """
     data = load_blue_marble(resolution="01m", region=[-10, -9, 3, 5])
-    assert data.shape == (3, 12, 3)
+    assert data.shape == (3, 120, 60)
+    assert data.dtype == "uint8"
     assert data.gmt.registration == 1
     assert data.gmt.gtype == 1
-    assert data.coords["lat"].data.min() == 3.0
-    assert data.coords["lat"].data.max() == 5.0
-    assert data.coords["lon"].data.min() == -10.0
-    assert data.coords["lon"].data.max() == -9.0
-    npt.assert_allclose(data.min(), 0, atol=1)
-    npt.assert_allclose(data.max(), 255, atol=1)
+    assert data.coords["y"].data.min() == 3.008333333333333
+    assert data.coords["y"].data.max() == 4.991666666666666
+    assert data.coords["x"].data.min() == -9.991666666666667
+    assert data.coords["x"].data.max() == -9.008333333333335
+    npt.assert_allclose(data.min(), 10, atol=1)
+    npt.assert_allclose(data.max(), 79, atol=1)

From a8600d35bc6db1e6f7add6a1bb2dca813421e8e5 Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Sat, 16 Mar 2024 18:31:29 +1300
Subject: [PATCH 15/31] Lint

---
 pygmt/datasets/earth_daynight.py            | 1 +
 pygmt/src/grdcut.py                         | 2 +-
 pygmt/tests/test_datasets_earth_daynight.py | 1 +
 3 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/pygmt/datasets/earth_daynight.py b/pygmt/datasets/earth_daynight.py
index 566472cf869..5ffab2b0ed7 100644
--- a/pygmt/datasets/earth_daynight.py
+++ b/pygmt/datasets/earth_daynight.py
@@ -4,6 +4,7 @@
 
 The grids are available in various resolutions.
 """
+
 from pygmt.datasets.load_remote_dataset import _load_remote_dataset
 from pygmt.helpers import kwargs_to_strings
 
diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py
index 26f29782878..2d2e106103a 100644
--- a/pygmt/src/grdcut.py
+++ b/pygmt/src/grdcut.py
@@ -27,7 +27,7 @@
     f="coltypes",
 )
 @kwargs_to_strings(R="sequence")
-def grdcut(grid, engine: str = None, **kwargs):
+def grdcut(grid, engine: str | None = None, **kwargs):
     r"""
     Extract subregion from a grid.
 
diff --git a/pygmt/tests/test_datasets_earth_daynight.py b/pygmt/tests/test_datasets_earth_daynight.py
index ff315837b03..de24753eff9 100644
--- a/pygmt/tests/test_datasets_earth_daynight.py
+++ b/pygmt/tests/test_datasets_earth_daynight.py
@@ -1,6 +1,7 @@
 """
 Test basic functionality for loading Blue and Black Marble datasets.
 """
+
 import numpy as np
 import numpy.testing as npt
 import pytest

From 6eabbf4fbab240843c0c804a2a6d16d09313c636 Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Tue, 18 Jun 2024 20:54:47 +1200
Subject: [PATCH 16/31] Handle image kind in _load_remote_dataset function

Special case `earth_day` to be loaded as GMT_IMAGE rather than GMT_GRID. Need to use `GMT_OUT|GMT_IS_REFERENCE` in virtualfile_out to avoid segfault, xref #3115.
---
 pygmt/clib/session.py                 |  9 +++++++--
 pygmt/datasets/load_remote_dataset.py | 12 +++++++++---
 2 files changed, 16 insertions(+), 5 deletions(-)

diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py
index b0aaff44ec3..e444abbbe04 100644
--- a/pygmt/clib/session.py
+++ b/pygmt/clib/session.py
@@ -1697,7 +1697,9 @@ def virtualfile_from_data(
 
     @contextlib.contextmanager
     def virtualfile_out(
-        self, kind: Literal["dataset", "grid"] = "dataset", fname: str | None = None
+        self,
+        kind: Literal["dataset", "grid", "image"] = "dataset",
+        fname: str | None = None,
     ):
         r"""
         Create a virtual file or an actual file for storing output data.
@@ -1754,8 +1756,11 @@ def virtualfile_out(
             family, geometry = {
                 "dataset": ("GMT_IS_DATASET", "GMT_IS_PLP"),
                 "grid": ("GMT_IS_GRID", "GMT_IS_SURFACE"),
+                "image": ("GMT_IS_IMAGE", "GMT_IS_SURFACE"),
             }[kind]
-            with self.open_virtualfile(family, geometry, "GMT_OUT", None) as vfile:
+            with self.open_virtualfile(
+                family, geometry, "GMT_OUT|GMT_IS_REFERENCE", None
+            ) as vfile:
                 yield vfile
 
     def inquire_virtualfile(self, vfname: str) -> int:
diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py
index 860d1ee1297..b191c024a6c 100644
--- a/pygmt/datasets/load_remote_dataset.py
+++ b/pygmt/datasets/load_remote_dataset.py
@@ -427,14 +427,20 @@ def _load_remote_dataset(
         )
 
     # Currently, only grids are supported. Will support images in the future.
-    kwdict = {"T": "g", "R": region}  # region can be None
+    kwdict = {"R": region}  # region can be None
+    if name in ("earth_day",):
+        kind = "image"
+        kwdict.update({"T": "i"})
+    else:
+        kind = "grid"
+        kwdict.update({"T": "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

From 54184fbbda204140fda2fe4690aca39188bc993f Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Sat, 28 Sep 2024 15:30:51 +1200
Subject: [PATCH 17/31] Update attrs (name, long_name, description) of
 blue_marble

Set long_name to 'blue_marble', check that name is 'z' and description is 'NASA Day Images'.
---
 pygmt/datasets/load_remote_dataset.py       | 2 +-
 pygmt/tests/test_datasets_earth_daynight.py | 5 +++--
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py
index 896942f7c4e..052b9734f2a 100644
--- a/pygmt/datasets/load_remote_dataset.py
+++ b/pygmt/datasets/load_remote_dataset.py
@@ -81,7 +81,7 @@ class GMTRemoteDataset(NamedTuple):
     "earth_day": GMTRemoteDataset(
         description="NASA Day Images",
         units=None,
-        extra_attributes={"horizontal_datum": "WGS84"},
+        extra_attributes={"long_name": "blue_marble", "horizontal_datum": "WGS84"},
         resolutions={
             "01d": Resolution("01d", registrations=["pixel"]),
             "30m": Resolution("30m", registrations=["pixel"]),
diff --git a/pygmt/tests/test_datasets_earth_daynight.py b/pygmt/tests/test_datasets_earth_daynight.py
index de24753eff9..7536560940f 100644
--- a/pygmt/tests/test_datasets_earth_daynight.py
+++ b/pygmt/tests/test_datasets_earth_daynight.py
@@ -15,9 +15,10 @@ def test_blue_marble_01d():
     Test some properties of the Blue Marble 01d data.
     """
     data = load_blue_marble(resolution="01d")
-    assert data.name == "blue_marble"
+    assert data.name == "z"
+    assert data.long_name == "blue_marble"
     assert data.attrs["horizontal_datum"] == "WGS84"
-    assert data.attrs["long_name"] == "NASA Day Images"
+    assert data.attrs["description"] == "NASA Day Images"
     assert data.shape == (3, 180, 360)
     assert data.dtype == "uint8"
     assert data.gmt.registration == 1

From e6c29645ab932af487f17e79df01fd087456b3ce Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Sat, 28 Sep 2024 15:43:48 +1200
Subject: [PATCH 18/31] Add type hints to resolution and region parameters

Following changes at https://github.com/GenericMappingTools/pygmt/pull/3260 and https://github.com/GenericMappingTools/pygmt/pull/3272
---
 pygmt/datasets/earth_daynight.py | 40 ++++++++++++++++++++++----------
 1 file changed, 28 insertions(+), 12 deletions(-)

diff --git a/pygmt/datasets/earth_daynight.py b/pygmt/datasets/earth_daynight.py
index f231d3f345c..fcb78ad445b 100644
--- a/pygmt/datasets/earth_daynight.py
+++ b/pygmt/datasets/earth_daynight.py
@@ -5,14 +5,32 @@
 The grids 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
-from pygmt.helpers import kwargs_to_strings
 
 __doctest_skip__ = ["load_blue_marble"]
 
 
-@kwargs_to_strings(region="sequence")
-def load_blue_marble(resolution="01d", region=None):
+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.
 
@@ -36,20 +54,18 @@ def load_blue_marble(resolution="01d", region=None):
 
     Parameters
     ----------
-    resolution : str
+    resolution
         The image resolution. The suffix ``d``, ``m``, and ``s`` stand for arc-degree,
-        arc-minute, and arc-second. It can be ``"01d"``, ``"30m"``, ``"20m"``,
-        ``"15m"``, ``"10m"``, ``"06m"``, ``"05m"``, ``"04m"``, ``"03m"``, ``"02m"``,
-        ``"01m"``, or ``"30s"``.
+        arc-minute, and arc-second.
 
-    region : str or list
-        The subregion of the image to load, in the form of a list [*xmin*, *xmax*,
-        *ymin*, *ymax*] or a string *xmin/xmax/ymin/ymax*. Required for images with
-        resolutions higher than 5 arc-minutes (i.e., ``"05m"``).
+    region
+        The subregion of the image to load, in the form of a sequence [*xmin*, *xmax*,
+        *ymin*, *ymax*]. Required for images with resolutions higher than 5 arc-minutes
+        (i.e., ``"05m"``).
 
     Returns
     -------
-    image : :class:`xarray.DataArray`
+    image
         The NASA Blue Marble image. Coordinates are latitude and longitude in degrees.
 
     Examples

From 01d1b8014e3f2abf4e349691b286328df65dd377 Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Mon, 30 Sep 2024 08:57:16 +1300
Subject: [PATCH 19/31] Partially revert "Allow grdcut to slice GeoTIFFs so
 load_blue_marble's region param works"

This partially reverts commit 767d7c7e4557799df39471d690831eff9d0dc16b.
---
 pygmt/src/grdcut.py | 13 +++----------
 1 file changed, 3 insertions(+), 10 deletions(-)

diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py
index be5d03b6851..e2e4bc8c91c 100644
--- a/pygmt/src/grdcut.py
+++ b/pygmt/src/grdcut.py
@@ -28,7 +28,7 @@
     f="coltypes",
 )
 @kwargs_to_strings(R="sequence")
-def grdcut(grid, engine: str | None = None, **kwargs) -> xr.DataArray | None:
+def grdcut(grid, **kwargs) -> xr.DataArray | None:
     r"""
     Extract subregion from a grid.
 
@@ -75,10 +75,6 @@ def grdcut(grid, engine: str | None = None, **kwargs) -> xr.DataArray | None:
         NaNs, append **+N** to strip off such columns before (optionally)
         considering the range of the core subset for further reduction of the
         area.
-    engine : str
-        Optional. Installed backend or subclass of xarray.backends.BackendEntrypoint.
-        Engine to use when reading files. If not provided, the default engine is chosen
-        based on available dependencies, with a preference for "netcdf4".
 
     {verbose}
     {coltypes}
@@ -104,8 +100,7 @@ def grdcut(grid, engine: str | None = None, **kwargs) -> xr.DataArray | None:
     >>> # 12° E to 15° E and a latitude range of 21° N to 24° N
     >>> new_grid = pygmt.grdcut(grid=grid, region=[12, 15, 21, 24])
     """
-    suffix = ".tif" if engine == "rasterio" else ".nc"
-    with GMTTempFile(suffix=suffix) as tmpfile:
+    with GMTTempFile(suffix=".nc") as tmpfile:
         with Session() as lib:
             with lib.virtualfile_in(check_kind="raster", data=grid) as vingrd:
                 if (outgrid := kwargs.get("G")) is None:
@@ -114,6 +109,4 @@ def grdcut(grid, engine: str | None = None, **kwargs) -> xr.DataArray | None:
                     module="grdcut", args=build_arg_list(kwargs, infile=vingrd)
                 )
 
-        return (
-            load_dataarray(outgrid, engine=engine) if outgrid == tmpfile.name else None
-        )
+        return load_dataarray(outgrid) if outgrid == tmpfile.name else None

From 0ca210a52f1951327b167f36c644753c3715a7a6 Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Mon, 30 Sep 2024 08:57:33 +1300
Subject: [PATCH 20/31] Revert "Use rioxarray.open_rasterio to open earth_day
 images"

This reverts commit aa8a47477a13f22889491e00093da9488f091557.
---
 pygmt/io.py | 13 +++----------
 1 file changed, 3 insertions(+), 10 deletions(-)

diff --git a/pygmt/io.py b/pygmt/io.py
index a6f368f5f72..33f3bfdbb12 100644
--- a/pygmt/io.py
+++ b/pygmt/io.py
@@ -40,15 +40,8 @@ def load_dataarray(filename_or_obj, **kwargs):
     if "cache" in kwargs:
         raise TypeError("cache has no effect in this context")
 
-    if kwargs.get("engine") == "rasterio":
-        import rioxarray
-
-        kwargs.pop("engine")
-        with rioxarray.open_rasterio(filename_or_obj, **kwargs) as dataarray:
-            result = dataarray.load()
-    else:
-        with xr.open_dataarray(filename_or_obj, **kwargs) as dataarray:
-            result = dataarray.load()
-            _ = result.gmt  # load GMTDataArray accessor information
+    with xr.open_dataarray(filename_or_obj, **kwargs) as dataarray:
+        result = dataarray.load()
+        _ = result.gmt  # load GMTDataArray accessor information
 
     return result

From 28d668d7ccf34f508a62ef2b2cc7a26a6ea0414b Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Mon, 30 Sep 2024 09:07:36 +1300
Subject: [PATCH 21/31] Use npt.assert_allclose in
 test_blue_marble_01m_default_registration

---
 pygmt/tests/test_datasets_earth_daynight.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/pygmt/tests/test_datasets_earth_daynight.py b/pygmt/tests/test_datasets_earth_daynight.py
index 7536560940f..215ff3390d8 100644
--- a/pygmt/tests/test_datasets_earth_daynight.py
+++ b/pygmt/tests/test_datasets_earth_daynight.py
@@ -54,9 +54,9 @@ def test_blue_marble_01m_default_registration():
     assert data.dtype == "uint8"
     assert data.gmt.registration == 1
     assert data.gmt.gtype == 1
-    assert data.coords["y"].data.min() == 3.008333333333333
-    assert data.coords["y"].data.max() == 4.991666666666666
-    assert data.coords["x"].data.min() == -9.991666666666667
-    assert data.coords["x"].data.max() == -9.008333333333335
+    npt.assert_allclose(data.coords["y"].data.min(), 3.0083333333333333)
+    npt.assert_allclose(data.coords["y"].data.max(), 4.991666666666666)
+    npt.assert_allclose(data.coords["x"].data.min(), -9.991666666666667)
+    npt.assert_allclose(data.coords["x"].data.max(), -9.008333333333335)
     npt.assert_allclose(data.min(), 10, atol=1)
     npt.assert_allclose(data.max(), 79, atol=1)

From a6a08e30cc99e33bba21f6854817cdf8bde4e3fb Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Mon, 30 Sep 2024 10:11:21 +1300
Subject: [PATCH 22/31] Plot xr_image in same fashion as @earth_day_01d using
 region="d"

GMT can automatically uses -JQ15c for `@earth_day_01d`, but not the xarray.DataArray, so need to set `-Rd` here.
---
 pygmt/tests/test_grdimage_image.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pygmt/tests/test_grdimage_image.py b/pygmt/tests/test_grdimage_image.py
index e2ec8980d18..426658a17e1 100644
--- a/pygmt/tests/test_grdimage_image.py
+++ b/pygmt/tests/test_grdimage_image.py
@@ -37,7 +37,7 @@ def test_grdimage_image_dataarray(xr_image):
     Plot a 3-band RGB image using xarray.DataArray input.
     """
     fig = Figure()
-    fig.grdimage(grid=xr_image)
+    fig.grdimage(grid=xr_image, region="d")
     return fig
 
 

From f13d12be5c064fb3c1b7423d495d96c0bad7fd9d Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Mon, 30 Sep 2024 12:54:38 +1300
Subject: [PATCH 23/31] Delete comment about images not being supported

Co-authored-by: Dongdong Tian <seisman.info@gmail.com>
---
 pygmt/datasets/load_remote_dataset.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py
index 052b9734f2a..c79e9936fef 100644
--- a/pygmt/datasets/load_remote_dataset.py
+++ b/pygmt/datasets/load_remote_dataset.py
@@ -428,7 +428,6 @@ 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 = {"R": region}  # region can be None
     if name in {"earth_day"}:
         kind = "image"

From 471c6e8b2364dbcbf9b3419133a62e8a6525796b Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Mon, 30 Sep 2024 13:03:10 +1300
Subject: [PATCH 24/31] Set OGC:CRS84 projection on blue_marble image

Set the projection system so that the xarray.DataArray image is plotted as a Geographic image rather than a Cartesian image by GMT.
---
 pygmt/datasets/earth_daynight.py   | 3 +++
 pygmt/tests/test_grdimage_image.py | 2 +-
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/pygmt/datasets/earth_daynight.py b/pygmt/datasets/earth_daynight.py
index fcb78ad445b..76bb3564e9e 100644
--- a/pygmt/datasets/earth_daynight.py
+++ b/pygmt/datasets/earth_daynight.py
@@ -82,4 +82,7 @@ def load_blue_marble(
         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/tests/test_grdimage_image.py b/pygmt/tests/test_grdimage_image.py
index 426658a17e1..e2ec8980d18 100644
--- a/pygmt/tests/test_grdimage_image.py
+++ b/pygmt/tests/test_grdimage_image.py
@@ -37,7 +37,7 @@ def test_grdimage_image_dataarray(xr_image):
     Plot a 3-band RGB image using xarray.DataArray input.
     """
     fig = Figure()
-    fig.grdimage(grid=xr_image, region="d")
+    fig.grdimage(grid=xr_image)
     return fig
 
 

From e886408f993a89d1ba9700fdcd9a41007ab878d3 Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Mon, 30 Sep 2024 14:22:17 +1300
Subject: [PATCH 25/31] Add note on registration/gtype

Xref #2384.
---
 pygmt/datasets/earth_daynight.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/pygmt/datasets/earth_daynight.py b/pygmt/datasets/earth_daynight.py
index 76bb3564e9e..55220f90e6d 100644
--- a/pygmt/datasets/earth_daynight.py
+++ b/pygmt/datasets/earth_daynight.py
@@ -68,6 +68,16 @@ def load_blue_marble(
     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 5b6bb42e457aa735ccf5cd2d93ea46d76b1fc321 Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Mon, 30 Sep 2024 14:25:39 +1300
Subject: [PATCH 26/31] Rename earth_daynight to earth_day

---
 pygmt/datasets/__init__.py                                      | 2 +-
 pygmt/datasets/{earth_daynight.py => earth_day.py}              | 0
 ...st_datasets_earth_daynight.py => test_datasets_earth_day.py} | 2 +-
 3 files changed, 2 insertions(+), 2 deletions(-)
 rename pygmt/datasets/{earth_daynight.py => earth_day.py} (100%)
 rename pygmt/tests/{test_datasets_earth_daynight.py => test_datasets_earth_day.py} (96%)

diff --git a/pygmt/datasets/__init__.py b/pygmt/datasets/__init__.py
index fe0636245fd..3aeb56d44fd 100644
--- a/pygmt/datasets/__init__.py
+++ b/pygmt/datasets/__init__.py
@@ -5,7 +5,7 @@
 """
 
 from pygmt.datasets.earth_age import load_earth_age
-from pygmt.datasets.earth_daynight import load_blue_marble
+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_daynight.py b/pygmt/datasets/earth_day.py
similarity index 100%
rename from pygmt/datasets/earth_daynight.py
rename to pygmt/datasets/earth_day.py
diff --git a/pygmt/tests/test_datasets_earth_daynight.py b/pygmt/tests/test_datasets_earth_day.py
similarity index 96%
rename from pygmt/tests/test_datasets_earth_daynight.py
rename to pygmt/tests/test_datasets_earth_day.py
index 215ff3390d8..b38d5d51396 100644
--- a/pygmt/tests/test_datasets_earth_daynight.py
+++ b/pygmt/tests/test_datasets_earth_day.py
@@ -1,5 +1,5 @@
 """
-Test basic functionality for loading Blue and Black Marble datasets.
+Test basic functionality for loading Blue Marble datasets.
 """
 
 import numpy as np

From 755275e215bdc84e230731bcb018586d903893de Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Mon, 30 Sep 2024 14:36:16 +1300
Subject: [PATCH 27/31] Remove test_blue_marble_01m_default_registration

The earth_day images are always pixel-registered anyway.
---
 pygmt/tests/test_datasets_earth_day.py | 18 ------------------
 1 file changed, 18 deletions(-)

diff --git a/pygmt/tests/test_datasets_earth_day.py b/pygmt/tests/test_datasets_earth_day.py
index b38d5d51396..2bcf226f72e 100644
--- a/pygmt/tests/test_datasets_earth_day.py
+++ b/pygmt/tests/test_datasets_earth_day.py
@@ -42,21 +42,3 @@ def test_blue_marble_01d_with_region():
     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)
-
-
-def test_blue_marble_01m_default_registration():
-    """
-    Test that the grid returned by default for the 1 arc-minute resolution has
-    a "pixel" registration.
-    """
-    data = load_blue_marble(resolution="01m", region=[-10, -9, 3, 5])
-    assert data.shape == (3, 120, 60)
-    assert data.dtype == "uint8"
-    assert data.gmt.registration == 1
-    assert data.gmt.gtype == 1
-    npt.assert_allclose(data.coords["y"].data.min(), 3.0083333333333333)
-    npt.assert_allclose(data.coords["y"].data.max(), 4.991666666666666)
-    npt.assert_allclose(data.coords["x"].data.min(), -9.991666666666667)
-    npt.assert_allclose(data.coords["x"].data.max(), -9.008333333333335)
-    npt.assert_allclose(data.min(), 10, atol=1)
-    npt.assert_allclose(data.max(), 79, atol=1)

From f7e381e2f481d4683de027b2801b3dee245c5688 Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Mon, 30 Sep 2024 14:39:17 +1300
Subject: [PATCH 28/31] Remove mention of tiled images for 5 arc-min or higher

The earth_day images are not tiled, but come as a single GeoTIFF file for all resolutions.
---
 pygmt/datasets/earth_day.py           |  9 ++++-----
 pygmt/datasets/load_remote_dataset.py | 12 ++++++------
 2 files changed, 10 insertions(+), 11 deletions(-)

diff --git a/pygmt/datasets/earth_day.py b/pygmt/datasets/earth_day.py
index 55220f90e6d..2e62d2a1b25 100644
--- a/pygmt/datasets/earth_day.py
+++ b/pygmt/datasets/earth_day.py
@@ -1,8 +1,8 @@
 """
-Function to download the NASA Blue Marble image datasets from the GMT data
-server, and load as :class:`xarray.DataArray`.
+Function to download the NASA Blue Marble image datasets from the GMT data server, and
+load as :class:`xarray.DataArray`.
 
-The grids are available in various resolutions.
+The images are available in various resolutions.
 """
 
 from collections.abc import Sequence
@@ -60,8 +60,7 @@ def load_blue_marble(
 
     region
         The subregion of the image to load, in the form of a sequence [*xmin*, *xmax*,
-        *ymin*, *ymax*]. Required for images with resolutions higher than 5 arc-minutes
-        (i.e., ``"05m"``).
+        *ymin*, *ymax*].
 
     Returns
     -------
diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py
index c79e9936fef..d27a92950ca 100644
--- a/pygmt/datasets/load_remote_dataset.py
+++ b/pygmt/datasets/load_remote_dataset.py
@@ -89,12 +89,12 @@ class GMTRemoteDataset(NamedTuple):
             "15m": Resolution("15m", registrations=["pixel"]),
             "10m": Resolution("10m", registrations=["pixel"]),
             "06m": Resolution("06m", registrations=["pixel"]),
-            "05m": Resolution("05m", registrations=["pixel"], tiled=True),
-            "04m": Resolution("04m", registrations=["pixel"], tiled=True),
-            "03m": Resolution("03m", registrations=["pixel"], tiled=True),
-            "02m": Resolution("02m", registrations=["pixel"], tiled=True),
-            "01m": Resolution("01m", registrations=["pixel"], tiled=True),
-            "30s": Resolution("30s", registrations=["pixel"], tiled=True),
+            "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(

From 9a5e36bde20d4648196fc3bf5d0e4d37cbc3968b Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Mon, 30 Sep 2024 14:58:07 +1300
Subject: [PATCH 29/31] Apply suggestions from code review

Co-Authored-By: Dongdong Tian <seisman.info@gmail.com>
---
 pygmt/datasets/load_remote_dataset.py  | 8 ++------
 pygmt/tests/test_datasets_earth_day.py | 3 ---
 2 files changed, 2 insertions(+), 9 deletions(-)

diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py
index d27a92950ca..49b28fa6a6f 100644
--- a/pygmt/datasets/load_remote_dataset.py
+++ b/pygmt/datasets/load_remote_dataset.py
@@ -429,12 +429,8 @@ def _load_remote_dataset(
         )
 
     kwdict = {"R": region}  # region can be None
-    if name in {"earth_day"}:
-        kind = "image"
-        kwdict.update({"T": "i"})
-    else:
-        kind = "grid"
-        kwdict.update({"T": "g"})
+    kind = "image" if name in {"earth_day"} else "grid"
+    kwdict["T"] = "i" if kind == "image" else "g"
     with Session() as lib:
         with lib.virtualfile_out(kind=kind) as voutgrd:
             lib.call_module(
diff --git a/pygmt/tests/test_datasets_earth_day.py b/pygmt/tests/test_datasets_earth_day.py
index 2bcf226f72e..f04d7f7f970 100644
--- a/pygmt/tests/test_datasets_earth_day.py
+++ b/pygmt/tests/test_datasets_earth_day.py
@@ -4,11 +4,8 @@
 
 import numpy as np
 import numpy.testing as npt
-import pytest
 from pygmt.datasets import load_blue_marble
 
-rioxarray = pytest.importorskip("rioxarray")
-
 
 def test_blue_marble_01d():
     """

From a0fa29b301e820eeb980d66cbcbfd2073a34c567 Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Mon, 30 Sep 2024 15:14:44 +1300
Subject: [PATCH 30/31] Create kwdict with R and T keys once

Co-authored-by: Dongdong Tian <seisman.info@gmail.com>
---
 pygmt/datasets/load_remote_dataset.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py
index 49b28fa6a6f..5d900e362fd 100644
--- a/pygmt/datasets/load_remote_dataset.py
+++ b/pygmt/datasets/load_remote_dataset.py
@@ -428,9 +428,11 @@ def _load_remote_dataset(
             f"'region' is required for {dataset.description} resolution '{resolution}'."
         )
 
-    kwdict = {"R": region}  # region can be None
     kind = "image" if name in {"earth_day"} else "grid"
-    kwdict["T"] = "i" if kind == "image" else "g"
+    kwdict = {
+        "R": region,  # region can be None
+        "T": "i" if kind == "image" else "g",
+    }  
     with Session() as lib:
         with lib.virtualfile_out(kind=kind) as voutgrd:
             lib.call_module(

From 70f60414e9e1c8e83256ae98da81e685074d59cd Mon Sep 17 00:00:00 2001
From: Wei Ji <23487320+weiji14@users.noreply.github.com>
Date: Mon, 30 Sep 2024 15:15:59 +1300
Subject: [PATCH 31/31] Lint

---
 pygmt/datasets/load_remote_dataset.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py
index 5d900e362fd..e1d12596f86 100644
--- a/pygmt/datasets/load_remote_dataset.py
+++ b/pygmt/datasets/load_remote_dataset.py
@@ -432,7 +432,7 @@ def _load_remote_dataset(
     kwdict = {
         "R": region,  # region can be None
         "T": "i" if kind == "image" else "g",
-    }  
+    }
     with Session() as lib:
         with lib.virtualfile_out(kind=kind) as voutgrd:
             lib.call_module(