Wrap GMT's standard data type GMT_GRID for grids

doc/api/index.rst

@@ -292,6 +292,7 @@ Python objects to and from GMT virtual files:
     clib.Session.virtualfile_in
     clib.Session.virtualfile_out
     clib.Session.virtualfile_to_dataset
+    clib.Session.virtualfile_to_raster
 
 Low level access (these are mostly used by the :mod:`pygmt.clib` package):
 

pygmt/clib/session.py

@@ -14,6 +14,7 @@
 
 import numpy as np
 import pandas as pd
+import xarray as xr
 from packaging.version import Version
 from pygmt.clib.conversion import (
     array_to_datetime,
@@ -1739,7 +1740,9 @@
         return c_inquire_virtualfile(self.session_pointer, vfname.encode())
 
     def read_virtualfile(
-        self, vfname: str, kind: Literal["dataset", "grid", None] = None
+        self,
+        vfname: str,
+        kind: Literal["dataset", "grid", "image", "cube", None] = None,
     ):
         """
         Read data from a virtual file and optionally cast into a GMT data container.
@@ -1798,6 +1801,8 @@
         # _GMT_DATASET).
         if kind is None:  # Return the ctypes void pointer
             return pointer
+        if kind in ["image", "cube"]:
+            raise NotImplementedError(f"kind={kind} is not supported yet.")
         dtype = {"dataset": _GMT_DATASET, "grid": _GMT_GRID}[kind]
         return ctp.cast(pointer, ctp.POINTER(dtype))
 
@@ -1946,6 +1951,70 @@
             return result.to_numpy()
         return result  # pandas.DataFrame output
 
+    def virtualfile_to_raster(
+        self,
+        vfname: str,
+        kind: Literal["grid", "image", "cube", None] = "grid",
+        outgrid: str | None = None,
+    ) -> xr.DataArray | None:
+        """
+        Output raster data stored in a virtual file to an :class:`xarray.DataArray`
+        object.
+
+        The raster data can be a grid, an image or a cube.
+
+        Parameters
+        ----------
+        vfname
+            The virtual file name that stores the result grid/image/cube.
+        kind
+            Type of the raster data. Valid values are ``"grid"``, ``"image"``,
+            ``"cube"`` or ``None``. If ``None``, will inquire the data type from the
+            virtual file name.
+        outgrid
+            Name of the output grid/image/cube. If specified, it means the raster data
+            was already saved into an actual file and will return ``None``.
+
+        Returns
+        -------
+        result
+            The result grid/image/cube. If ``outgrid`` is specified, return ``None``.
+
+        Examples
+        --------
+        >>> from pathlib import Path
+        >>> from pygmt.clib import Session
+        >>> from pygmt.helpers import GMTTempFile
+        >>> with Session() as lib:
+        ...     # file output
+        ...     with GMTTempFile(suffix=".nc") as tmpfile:
+        ...         outgrid = tmpfile.name
+        ...         with lib.virtualfile_out(kind="grid", fname=outgrid) as voutgrd:
+        ...             lib.call_module("read", f"@earth_relief_01d_g {voutgrd} -Tg")
+        ...             result = lib.virtualfile_to_raster(
+        ...                 vfname=voutgrd, outgrid=outgrid
+        ...             )
+        ...             assert result == None
+        ...             assert Path(outgrid).stat().st_size > 0
+        ...
+        ...     # xarray.DataArray output
+        ...     outgrid = None
+        ...     with lib.virtualfile_out(kind="grid", fname=outgrid) as voutgrd:
+        ...         lib.call_module("read", f"@earth_relief_01d_g {voutgrd} -Tg")
+        ...         result = lib.virtualfile_to_raster(vfname=voutgrd, outgrid=outgrid)
+        ...         assert isinstance(result, xr.DataArray)
+        """
+        if outgrid is not None:
+            return None
+        if kind is None:  # Inquire the data family from the virtualfile
+            family = self.inquire_virtualfile(vfname)
+            kind = {  # type: ignore[assignment]
+                self["GMT_IS_GRID"]: "grid",
+                self["GMT_IS_IMAGE"]: "image",
+                self["GMT_IS_CUBE"]: "cube",
+            }[family]
+        return self.read_virtualfile(vfname, kind=kind).contents.to_dataarray()
+
     def extract_region(self):
         """
         Extract the WESN bounding box of the currently active figure.

pygmt/datatypes/grid.py

@@ -3,7 +3,197 @@
 """
 
 import ctypes as ctp
+from typing import ClassVar
+
+import numpy as np
+import xarray as xr
+from pygmt.datatypes.header import _GMT_GRID_HEADER, gmt_grdfloat
 
 
 class _GMT_GRID(ctp.Structure):  # noqa: N801
-    pass
+    """
+    GMT grid structure for holding a grid and its header.
+
+    This class is only meant for internal use and is not exposed to users. See the GMT
+    source code gmt_resources.h for the original C structure definitions.
+
+    Examples
+    --------
+    >>> from pygmt.clib import Session
+    >>> with Session() as lib:
+    ...     with lib.virtualfile_out(kind="grid") as voutgrd:
+    ...         lib.call_module("read", f"@static_earth_relief.nc {voutgrd} -Tg")
+    ...         # Read the grid from the virtual file
+    ...         grid = lib.read_virtualfile(voutgrd, kind="grid").contents
+    ...         # The grid header
+    ...         header = grid.header.contents
+    ...         # Access the header properties
+    ...         print(header.n_rows, header.n_columns, header.registration)
+    ...         print(header.wesn[:], header.z_min, header.z_max, header.inc[:])
+    ...         print(header.z_scale_factor, header.z_add_offset)
+    ...         print(header.x_units, header.y_units, header.z_units)
+    ...         print(header.title)
+    ...         print(header.command)
+    ...         print(header.remark)
+    ...         print(header.nm, header.size, header.complex_mode)
+    ...         print(header.type, header.n_bands, header.mx, header.my)
+    ...         print(header.pad[:])
+    ...         print(header.mem_layout, header.nan_value, header.xy_off)
+    ...         # The x and y coordinates
+    ...         print(grid.x[: header.n_columns])
+    ...         print(grid.y[: header.n_rows])
+    ...         # The data array (with paddings)
+    ...         data = np.reshape(
+    ...             grid.data[: header.mx * header.my], (header.my, header.mx)
+    ...         )
+    ...         # The data array (without paddings)
+    ...         pad = header.pad[:]
+    ...         data = data[pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1]]
+    ...         print(data)
+    14 8 1
+    [-55.0, -47.0, -24.0, -10.0] 190.0 981.0 [1.0, 1.0]
+    1.0 0.0
+    b'longitude [degrees_east]' b'latitude [degrees_north]' b'elevation (m)'
+    b'Produced by grdcut'
+    b'grdcut @earth_relief_01d_p -R-55/-47/-24/-10 -Gstatic_earth_relief.nc'
+    b'Reduced by Gaussian Cartesian filtering (111.2 km fullwidth) from ...'
+    112 216 0
+    18 1 12 18
+    [2, 2, 2, 2]
+    b'' nan 0.5
+    [-54.5, -53.5, -52.5, -51.5, -50.5, -49.5, -48.5, -47.5]
+    [-10.5, -11.5, -12.5, -13.5, -14.5, -15.5, ..., -22.5, -23.5]
+    [[347.5 331.5 309.  282.  190.  208.  299.5 348. ]
+    [349.  313.  325.5 247.  191.  225.  260.  452.5]
+    [345.5 320.  335.  292.  207.5 247.  325.  346.5]
+    [450.5 395.5 366.  248.  250.  354.5 550.  797.5]
+    [494.5 488.5 357.  254.5 286.  484.5 653.5 930. ]
+    [601.  526.5 535.  299.  398.5 645.  797.5 964. ]
+    [308.  595.5 555.5 556.  580.  770.  927.  920. ]
+    [521.5 682.5 796.  886.  571.5 638.5 739.5 881.5]
+    [310.  521.5 757.  570.5 538.5 524.  686.5 794. ]
+    [561.5 539.  446.5 481.5 439.5 553.  726.5 981. ]
+    [557.  435.  385.5 345.5 413.5 496.  519.5 833.5]
+    [373.  367.5 349.  352.5 419.5 428.  570.  667.5]
+    [383.  284.5 344.5 394.  491.  556.5 578.5 618.5]
+    [347.5 344.5 386.  640.5 617.  579.  646.5 671. ]]
+    """
+
+    _fields_: ClassVar = [
+        # Pointer to full GMT header for grid
+        ("header", ctp.POINTER(_GMT_GRID_HEADER)),
+        # Pointer to grid data
+        ("data", ctp.POINTER(gmt_grdfloat)),
+        # Pointer to x coordinate vector
+        ("x", ctp.POINTER(ctp.c_double)),
+        # Pointer to y coordinate vector
+        ("y", ctp.POINTER(ctp.c_double)),
+        # Low-level information for GMT use only
+        ("hidden", ctp.c_void_p),
+    ]
+
+    def to_dataarray(self) -> xr.DataArray:
+        """
+        Convert a _GMT_GRID object to a :class:`xarray.DataArray` object.
+
+        Returns
+        -------
+        dataarray
+            A :class:`xr.DataArray` object.
+
+        Examples
+        --------
+        >>> from pygmt.clib import Session
+        >>> with Session() as lib:
+        ...     with lib.virtualfile_out(kind="grid") as voutgrd:
+        ...         lib.call_module("read", f"@static_earth_relief.nc {voutgrd} -Tg")
+        ...         # Read the grid from the virtual file
+        ...         grid = lib.read_virtualfile(voutgrd, kind="grid")
+        ...         # Convert to xarray.DataArray and use it later
+        ...         da = grid.contents.to_dataarray()
+        >>> da  # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS
+        <xarray.DataArray 'z' (lat: 14, lon: 8)>...
+        array([[347.5, 344.5, 386. , 640.5, 617. , 579. , 646.5, 671. ],
+               [383. , 284.5, 344.5, 394. , 491. , 556.5, 578.5, 618.5],
+               [373. , 367.5, 349. , 352.5, 419.5, 428. , 570. , 667.5],
+               [557. , 435. , 385.5, 345.5, 413.5, 496. , 519.5, 833.5],
+               [561.5, 539. , 446.5, 481.5, 439.5, 553. , 726.5, 981. ],
+               [310. , 521.5, 757. , 570.5, 538.5, 524. , 686.5, 794. ],
+               [521.5, 682.5, 796. , 886. , 571.5, 638.5, 739.5, 881.5],
+               [308. , 595.5, 555.5, 556. , 580. , 770. , 927. , 920. ],
+               [601. , 526.5, 535. , 299. , 398.5, 645. , 797.5, 964. ],
+               [494.5, 488.5, 357. , 254.5, 286. , 484.5, 653.5, 930. ],
+               [450.5, 395.5, 366. , 248. , 250. , 354.5, 550. , 797.5],
+               [345.5, 320. , 335. , 292. , 207.5, 247. , 325. , 346.5],
+               [349. , 313. , 325.5, 247. , 191. , 225. , 260. , 452.5],
+               [347.5, 331.5, 309. , 282. , 190. , 208. , 299.5, 348. ]])
+        Coordinates:
+          * lat      (lat) float64... -23.5 -22.5 -21.5 -20.5 ... -12.5 -11.5 -10.5
+          * lon      (lon) float64... -54.5 -53.5 -52.5 -51.5 -50.5 -49.5 -48.5 -47.5
+        Attributes:
+            Conventions:   CF-1.7
+            title:         Produced by grdcut
+            history:       grdcut @earth_relief_01d_p -R-55/-47/-24/-10 -Gstatic_ea...
+            description:   Reduced by Gaussian Cartesian filtering (111.2 km fullwi...
+            long_name:     elevation (m)
+            actual_range:  [190. 981.]
+        >>> da.coords["lon"]  # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS
+        <xarray.DataArray 'lon' (lon: 8)>...
+        array([-54.5, -53.5, -52.5, -51.5, -50.5, -49.5, -48.5, -47.5])
+        Coordinates:
+          * lon      (lon) float64... -54.5 -53.5 -52.5 -51.5 -50.5 -49.5 -48.5 -47.5
+        Attributes:
+            long_name:      longitude
+            units:          degrees_east
+            standard_name:  longitude
+            axis:           X
+            actual_range:   [-55. -47.]
+        >>> da.coords["lat"]  # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS
+        <xarray.DataArray 'lat' (lat: 14)>...
+        array([-23.5, -22.5, -21.5, -20.5, -19.5, -18.5, -17.5, -16.5, -15.5, -14.5,
+            -13.5, -12.5, -11.5, -10.5])
+        Coordinates:
+          * lat      (lat) float64... -23.5 -22.5 -21.5 -20.5 ... -12.5 -11.5 -10.5
+        Attributes:
+            long_name:      latitude
+            units:          degrees_north
+            standard_name:  latitude
+            axis:           Y
+            actual_range:   [-24. -10.]
+        >>> da.gmt.registration, da.gmt.gtype
+        (1, 1)
+        """
+        # The grid header
+        header = self.header.contents
+
+        # Get dimensions and their attributes from the header.
+        dims, dim_attrs = header.dims, header.dim_attrs
+        # The coordinates, given as a tuple of the form (dims, data, attrs)
+        coords = [
+            (dims[0], self.y[: header.n_rows], dim_attrs[0]),
+            (dims[1], self.x[: header.n_columns], dim_attrs[1]),
+        ]
+
+        # The data array without paddings
+        pad = header.pad[:]
+        data = np.reshape(self.data[: header.mx * header.my], (header.my, header.mx))[
+            pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1]
+        ]
+
+        # Create the xarray.DataArray object
+        grid = xr.DataArray(
+            data, coords=coords, name=header.name, attrs=header.data_attrs
+        )
+
+        # Flip the coordinates and data if necessary so that coordinates are ascending.
+        # `grid.sortby(list(grid.dims))` sometimes causes crashes.
+        # The solution comes from https://github.com/pydata/xarray/discussions/6695.
+        for dim in grid.dims:
+            if grid[dim][0] > grid[dim][1]:
+                grid = grid.isel({dim: slice(None, None, -1)})
+
+        # Set GMT accessors.
+        # Must put at the end, otherwise info gets lost after certain grid operations.
+        grid.gmt.registration = header.registration
+        grid.gmt.gtype = header.gtype
+        return grid

pygmt/helpers/decorators.py

@@ -267,10 +267,12 @@
             - ``file`` will save the result to the file specified by the ``outfile``
               parameter.""",
     "outgrid": """
-        outgrid : str or None
-            Name of the output netCDF grid file. For writing a specific grid
-            file format or applying basic data operations to the output grid,
-            see :gmt-docs:`gmt.html#grd-inout-full` for the available modifiers.""",
+        outgrid
+            Name of the output netCDF grid file. If not specified, will return an
+            :class:`xarray.DataArray` object. For writing a specific grid file format or
+            applying basic data operations to the output grid, see
+            :gmt-docs:`gmt.html#grd-inout-full` for the available modifiers.
+        """,
     "panel": r"""
         panel : bool, int, or list
             [*row,col*\|\ *index*].

pygmt/src/binstats.py

@@ -3,21 +3,13 @@
 """
 
 from pygmt.clib import Session
-from pygmt.helpers import (
-    GMTTempFile,
-    build_arg_string,
-    fmt_docstring,
-    kwargs_to_strings,
-    use_alias,
-)
-from pygmt.io import load_dataarray
+from pygmt.helpers import build_arg_string, fmt_docstring, kwargs_to_strings, use_alias
 
 
 @fmt_docstring
 @use_alias(
     C="statistic",
     E="empty",
-    G="outgrid",
     I="spacing",
     N="normalize",
     R="region",
@@ -31,7 +23,7 @@
     r="registration",
 )
 @kwargs_to_strings(I="sequence", R="sequence", i="sequence_comma")
-def binstats(data, **kwargs):
+def binstats(data, outgrid: str | None = None, **kwargs):
     r"""
     Bin spatial data and determine statistics per bin.
 
@@ -110,13 +102,13 @@ def binstats(data, **kwargs):
         - None if ``outgrid`` is set (grid output will be stored in file set by
           ``outgrid``)
     """
-    with GMTTempFile(suffix=".nc") as tmpfile:
-        with Session() as lib:
-            with lib.virtualfile_in(check_kind="vector", data=data) as vintbl:
-                if (outgrid := kwargs.get("G")) is None:
-                    kwargs["G"] = outgrid = tmpfile.name  # output to tmpfile
-                lib.call_module(
-                    module="binstats", args=build_arg_string(kwargs, infile=vintbl)
-                )
-
-        return load_dataarray(outgrid) if outgrid == tmpfile.name else None
+    with Session() as lib:
+        with (
+            lib.virtualfile_in(check_kind="vector", data=data) as vintbl,
+            lib.virtualfile_out(kind="grid", fname=outgrid) as voutgrd,
+        ):
+            kwargs["G"] = voutgrd
+            lib.call_module(
+                module="binstats", args=build_arg_string(kwargs, infile=vintbl)
+            )
+            return lib.virtualfile_to_raster(vfname=voutgrd, outgrid=outgrid)