Deprecate pygmt.io.load_dataarray, use xarray.load_dataarray instead (#3922)

weiji14 · web-flow · commit ae82fd646df4 · 2025-04-29T16:23:09.000+12:00
Deprecate usage of `pygmt.io.load_dataarray` for loading NetCDF
and other grids, use `xr.load_dataarray(..., engine="gmt")` instead.

* Switch from pygmt.io.load_dataarray to xr.load_dataarray

Specifically with arguments (..., engine="gmt", decode_kind="grid").

* Suggest passing raster_kind="grid" or "image" in FutureWarning message
* Move test_io_load_dataarray benchmark test to test_xarray_backend
* Replace pygmt.io.load_dataarray with xr.load_dataarray in tests
* Replace xr.open_dataarray with xr.load_dataarray in tests
diff --git a/pygmt/datasets/samples.py b/pygmt/datasets/samples.py
@@ -8,7 +8,6 @@
 import pandas as pd
 import xarray as xr
 from pygmt.exceptions import GMTInvalidInput
-from pygmt.io import load_dataarray
 from pygmt.src import which
 
 
@@ -204,7 +203,7 @@ def _load_earth_relief_holes() -> xr.DataArray:
         is in meters.
     """
     fname = which("@earth_relief_20m_holes.grd", download="c")
-    return load_dataarray(fname, engine="netcdf4")
+    return xr.load_dataarray(fname, engine="gmt", raster_kind="grid")
 
 
 class GMTSampleData(NamedTuple):
diff --git a/pygmt/helpers/testing.py b/pygmt/helpers/testing.py
@@ -7,8 +7,8 @@
 import string
 from pathlib import Path
 
+import xarray as xr
 from pygmt.exceptions import GMTImageComparisonFailure
-from pygmt.io import load_dataarray
 from pygmt.src import which
 
 
@@ -154,7 +154,7 @@ def load_static_earth_relief():
         A grid of Earth relief for internal tests.
     """
     fname = which("@static_earth_relief.nc", download="c")
-    return load_dataarray(fname)
+    return xr.load_dataarray(fname, engine="gmt", raster_kind="grid")
 
 
 def skip_if_no(package):
diff --git a/pygmt/io.py b/pygmt/io.py
@@ -2,9 +2,12 @@
 PyGMT input/output (I/O) utilities.
 """
 
+import warnings
+
 import xarray as xr
 
 
+# TODO(PyGMT>=0.20.0): Remove pygmt.io.load_dataarray
 def load_dataarray(filename_or_obj, **kwargs):
     """
     Open, load into memory, and close a DataArray from a file or file-like object
@@ -19,6 +22,12 @@ def load_dataarray(filename_or_obj, **kwargs):
     to :py:func:`xarray.open_dataarray`. See that documentation for further
     details.
 
+    .. deprecated:: v0.16.0
+       The 'pygmt.io.load_dataarray' function will be removed in v0.20.0. Please use
+       `xarray.load_dataarray(..., engine='gmt', raster_kind='grid')` instead if you
+       were reading grids using the engine='netcdf'; otherwise use `raster_kind='image'`
+       if you were reading multi-band images using engine='rasterio'.
+
     Parameters
     ----------
     filename_or_obj : str or pathlib.Path or file-like or DataStore
@@ -37,6 +46,15 @@ def load_dataarray(filename_or_obj, **kwargs):
     --------
     xarray.open_dataarray
     """
+    msg = (
+        "The 'pygmt.io.load_dataarray' function will be removed in v0.20.0. Please use "
+        "`xarray.load_dataarray(..., engine='gmt', raster_kind='grid')` instead if you "
+        "were reading grids using the engine='netcdf'; otherwise use "
+        "`raster_kind='image'` if you were reading multi-band images using "
+        "engine='rasterio'."
+    )
+    warnings.warn(message=msg, category=FutureWarning, stacklevel=1)
+
     if "cache" in kwargs:
         msg = "'cache' has no effect in this context."
         raise TypeError(msg)
diff --git a/pygmt/tests/test_clib_put_matrix.py b/pygmt/tests/test_clib_put_matrix.py
@@ -111,12 +111,14 @@ def test_put_matrix_grid(dtypes):
                     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])
-                    )
+                dataarray = xr.load_dataarray(
+                    tmp_grid.name, engine="gmt", raster_kind="grid"
+                )
+                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])
+                )
diff --git a/pygmt/tests/test_clib_read_data.py b/pygmt/tests/test_clib_read_data.py
@@ -11,7 +11,6 @@
 from pygmt.clib import Session
 from pygmt.exceptions import GMTCLibError
 from pygmt.helpers import GMTTempFile
-from pygmt.io import load_dataarray
 from pygmt.src import which
 
 try:
@@ -27,7 +26,9 @@ def fixture_expected_xrgrid():
     """
     The expected xr.DataArray object for the static_earth_relief.nc file.
     """
-    return load_dataarray(which("@static_earth_relief.nc"))
+    return xr.load_dataarray(
+        "@static_earth_relief.nc", engine="gmt", raster_kind="grid"
+    )
 
 
 @pytest.fixture(scope="module", name="expected_xrimage")
diff --git a/pygmt/tests/test_dimfilter.py b/pygmt/tests/test_dimfilter.py
@@ -6,7 +6,7 @@
 
 import pytest
 import xarray as xr
-from pygmt import dimfilter, load_dataarray
+from pygmt import dimfilter
 from pygmt.enums import GridRegistration, GridType
 from pygmt.exceptions import GMTInvalidInput
 from pygmt.helpers import GMTTempFile
@@ -57,7 +57,7 @@ def test_dimfilter_outgrid(grid, expected_grid):
         )
         assert result is None  # return value is None
         assert Path(tmpfile.name).stat().st_size > 0  # check that outgrid exists
-        temp_grid = load_dataarray(tmpfile.name)
+        temp_grid = xr.load_dataarray(tmpfile.name, engine="gmt", raster_kind="grid")
         xr.testing.assert_allclose(a=temp_grid, b=expected_grid)
 
 
diff --git a/pygmt/tests/test_grdclip.py b/pygmt/tests/test_grdclip.py
@@ -8,7 +8,7 @@
 import numpy.testing as npt
 import pytest
 import xarray as xr
-from pygmt import grdclip, load_dataarray
+from pygmt import grdclip
 from pygmt.datasets import load_earth_mask
 from pygmt.enums import GridRegistration, GridType
 from pygmt.exceptions import GMTInvalidInput
@@ -54,7 +54,7 @@ def test_grdclip_outgrid(grid, expected_grid):
         )
         assert result is None  # return value is None
         assert Path(tmpfile.name).stat().st_size > 0  # check that outgrid exists
-        temp_grid = load_dataarray(tmpfile.name)
+        temp_grid = xr.load_dataarray(tmpfile.name, engine="gmt", raster_kind="grid")
         assert temp_grid.dims == ("lat", "lon")
         assert temp_grid.gmt.gtype == GridType.GEOGRAPHIC
         assert temp_grid.gmt.registration == GridRegistration.PIXEL
diff --git a/pygmt/tests/test_grdcut.py b/pygmt/tests/test_grdcut.py
@@ -5,7 +5,7 @@
 import numpy as np
 import pytest
 import xarray as xr
-from pygmt import grdcut, load_dataarray
+from pygmt import grdcut
 from pygmt.exceptions import GMTInvalidInput
 from pygmt.helpers import GMTTempFile
 from pygmt.helpers.testing import load_static_earth_relief
@@ -50,7 +50,7 @@ def test_grdcut_dataarray_in_file_out(grid, expected_grid, region):
     with GMTTempFile(suffix=".nc") as tmpfile:
         result = grdcut(grid, outgrid=tmpfile.name, region=region)
         assert result is None  # grdcut returns None if output to a file
-        temp_grid = load_dataarray(tmpfile.name)
+        temp_grid = xr.load_dataarray(tmpfile.name, engine="gmt", raster_kind="grid")
         xr.testing.assert_allclose(a=temp_grid, b=expected_grid)
 
 
diff --git a/pygmt/tests/test_grdfill.py b/pygmt/tests/test_grdfill.py
@@ -8,7 +8,7 @@
 import numpy.testing as npt
 import pytest
 import xarray as xr
-from pygmt import grdfill, load_dataarray
+from pygmt import grdfill
 from pygmt.enums import GridRegistration, GridType
 from pygmt.exceptions import GMTInvalidInput
 from pygmt.helpers import GMTTempFile
@@ -96,7 +96,7 @@ def test_grdfill_file_out(grid, expected_grid):
         result = grdfill(grid=grid, constantfill=20, outgrid=tmpfile.name)
         assert result is None  # return value is None
         assert Path(tmpfile.name).stat().st_size > 0  # check that outfile exists
-        temp_grid = load_dataarray(tmpfile.name)
+        temp_grid = xr.load_dataarray(tmpfile.name, engine="gmt", raster_kind="grid")
         xr.testing.assert_allclose(a=temp_grid, b=expected_grid)
 
 
diff --git a/pygmt/tests/test_grdfilter.py b/pygmt/tests/test_grdfilter.py
@@ -7,7 +7,7 @@
 import numpy as np
 import pytest
 import xarray as xr
-from pygmt import grdfilter, load_dataarray
+from pygmt import grdfilter
 from pygmt.enums import GridRegistration, GridType
 from pygmt.exceptions import GMTInvalidInput
 from pygmt.helpers import GMTTempFile
@@ -71,7 +71,7 @@ def test_grdfilter_dataarray_in_file_out(grid, expected_grid):
         )
         assert result is None  # return value is None
         assert Path(tmpfile.name).stat().st_size > 0  # check that outgrid exists
-        temp_grid = load_dataarray(tmpfile.name)
+        temp_grid = xr.load_dataarray(tmpfile.name, engine="gmt", raster_kind="grid")
         xr.testing.assert_allclose(a=temp_grid, b=expected_grid)
 
 
diff --git a/pygmt/tests/test_grdgradient.py b/pygmt/tests/test_grdgradient.py
@@ -6,7 +6,7 @@
 
 import pytest
 import xarray as xr
-from pygmt import grdgradient, load_dataarray
+from pygmt import grdgradient
 from pygmt.enums import GridRegistration, GridType
 from pygmt.exceptions import GMTInvalidInput
 from pygmt.helpers import GMTTempFile
@@ -50,7 +50,7 @@ def test_grdgradient_outgrid(grid, expected_grid):
         )
         assert result is None  # return value is None
         assert Path(tmpfile.name).stat().st_size > 0  # check that outgrid exists
-        temp_grid = load_dataarray(tmpfile.name)
+        temp_grid = xr.load_dataarray(tmpfile.name, engine="gmt", raster_kind="grid")
         xr.testing.assert_allclose(a=temp_grid, b=expected_grid)
 
 
diff --git a/pygmt/tests/test_grdhisteq.py b/pygmt/tests/test_grdhisteq.py
@@ -8,7 +8,7 @@
 import pandas as pd
 import pytest
 import xarray as xr
-from pygmt import grdhisteq, load_dataarray
+from pygmt import grdhisteq
 from pygmt.enums import GridRegistration, GridType
 from pygmt.exceptions import GMTInvalidInput
 from pygmt.helpers import GMTTempFile
@@ -67,7 +67,7 @@ def test_equalize_grid_outgrid_file(grid, expected_grid, region):
         )
         assert result is None  # return value is None
         assert Path(tmpfile.name).stat().st_size > 0  # check that outgrid exists
-        temp_grid = load_dataarray(tmpfile.name)
+        temp_grid = xr.load_dataarray(tmpfile.name, engine="gmt", raster_kind="grid")
         xr.testing.assert_allclose(a=temp_grid, b=expected_grid)
 
 
diff --git a/pygmt/tests/test_grdlandmask.py b/pygmt/tests/test_grdlandmask.py
@@ -6,7 +6,7 @@
 
 import pytest
 import xarray as xr
-from pygmt import grdlandmask, load_dataarray
+from pygmt import grdlandmask
 from pygmt.enums import GridRegistration, GridType
 from pygmt.exceptions import GMTInvalidInput
 from pygmt.helpers import GMTTempFile
@@ -42,7 +42,7 @@ def test_grdlandmask_outgrid(expected_grid):
         result = grdlandmask(outgrid=tmpfile.name, spacing=1, region=[125, 130, 30, 35])
         assert result is None  # return value is None
         assert Path(tmpfile.name).stat().st_size > 0  # check that outgrid exists
-        temp_grid = load_dataarray(tmpfile.name)
+        temp_grid = xr.load_dataarray(tmpfile.name, engine="gmt", raster_kind="grid")
         xr.testing.assert_allclose(a=temp_grid, b=expected_grid)
 
 
diff --git a/pygmt/tests/test_grdproject.py b/pygmt/tests/test_grdproject.py
@@ -6,7 +6,7 @@
 
 import pytest
 import xarray as xr
-from pygmt import grdproject, load_dataarray
+from pygmt import grdproject
 from pygmt.enums import GridRegistration, GridType
 from pygmt.exceptions import GMTInvalidInput
 from pygmt.helpers import GMTTempFile
@@ -56,7 +56,7 @@ def test_grdproject_file_out(grid, expected_grid):
         )
         assert result is None  # return value is None
         assert Path(tmpfile.name).stat().st_size > 0  # check that outgrid exists
-        temp_grid = load_dataarray(tmpfile.name)
+        temp_grid = xr.load_dataarray(tmpfile.name, engine="gmt", raster_kind="grid")
         xr.testing.assert_allclose(a=temp_grid, b=expected_grid)
 
 
diff --git a/pygmt/tests/test_grdsample.py b/pygmt/tests/test_grdsample.py
@@ -6,7 +6,7 @@
 
 import pytest
 import xarray as xr
-from pygmt import grdsample, load_dataarray
+from pygmt import grdsample
 from pygmt.enums import GridRegistration, GridType
 from pygmt.helpers import GMTTempFile
 from pygmt.helpers.testing import load_static_earth_relief
@@ -67,7 +67,7 @@ def test_grdsample_file_out(grid, expected_grid, region, spacing):
         )
         assert result is None  # return value is None
         assert Path(tmpfile.name).stat().st_size > 0  # check that outgrid exists
-        temp_grid = load_dataarray(tmpfile.name)
+        temp_grid = xr.load_dataarray(tmpfile.name, engine="gmt", raster_kind="grid")
         xr.testing.assert_allclose(a=temp_grid, b=expected_grid)
 
 
diff --git a/pygmt/tests/test_io.py b/pygmt/tests/test_io.py
@@ -10,7 +10,7 @@
 from pygmt.io import load_dataarray
 
 
-@pytest.mark.benchmark
+# TODO(PyGMT>=0.20.0): Remove test_io_load_dataarray
 def test_io_load_dataarray():
     """
     Check that load_dataarray works to read a netCDF grid with GMTDataArrayAccessor
@@ -22,7 +22,10 @@ def test_io_load_dataarray():
             data=rng.random((2, 2)), coords=[[0.1, 0.2], [0.3, 0.4]], dims=("x", "y")
         )
         grid.to_netcdf(tmpfile.name)
