GMT_IMAGE: Implement the GMT_IMAGE.to_dataarray method for 3-band images

pygmt/clib/session.py

@@ -25,7 +25,7 @@
     vectors_to_arrays,
 )
 from pygmt.clib.loading import load_libgmt
-from pygmt.datatypes import _GMT_DATASET, _GMT_GRID
+from pygmt.datatypes import _GMT_DATASET, _GMT_GRID, _GMT_IMAGE
 from pygmt.exceptions import (
     GMTCLibError,
     GMTCLibNoSessionError,
@@ -1697,7 +1697,9 @@ def virtualfile_from_data(
 
     @contextlib.contextmanager
     def virtualfile_out(
-        self, kind: Literal["dataset", "grid"] = "dataset", fname: str | None = None
+        self,
+        kind: Literal["dataset", "grid", "image"] = "dataset",
+        fname: str | None = None,
     ):
         r"""
         Create a virtual file or an actual file for storing output data.
@@ -1710,8 +1712,8 @@ def virtualfile_out(
         Parameters
         ----------
         kind
-            The data kind of the virtual file to create. Valid values are ``"dataset"``
-            and ``"grid"``. Ignored if ``fname`` is specified.
+            The data kind of the virtual file to create. Valid values are ``"dataset"``,
+            ``"grid"``, and ``"image"``. Ignored if ``fname`` is specified.
         fname
             The name of the actual file to write the output data. No virtual file will
             be created.
@@ -1754,8 +1756,11 @@ def virtualfile_out(
             family, geometry = {
                 "dataset": ("GMT_IS_DATASET", "GMT_IS_PLP"),
                 "grid": ("GMT_IS_GRID", "GMT_IS_SURFACE"),
+                "image": ("GMT_IS_IMAGE", "GMT_IS_SURFACE"),
             }[kind]
-            with self.open_virtualfile(family, geometry, "GMT_OUT", None) as vfile:
+            with self.open_virtualfile(
+                family, geometry, "GMT_OUT|GMT_IS_REFERENCE", None
+            ) as vfile:
                 yield vfile
 
     def inquire_virtualfile(self, vfname: str) -> int:
@@ -1801,7 +1806,8 @@ def read_virtualfile(
             Name of the virtual file to read.
         kind
             Cast the data into a GMT data container. Valid values are ``"dataset"``,
-            ``"grid"`` and ``None``. If ``None``, will return a ctypes void pointer.
+            ``"grid"``, ``"image"`` and ``None``. If ``None``, will return a ctypes void
+            pointer.
 
         Examples
         --------
@@ -1849,9 +1855,9 @@ def read_virtualfile(
         # _GMT_DATASET).
         if kind is None:  # Return the ctypes void pointer
             return pointer
-        if kind in ["image", "cube"]:
+        if kind == "cube":
             raise NotImplementedError(f"kind={kind} is not supported yet.")
-        dtype = {"dataset": _GMT_DATASET, "grid": _GMT_GRID}[kind]
+        dtype = {"dataset": _GMT_DATASET, "grid": _GMT_GRID, "image": _GMT_IMAGE}[kind]
         return ctp.cast(pointer, ctp.POINTER(dtype))
 
     def virtualfile_to_dataset(

pygmt/datatypes/__init__.py

@@ -4,3 +4,4 @@
 
 from pygmt.datatypes.dataset import _GMT_DATASET
 from pygmt.datatypes.grid import _GMT_GRID
+from pygmt.datatypes.image import _GMT_IMAGE

pygmt/datatypes/image.py

@@ -0,0 +1,108 @@
+"""
+Wrapper for the GMT_IMAGE data type.
+"""
+
+import ctypes as ctp
+from typing import ClassVar
+
+import numpy as np
+import xarray as xr
+from pygmt.datatypes.header import _GMT_GRID_HEADER
+
+
+class _GMT_IMAGE(ctp.Structure):  # noqa: N801
+    """
+    GMT image data structure.
+
+    Examples
+    --------
+    >>> from pygmt.clib import Session
+    >>> import numpy as np
+    >>> import xarray as xr
+    >>> import rioxarray
+
+    >>> with Session() as lib:
+    ...     with lib.virtualfile_out(kind="image") as voutimg:
+    ...         lib.call_module("read", f"@earth_day_01d {voutimg} -Ti")
+    ...         ds = lib.read_virtualfile(vfname=voutimg, kind="image").contents
+    ...         header = ds.header.contents
+    ...         pad = header.pad[:]
+    ...         print(ds.type, header.n_bands, header.n_rows, header.n_columns)
+    ...         print(header.pad[:])
+    ...         data = np.reshape(
+    ...             ds.data[: header.n_bands * header.mx * header.my],
+    ...             (header.my, header.mx, header.n_bands),
+    ...         )
+    ...         data = data[pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1], :]
+    ...         x = ds.x[: header.n_columns]
+    ...         y = ds.y[: header.n_rows]
+    >>> da = xr.DataArray(
+    ...     data=data,
+    ...     dims=["y", "x", "band"],
+    ...     coords={"y": y, "x": x, "band": [1, 2, 3]},
+    ... )
+    >>> da = da.transpose("band", "y", "x")
+    >>> da = da.sortby(list(data.dims))
+    >>> da.plot.imshow()
+    """
+
+    _fields_: ClassVar = [
+        # Data type, e.g. GMT_FLOAT
+        ("type", ctp.c_int),
+        # Array with color lookup values
+        ("colormap", ctp.POINTER(ctp.c_int)),
+        # Number of colors in a paletted image
+        ("n_indexed_colors", ctp.c_int),
+        # Pointer to full GMT header for the image
+        ("header", ctp.POINTER(_GMT_GRID_HEADER)),
+        # Pointer to actual image
+        ("data", ctp.POINTER(ctp.c_ubyte)),
+        # Pointer to an optional transparency layer stored in a separate variable
+        ("alpha", ctp.POINTER(ctp.c_ubyte)),
+        # Color interpolation
+        ("color_interp", ctp.c_char_p),
+        # Pointer to the x-coordinate vector
+        ("x", ctp.POINTER(ctp.c_double)),
+        # Pointer to the y-coordinate vector
+        ("y", ctp.POINTER(ctp.c_double)),
+        # Book-keeping variables "hidden" from the API
+        ("hidden", ctp.c_void_p),
+    ]
+
+    def to_dataarray(self) -> xr.DataArray:
+        """
+        Convert a _GMT_IMAGE object to an :class:`xarray.DataArray` object.
+
+        Returns
+        -------
+        dataarray
+            A :class:`xarray.DataArray` object.
+        """
+
+        # Get image header
+        header: _GMT_GRID_HEADER = self.header.contents
+
+        # Get DataArray without padding
+        pad = header.pad[:]
+        data: np.ndarray = np.reshape(
+            a=self.data[: header.n_bands * header.mx * header.my],
+            newshape=(header.my, header.mx, header.n_bands),
+        )[pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1], :]
+
+        # Get x and y coordinates
+        coords: dict[str, list | np.ndarray] = {
+            "x": self.x[: header.n_columns],
+            "y": self.y[: header.n_rows],
+            "band": np.array([0, 1, 2], dtype=np.uint8),
+        }
+
+        # Create the xarray.DataArray object
+        image = xr.DataArray(
+            data=data,
+            coords=coords,
+            dims=("y", "x", "band"),
+            name=header.name,
+            attrs=header.data_attrs,
+        )
+
+        return image