diff --git a/pygmt/tests/test_clib_put_matrix.py b/pygmt/tests/test_clib_put_matrix.py
index 2dede298fbc..c4ec764396e 100644
--- a/pygmt/tests/test_clib_put_matrix.py
+++ b/pygmt/tests/test_clib_put_matrix.py
@@ -4,6 +4,7 @@
 import numpy as np
 import numpy.testing as npt
 import pytest
+import xarray as xr
 
 
 from .test_clib import mock
@@ -55,7 +56,7 @@ def test_put_matrix_fails():
 
 
 def test_put_matrix_grid():
-    "Check that assigning a numpy 2d array to a grid works"
+    "Check that assigning a numpy 2d array to an ASCII and NetCDF grid works"
     dtypes = "float32 float64 int32 int64 uint32 uint64".split()
     wesn = [10, 15, 30, 40, 0, 0]
     inc = [1, 1]
@@ -85,3 +86,23 @@ def test_put_matrix_grid():
                 # Load the data and check that it's correct
                 newdata = tmp_file.loadtxt(dtype=dtype)
                 npt.assert_allclose(newdata, data)
+
+            # Save the data to a netCDF grid and check that xarray can load it
+            with GMTTempFile() as tmp_grid:
+                lib.write_data(
+                    "GMT_IS_MATRIX",
+                    "GMT_IS_SURFACE",
+                    "GMT_CONTAINER_AND_DATA",
+                    wesn,
+                    tmp_grid.name,
+                    grid,
+                )
+                with xr.open_dataarray(tmp_grid.name) as dataarray:
+                    assert dataarray.shape == shape
+                    npt.assert_allclose(dataarray.data, np.flipud(data))
+                    npt.assert_allclose(
+                        dataarray.coords["x"].actual_range, np.array(wesn[0:2])
+                    )
+                    npt.assert_allclose(
+                        dataarray.coords["y"].actual_range, np.array(wesn[2:4])
+                    )