-        dataarray = load_dataarray(tmpfile.name)
+
+        with pytest.warns(FutureWarning):
+            dataarray = load_dataarray(tmpfile.name)
+
         assert dataarray.gmt.gtype == GridType.CARTESIAN
         assert dataarray.gmt.registration == GridRegistration.PIXEL
         # this would fail if we used xr.open_dataarray instead of load_dataarray
diff --git a/pygmt/tests/test_nearneighbor.py b/pygmt/tests/test_nearneighbor.py
@@ -82,7 +82,7 @@ def test_nearneighbor_with_outgrid_param(ship_data):
         )
         assert output is None  # check that output is None since outgrid is set
         assert Path(tmpfile.name).stat().st_size > 0  # check that outgrid exists
-        with xr.open_dataarray(tmpfile.name) as grid:
-            assert isinstance(grid, xr.DataArray)  # ensure netCDF grid loads ok
-            assert grid.shape == (121, 121)
-            npt.assert_allclose(grid.mean(), -2378.2385)
+        grid = xr.load_dataarray(tmpfile.name, engine="gmt", raster_kind="grid")
+        assert isinstance(grid, xr.DataArray)  # ensure netCDF grid loads ok
+        assert grid.shape == (121, 121)
+        npt.assert_allclose(grid.mean(), -2378.2385)
diff --git a/pygmt/tests/test_surface.py b/pygmt/tests/test_surface.py
@@ -145,5 +145,5 @@ def test_surface_with_outgrid_param(data, region, spacing, expected_grid):
         )
         assert output is None  # check that output is None since outgrid is set
         assert Path(tmpfile.name).stat().st_size > 0  # check that outgrid exists
