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