diff --git a/pygmt/tests/baseline/test_grdimage_grid_no_redunant_360.png.dvc b/pygmt/tests/baseline/test_grdimage_grid_no_redunant_360.png.dvc
new file mode 100644
index 00000000000..55edc71dbfe
--- /dev/null
+++ b/pygmt/tests/baseline/test_grdimage_grid_no_redunant_360.png.dvc
@@ -0,0 +1,5 @@
+outs:
+- md5: 6145622653eaedd2b4845400aa09ac75
+  size: 239431
+  hash: md5
+  path: test_grdimage_grid_no_redunant_360.png
diff --git a/pygmt/tests/test_grdimage.py b/pygmt/tests/test_grdimage.py
index 8f0e538d270..ee48875240a 100644
--- a/pygmt/tests/test_grdimage.py
+++ b/pygmt/tests/test_grdimage.py
@@ -5,7 +5,9 @@
 import numpy as np
 import pytest
 import xarray as xr
+from packaging.version import Version
 from pygmt import Figure
+from pygmt.clib import __gmt_version__
 from pygmt.datasets import load_earth_relief
 from pygmt.exceptions import GMTInvalidInput
 from pygmt.helpers.testing import check_figures_equal
@@ -252,3 +254,41 @@ def test_grdimage_imgout_fails(grid):
         fig.grdimage(grid, img_out="out.png")
     with pytest.raises(GMTInvalidInput):
         fig.grdimage(grid, A="out.png")
+
+
+@pytest.mark.xfail(
+    condition=Version(__gmt_version__) <= Version("6.5.0"),
+    reason="Upstream bug fixed in https://github.com/GenericMappingTools/gmt/pull/8554",
+)
+@pytest.mark.mpl_image_compare()
+def test_grdimage_grid_no_redunant_360():
+    """
+    Test that global grids with and without redundant 360/0 longitude values work.
+
+    Test for https://github.com/GenericMappingTools/pygmt/issues/3331.
+    """
+    # Global grid [-180, 180, -90, 90] with redundant longitude at 180/-180
+    da1 = load_earth_relief(region=[-180, 180, -90, 90])
+    # Global grid [0, 360, -90, 90] with redundant longitude at 360/0
+    da2 = load_earth_relief(region=[0, 360, -90, 90])
+
+    # Global grid [-180, 179, -90, 90] without redundant longitude at 180/-180
+    da3 = da1[:, 0:360]
+    da3.gmt.registration, da3.gmt.gtype = 0, 1
+    assert da3.shape == (181, 360)
+    assert da3.lon.to_numpy().min() == -180.0
+    assert da3.lon.to_numpy().max() == 179.0
+
+    # Global grid [0, 359, -90, 90] without redundant longitude at 360/0
+    da4 = da2[:, 0:360]
+    da4.gmt.registration, da4.gmt.gtype = 0, 1
+    assert da4.shape == (181, 360)
+    assert da4.lon.to_numpy().min() == 0.0
+    assert da4.lon.to_numpy().max() == 359.0
+
+    fig = Figure()
+    kwdict = {"projection": "W120/10c", "region": "g", "frame": "+tlon=120"}
+    fig.grdimage(da3, **kwdict)
+    fig.shift_origin(xshift="w+2c")
+    fig.grdimage(da4, **kwdict)
+    return fig