Allow GMTDataArrayAccessor to work on sliced datacubes (GenericMappin…

…gTools#1581)

If `grdinfo` can't read the source NetCDF file as it
is an n-dimensional datacube (instead of a 2D grid),
we fallback to setting the default gridline registration
and Cartesian grid type.

* Use with statement to ensure the grid is properly closed

So that the file can be deleted on Windows.

Co-authored-by: Will Schlitzer <[email protected]>

pygmt/accessors.py

@@ -34,7 +34,7 @@ def __init__(self, xarray_obj):
             self._registration, self._gtype = map(
                 int, grdinfo(self._source, per_column="n", o="10,11").split()
             )
-        except KeyError:
+        except (KeyError, ValueError):
             self._registration = 0  # Default to Gridline registration
             self._gtype = 0  # Default to Cartesian grid type
 

pygmt/tests/test_accessor.py

@@ -1,6 +1,8 @@
 """
 Test the behaviour of the GMTDataArrayAccessor class.
 """
+import os
+
 import pytest
 import xarray as xr
 from pygmt import which
@@ -66,3 +68,25 @@ def test_accessor_set_non_boolean():
 
     with pytest.raises(GMTInvalidInput):
         grid.gmt.gtype = 2
+
+
+def test_accessor_sliced_datacube():
+    """
+    Check that a 2D grid which is sliced from an n-dimensional datacube works
+    with accessor methods.
+
+    This is a regression test for
+    https://github.com/GenericMappingTools/pygmt/issues/1578.
+    """
+    try:
+        fname = which(
+            "https://github.com/pydata/xarray-data/raw/master/eraint_uvz.nc",
+            download="u",
+        )
+        with xr.open_dataset(fname) as dataset:
+            grid = dataset.sel(level=500, month=1, drop=True).z
+
+        assert grid.gmt.registration == 0  # gridline registration
+        assert grid.gmt.gtype == 0  # cartesian coordinate type
+    finally:
+        os.remove(fname)