**BREAKING** pygmt.grdcut: Refactor to store output in virtualfiles for grids

pygmt/clib/session.py

@@ -1942,6 +1942,7 @@ def virtualfile_out(
                 "grid": ("GMT_IS_GRID", "GMT_IS_SURFACE"),
                 "image": ("GMT_IS_IMAGE", "GMT_IS_SURFACE"),
             }[kind]
+            # For unknown reasons, 'GMT_OUT' crashes for 'image' kind.
             direction = "GMT_OUT|GMT_IS_REFERENCE" if kind == "image" else "GMT_OUT"
             with self.open_virtualfile(family, geometry, direction, None) as vfile:
                 yield vfile
@@ -2311,3 +2312,35 @@ def extract_region(self) -> np.ndarray:
         if status != 0:
             raise GMTCLibError("Failed to extract region from current figure.")
         return region
+
+
+def _raster_kind(raster: str) -> Literal["grid", "image"]:
+    """
+    Determine the raster kind.
+
+    Examples
+    --------
+    >>> _raster_kind("@earth_relief_01d")
+    'grid'
+    >>> _raster_kind("@static_earth_relief.nc")
+    'grid'
+    >>> _raster_kind("@earth_day_01d")
+    'image'
+    >>> _raster_kind("@hotspots.txt")
+    'grid'
+    """
+    # The logic here is because: an image can be read into a grid container, but a grid
+    # can't be read into an image container. So, try to read the file as an image first.
+    # If fails, try to read it as a grid.
+    with Session() as lib:
+        try:
+            img = lib.read_data(infile=raster, kind="image", mode="GMT_CONTAINER_ONLY")
+            return "image" if img.contents.header.contents.n_bands == 3 else "grid"
+        except GMTCLibError:
+            pass
+        try:
+            _ = lib.read_data(infile=raster, kind="grid", mode="GMT_CONTAINER_ONLY")
+            return "grid"
+        except GMTCLibError:
+            pass
+    return "grid"  # Fallback to "grid" and let GMT determine the type.

pygmt/src/grdcut.py

@@ -4,21 +4,21 @@
 
 import xarray as xr
 from pygmt.clib import Session
+from pygmt.clib.session import _raster_kind
+from pygmt.exceptions import GMTInvalidInput
 from pygmt.helpers import (
-    GMTTempFile,
     build_arg_list,
+    data_kind,
     fmt_docstring,
     kwargs_to_strings,
     use_alias,
 )
-from pygmt.io import load_dataarray
 
 __doctest_skip__ = ["grdcut"]
 
 
 @fmt_docstring
 @use_alias(
-    G="outgrid",
     R="region",
     J="projection",
     N="extend",
@@ -28,9 +28,9 @@
     f="coltypes",
 )
 @kwargs_to_strings(R="sequence")
-def grdcut(grid, **kwargs) -> xr.DataArray | None:
+def grdcut(grid, outgrid: str | None = None, **kwargs) -> xr.DataArray | None:
     r"""
-    Extract subregion from a grid.
+    Extract subregion from a grid or image.
 
     Produce a new ``outgrid`` file which is a subregion of ``grid``. The
     subregion is specified with ``region``; the specified range must not exceed
@@ -100,13 +100,24 @@
     >>> # 12° E to 15° E and a latitude range of 21° N to 24° N
     >>> new_grid = pygmt.grdcut(grid=grid, region=[12, 15, 21, 24])
     """
-    with GMTTempFile(suffix=".nc") as tmpfile:
-        with Session() as lib:
-            with lib.virtualfile_in(check_kind="raster", data=grid) as vingrd:
-                if (outgrid := kwargs.get("G")) is None:
-                    kwargs["G"] = outgrid = tmpfile.name  # output to tmpfile
-                lib.call_module(
-                    module="grdcut", args=build_arg_list(kwargs, infile=vingrd)
-                )
+    # Determine the output data kind based on the input data kind.
+    inkind = data_kind(grid)
+    match inkind:
+        case "image" | "grid":
+            outkind = inkind
+        case "file":
+            outkind = _raster_kind(grid)
+        case "_":
+            msg = f"Unsupported data type {type(grid)}."
+            raise GMTInvalidInput(msg)
 
-        return load_dataarray(outgrid) if outgrid == tmpfile.name else None
+    with Session() as lib:
+        with (
+            lib.virtualfile_in(check_kind="raster", data=grid) as vingrd,
+            lib.virtualfile_out(kind=outkind, fname=outgrid) as voutgrd,
+        ):
+            kwargs["G"] = voutgrd
+            lib.call_module(module="grdcut", args=build_arg_list(kwargs, infile=vingrd))
+            return lib.virtualfile_to_raster(
+                vfname=voutgrd, kind=outkind, outgrid=outgrid
+            )

pygmt/tests/test_grdcut_image.py

@@ -0,0 +1,85 @@
+"""
+Test pygmt.grdcut on images.
+"""
+
+from pathlib import Path
+
+import numpy as np
+import pytest
+import xarray as xr
+from pygmt import grdcut, which
+from pygmt.helpers import GMTTempFile
+
+try:
+    import rioxarray
+
+    _HAS_RIOXARRAY = True
+except ImportError:
+    _HAS_RIOXARRAY = False
+
+
+@pytest.fixture(scope="module", name="region")
+def fixture_region():
+    """
+    Set the data region.
+    """
+    return [-53, -49, -20, -17]
+
+
+@pytest.fixture(scope="module", name="expected_image")
+def fixture_expected_image():
+    """
+    Load the expected grdcut image result.
+    """
+    return xr.DataArray(
+        data=np.array(
+            [
+                [[90, 93, 95, 90], [91, 90, 91, 91], [91, 90, 89, 90]],
+                [[87, 88, 88, 89], [88, 87, 86, 85], [90, 90, 89, 88]],
+                [[48, 49, 49, 45], [49, 48, 47, 45], [48, 47, 48, 46]],
+            ],
+            dtype=np.uint8,
+        ),
+        coords={
+            "band": [1, 2, 3],
+            "x": [-52.5, -51.5, -50.5, -49.5],
+            "y": [-17.5, -18.5, -19.5],
+        },
+        dims=["band", "y", "x"],
+        attrs={
+            "scale_factor": 1.0,
+            "add_offset": 0.0,
+        },
+    )
+
+
+@pytest.mark.benchmark
+def test_grdcut_image_file(region, expected_image):
+    """
+    Test grdcut on an input image file.
+    """
+    result = grdcut("@earth_day_01d_p", region=region)
+    xr.testing.assert_allclose(a=result, b=expected_image)
+
+
+@pytest.mark.skipif(not _HAS_RIOXARRAY, reason="rioxarray is not installed")
+def test_grdcut_image_dataarray(region, expected_image):
+    """
+    Test grdcut on an input xarray.DataArray object.
+    """
+    raster = rioxarray.open_rasterio(which("@earth_day_01d", download="a")).load()
+    result = grdcut(raster, region=region)
+    xr.testing.assert_allclose(a=result, b=expected_image)
+
+
+def test_grdcut_image_file_in_file_out(region, expected_image):
+    """
+    Test grdcut on an input image file and outputs to another image file.
+    """
+    with GMTTempFile(suffix=".tif") as tmp:
+        result = grdcut("@earth_day_01d_p", region=region, outgrid=tmp.name)
+        assert result is None
+        assert Path(tmp.name).stat().st_size > 0
+    if _HAS_RIOXARRAY:
+        raster = rioxarray.open_rasterio(which("@earth_day_01d", download="a")).load()
+        xr.testing.assert_allclose(a=raster, b=expected_image)