-        with xr.open_dataarray(tmpfile.name) as grid:
-            check_values(grid, expected_grid)
+        grid = xr.load_dataarray(tmpfile.name, engine="gmt", raster_kind="grid")
+        check_values(grid, expected_grid)
diff --git a/pygmt/tests/test_triangulate.py b/pygmt/tests/test_triangulate.py
@@ -156,8 +156,8 @@ def test_regular_grid_with_outgrid_param(dataframe, expected_grid):
         )
         assert output is None  # check that output is None since outgrid is set
         assert Path(tmpfile.name).stat().st_size > 0  # check that outgrid exists
-        with xr.open_dataarray(tmpfile.name) as grid:
-            assert isinstance(grid, xr.DataArray)
-            assert grid.gmt.registration == GridRegistration.GRIDLINE
-            assert grid.gmt.gtype == GridType.CARTESIAN
-            xr.testing.assert_allclose(a=grid, b=expected_grid)
+        grid = xr.load_dataarray(tmpfile.name, engine="gmt", raster_kind="grid")
+        assert isinstance(grid, xr.DataArray)
+        assert grid.gmt.registration == GridRegistration.GRIDLINE
+        assert grid.gmt.gtype == GridType.CARTESIAN
+        xr.testing.assert_allclose(a=grid, b=expected_grid)
diff --git a/pygmt/tests/test_xarray_accessor.py b/pygmt/tests/test_xarray_accessor.py
@@ -20,8 +20,7 @@ def test_xarray_accessor_gridline_cartesian():
     Check that the accessor returns the correct registration and gtype values for a
     Cartesian, gridline-registered grid.
     """
-    fname = which(fname="@test.dat.nc", download="a")
-    grid = xr.open_dataarray(fname, engine="netcdf4")
+    grid = xr.load_dataarray("@test.dat.nc", engine="gmt", raster_kind="grid")
     assert grid.gmt.registration == GridRegistration.GRIDLINE
     assert grid.gmt.gtype == GridType.CARTESIAN
 
@@ -31,8 +30,7 @@ def test_xarray_accessor_pixel_geographic():
     Check that the accessor returns the correct registration and gtype values for a
     geographic, pixel-registered grid.
     """
-    fname = which(fname="@earth_relief_01d_p", download="a")
-    grid = xr.open_dataarray(fname, engine="netcdf4")
+    grid = xr.load_dataarray("@earth_relief_01d_p", engine="gmt", raster_kind="grid")
     assert grid.gmt.registration == GridRegistration.PIXEL
     assert grid.gmt.gtype == GridType.GEOGRAPHIC
 
diff --git a/pygmt/tests/test_xarray_backend.py b/pygmt/tests/test_xarray_backend.py
diff --git a/pygmt/tests/test_xyz2grd.py b/pygmt/tests/test_xyz2grd.py
