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