diff --git a/gmt/base_plotting.py b/gmt/base_plotting.py
index 36611978a70..e324785e90d 100644
--- a/gmt/base_plotting.py
+++ b/gmt/base_plotting.py
@@ -150,9 +150,7 @@ def grdimage(self, grid, **kwargs):
             if kind == "file":
                 file_context = dummy_context(grid)
             elif kind == "grid":
-                raise NotImplementedError(
-                    "Sorry, DataArray support is not yet functional."
-                )
+                file_context = lib.grid_to_vfile(grid)
             else:
                 raise GMTInvalidInput("Unrecognized data type: {}".format(type(grid)))
             with file_context as fname:
diff --git a/gmt/tests/baseline/test_grdimage.png b/gmt/tests/baseline/test_grdimage.png
new file mode 100644
index 00000000000..8a292a72cc0
Binary files /dev/null and b/gmt/tests/baseline/test_grdimage.png differ
diff --git a/gmt/tests/baseline/test_grdimage_slice.png b/gmt/tests/baseline/test_grdimage_slice.png
new file mode 100644
index 00000000000..f8ab2fa2ff3
Binary files /dev/null and b/gmt/tests/baseline/test_grdimage_slice.png differ
diff --git a/gmt/tests/test_grdimage.py b/gmt/tests/test_grdimage.py
index 7c5417850c5..620341612c2 100644
--- a/gmt/tests/test_grdimage.py
+++ b/gmt/tests/test_grdimage.py
@@ -9,6 +9,24 @@
 from ..datasets import load_earth_relief
 
 
+@pytest.mark.mpl_image_compare
+def test_grdimage():
+    "Plot an image using an xarray grid"
+    grid = load_earth_relief()
+    fig = Figure()
+    fig.grdimage(grid, cmap="earth", projection="W0/6i")
+    return fig
+
+
+@pytest.mark.mpl_image_compare
+def test_grdimage_slice():
+    "Plot an image using an xarray grid that has been sliced"
+    grid = load_earth_relief().sel(lat=slice(-30, 30))
+    fig = Figure()
+    fig.grdimage(grid, cmap="earth", projection="M6i")
+    return fig
+
+
 @pytest.mark.mpl_image_compare
 def test_grdimage_file():
     "Plot an image using file input"
@@ -16,7 +34,7 @@ def test_grdimage_file():
     fig.grdimage(
         "@earth_relief_60m",
         cmap="ocean",
-        region="-180/180/-70/70",
+        region=[-180, 180, -70, 70],
         projection="W0/10i",
         shading=True,
     )
@@ -28,6 +46,3 @@ def test_grdimage_fails():
     fig = Figure()
     with pytest.raises(GMTInvalidInput):
         fig.grdimage(np.arange(20).reshape((4, 5)))
-    grid = load_earth_relief()
-    with pytest.raises(NotImplementedError):
-        fig.grdimage(grid)