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)