From b59efe01e2885d02f5c9abc0be08713821c9b265 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sat, 23 Mar 2024 10:53:25 +0800 Subject: [PATCH 01/22] Wrap GMT's standard data type GMT_CUBE for cubes --- pygmt/clib/session.py | 15 +++--- pygmt/datatypes/__init__.py | 1 + pygmt/datatypes/cube.py | 93 +++++++++++++++++++++++++++++++++++++ 3 files changed, 102 insertions(+), 7 deletions(-) create mode 100644 pygmt/datatypes/cube.py diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 1b8b5483a28..de162d124c0 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -26,7 +26,7 @@ vectors_to_arrays, ) from pygmt.clib.loading import load_libgmt -from pygmt.datatypes import _GMT_DATASET, _GMT_GRID +from pygmt.datatypes import _GMT_CUBE, _GMT_DATASET, _GMT_GRID from pygmt.exceptions import ( GMTCLibError, GMTCLibNoSessionError, @@ -1789,7 +1789,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", "cube"] = "dataset", + fname: str | None = None, ) -> Generator[str, None, None]: r""" Create a virtual file or an actual file for storing output data. @@ -1846,6 +1848,7 @@ def virtualfile_out( family, geometry = { "dataset": ("GMT_IS_DATASET", "GMT_IS_PLP"), "grid": ("GMT_IS_GRID", "GMT_IS_SURFACE"), + "cube": ("GMT_IS_CUBE", "GMT_IS_VOLUME"), }[kind] with self.open_virtualfile(family, geometry, "GMT_OUT", None) as vfile: yield vfile @@ -1880,9 +1883,7 @@ def inquire_virtualfile(self, vfname: str) -> int: return c_inquire_virtualfile(self.session_pointer, vfname.encode()) def read_virtualfile( - self, - vfname: str, - kind: Literal["dataset", "grid", "image", "cube", None] = None, + self, vfname: str, kind: Literal["dataset", "grid", "cube", None] = None ): """ Read data from a virtual file and optionally cast into a GMT data container. @@ -1943,9 +1944,9 @@ def read_virtualfile( # _GMT_DATASET). if kind is None: # Return the ctypes void pointer return pointer - if kind in {"image", "cube"}: + if kind == "image": raise NotImplementedError(f"kind={kind} is not supported yet.") - dtype = {"dataset": _GMT_DATASET, "grid": _GMT_GRID}[kind] + dtype = {"dataset": _GMT_DATASET, "grid": _GMT_GRID, "cube": _GMT_CUBE}[kind] return ctp.cast(pointer, ctp.POINTER(dtype)) def virtualfile_to_dataset( diff --git a/pygmt/datatypes/__init__.py b/pygmt/datatypes/__init__.py index 237a050a9f7..16627ed5798 100644 --- a/pygmt/datatypes/__init__.py +++ b/pygmt/datatypes/__init__.py @@ -2,5 +2,6 @@ Wrappers for GMT data types. """ +from pygmt.datatypes.cube import _GMT_CUBE from pygmt.datatypes.dataset import _GMT_DATASET from pygmt.datatypes.grid import _GMT_GRID diff --git a/pygmt/datatypes/cube.py b/pygmt/datatypes/cube.py new file mode 100644 index 00000000000..b51341c36ab --- /dev/null +++ b/pygmt/datatypes/cube.py @@ -0,0 +1,93 @@ +""" +Wrapper for the GMT_CUBE 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, + GMT_GRID_UNIT_LEN80, + GMT_GRID_VARNAME_LEN80, + _parse_nameunits, + gmt_grdfloat, +) + + +class _GMT_CUBE(ctp.Structure): # noqa: N801 + """ + GMT cube data structure for 3D data. + """ + + _fields_: ClassVar = [ + # Pointer to full GMT 2-D header for a layer (common to all layers) + ("header", ctp.POINTER(_GMT_GRID_HEADER)), + # Pointer to the gmt_grdfloat 3-D cube - a stack of 2-D padded grids + ("data", ctp.POINTER(gmt_grdfloat)), + # Vector of x coordinates common to all layers + ("x", ctp.POINTER(ctp.c_double)), + # Vector of y coordinates common to all layers + ("y", ctp.POINTER(ctp.c_double)), + # Low-level information for GMT use only + ("hidden", ctp.c_void_p), + # GMT_CUBE_IS_STACK if input dataset was a list of 2-D grids rather than a + # single cube + ("mode", ctp.c_uint), + # Minimum/max z values (complements header->wesn[4]) + ("z_range", ctp.c_double * 2), + # z increment (complements inc[2]) (0 if variable z spacing) + ("z_inc", ctp.c_double), + # Array of z values (complements x, y) + ("z", ctp.POINTER(ctp.c_double)), + # Name of the 3-D variable, if read from file (or empty if just one) + ("name", ctp.c_char * GMT_GRID_VARNAME_LEN80), + # Units in 3rd direction (complements x_units, y_units, z_units) + ("units", ctp.c_char * GMT_GRID_UNIT_LEN80), + ] + + def to_dataarray(self): + """ + Convert the GMT_CUBE to an xarray.DataArray. + + Returns + ------- + xarray.DataArray: The data array representation of the GMT_CUBE. + """ + # The grid header + header = self.header.contents + + name = "cube" + # Dimensions and attributes + dims = header.dims + dim_attrs = header.dim_attrs + + # Patch for the 3rd dimension + dims.append("z") + z_attrs = {"actual_range": np.array(self.z_range[:]), "axis": "Z"} + long_name, units = _parse_nameunits(self.units.decode()) + if long_name: + z_attrs["long_name"] = long_name + if units: + z_attrs["units"] = units + dim_attrs.append(z_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]), + # header->n_bands is used for the number of layers for 3-D cubes + (dims[2], self.z[: header.n_bands], dim_attrs[1]), + ] + + # The data array without paddings + pad = header.pad[:] + data = np.reshape( + self.data[: header.mx * header.my * header.n_bands], + (header.my, header.mx, header.n_bands), + )[pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1], :] + + # Create the xarray.DataArray object + grid = xr.DataArray(data, coords=coords, name=name, attrs=header.dataA_attrs) + return grid From 277a30d4566df1859a766e8c855e0dcc8bd41891 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 19 Jul 2024 16:52:35 +0800 Subject: [PATCH 02/22] Remove the to_dataarray method --- pygmt/datatypes/cube.py | 48 ----------------------------------------- 1 file changed, 48 deletions(-) diff --git a/pygmt/datatypes/cube.py b/pygmt/datatypes/cube.py index b51341c36ab..15ca63adec5 100644 --- a/pygmt/datatypes/cube.py +++ b/pygmt/datatypes/cube.py @@ -5,13 +5,10 @@ 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_GRID_UNIT_LEN80, GMT_GRID_VARNAME_LEN80, - _parse_nameunits, gmt_grdfloat, ) @@ -46,48 +43,3 @@ class _GMT_CUBE(ctp.Structure): # noqa: N801 # Units in 3rd direction (complements x_units, y_units, z_units) ("units", ctp.c_char * GMT_GRID_UNIT_LEN80), ] - - def to_dataarray(self): - """ - Convert the GMT_CUBE to an xarray.DataArray. - - Returns - ------- - xarray.DataArray: The data array representation of the GMT_CUBE. - """ - # The grid header - header = self.header.contents - - name = "cube" - # Dimensions and attributes - dims = header.dims - dim_attrs = header.dim_attrs - - # Patch for the 3rd dimension - dims.append("z") - z_attrs = {"actual_range": np.array(self.z_range[:]), "axis": "Z"} - long_name, units = _parse_nameunits(self.units.decode()) - if long_name: - z_attrs["long_name"] = long_name - if units: - z_attrs["units"] = units - dim_attrs.append(z_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]), - # header->n_bands is used for the number of layers for 3-D cubes - (dims[2], self.z[: header.n_bands], dim_attrs[1]), - ] - - # The data array without paddings - pad = header.pad[:] - data = np.reshape( - self.data[: header.mx * header.my * header.n_bands], - (header.my, header.mx, header.n_bands), - )[pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1], :] - - # Create the xarray.DataArray object - grid = xr.DataArray(data, coords=coords, name=name, attrs=header.dataA_attrs) - return grid From 757321a263c89b39622732f0200823014f16aa5d Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 19 Jul 2024 16:58:41 +0800 Subject: [PATCH 03/22] Finalize the GMT_CUBE wrapper --- pygmt/datatypes/cube.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/pygmt/datatypes/cube.py b/pygmt/datatypes/cube.py index 15ca63adec5..8bac0d14659 100644 --- a/pygmt/datatypes/cube.py +++ b/pygmt/datatypes/cube.py @@ -15,7 +15,12 @@ class _GMT_CUBE(ctp.Structure): # noqa: N801 """ - GMT cube data structure for 3D data. + GMT cube data structure for 3-D data. + + The GMT_CUBE structure is a extension of the GMT_GRID structure to handle 3-D data + cubes. It requires a 2-D grid header and extended parameters for the 3rd dimension. + + header->n_bands is used for the number of layers in 3-D cubes. """ _fields_: ClassVar = [ @@ -29,10 +34,10 @@ class _GMT_CUBE(ctp.Structure): # noqa: N801 ("y", ctp.POINTER(ctp.c_double)), # Low-level information for GMT use only ("hidden", ctp.c_void_p), - # GMT_CUBE_IS_STACK if input dataset was a list of 2-D grids rather than a - # single cube + # mode=GMT_CUBE_IS_STACK means the input dataset was a list of 2-D grids, rather + # than a single cube. ("mode", ctp.c_uint), - # Minimum/max z values (complements header->wesn[4]) + # Minimum/maximum z values (complements header->wesn[4]) ("z_range", ctp.c_double * 2), # z increment (complements inc[2]) (0 if variable z spacing) ("z_inc", ctp.c_double), From f77d4128f8d8e4e13d8926c781a03abffea06567 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 19 Jul 2024 17:06:10 +0800 Subject: [PATCH 04/22] Update the docstrings --- pygmt/clib/session.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index de162d124c0..ee55359c016 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -1804,8 +1804,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 ``"cube"``. Ignored if ``fname`` is specified. fname The name of the actual file to write the output data. No virtual file will be created. @@ -1894,7 +1894,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"``, ``"cube"`` and ``None``. If ``None``, will return a ctypes void + pointer. Returns ------- From ecf57116d9d2b680fd584057feff85d6bcb9660f Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 19 Jul 2024 17:10:59 +0800 Subject: [PATCH 05/22] Add cube support in Session.read_data --- pygmt/clib/session.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index ee55359c016..6ed97ddf1fd 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -1070,7 +1070,7 @@ def put_matrix(self, dataset, matrix, pad=0): def read_data( self, infile: str, - kind: Literal["dataset", "grid"], + kind: Literal["dataset", "grid", "cube"], family: str | None = None, geometry: str | None = None, mode: str = "GMT_READ_NORMAL", @@ -1088,8 +1088,8 @@ def read_data( infile The input file name. kind - The data kind of the input file. Valid values are ``"dataset"`` and - ``"grid"``. + The data kind of the input file. Valid values are ``"dataset"``, ``"grid"`` + and ``"cube"``. family A valid GMT data family name (e.g., ``"GMT_IS_DATASET"``). See the ``FAMILIES`` attribute for valid names. If ``None``, will determine the data @@ -1140,6 +1140,7 @@ def read_data( _family, _geometry, dtype = { "dataset": ("GMT_IS_DATASET", "GMT_IS_PLP", _GMT_DATASET), "grid": ("GMT_IS_GRID", "GMT_IS_SURFACE", _GMT_GRID), + "cube": ("GMT_IS_CUBE", "GMT_IS_VOLUME", _GMT_CUBE), }[kind] if family is None: family = _family From 810fb12a248333ecad75bf34b1ffaeab967447e8 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 19 Jul 2024 18:16:03 +0800 Subject: [PATCH 06/22] Add two tests for reading grid/image as GMT_CUBE --- pygmt/tests/test_clib_read_data.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/pygmt/tests/test_clib_read_data.py b/pygmt/tests/test_clib_read_data.py index 43978b291c2..57660c74b53 100644 --- a/pygmt/tests/test_clib_read_data.py +++ b/pygmt/tests/test_clib_read_data.py @@ -132,6 +132,24 @@ def test_clib_read_data_grid_actual_image(): ) +def test_clib_read_data_cube_actual_grid(): + """ + Test the Session.read_data method for cube, but actually the file is a grid. + """ + with Session() as lib: + with pytest.raises(GMTCLibError): + lib.read_data("@earth_relief_01d_p", kind="cube", mode="GMT_CONTAINER_ONLY") + + +def test_clib_read_data_cube_actual_image(): + """ + Test the Session.read_data method for cube, but actually the file is an image. + """ + with Session() as lib: + with pytest.raises(GMTCLibError): + lib.read_data("@earth_day_01d_p", kind="cube", mode="GMT_CONTAINER_ONLY") + + def test_clib_read_data_fails(): """ Test that the Session.read_data method raises an exception if there are errors. From 2823168a6f1bc1b03ef7cd38587e0e552f770c95 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sat, 27 Jul 2024 16:42:55 +0800 Subject: [PATCH 07/22] Fix some typos --- pygmt/clib/session.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 757a7552703..77ef28e7f75 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -1090,7 +1090,7 @@ def read_data( The input file name. kind The data kind of the input file. Valid values are ``"dataset"``, ``"grid"``, - ``"image"`` and ``"cube"``, + ``"image"`` and ``"cube"``. family A valid GMT data family name (e.g., ``"GMT_IS_DATASET"``). See the ``FAMILIES`` attribute for valid names. If ``None``, will determine the data @@ -1895,7 +1895,9 @@ def inquire_virtualfile(self, vfname: str) -> int: return c_inquire_virtualfile(self.session_pointer, vfname.encode()) def read_virtualfile( - self, vfname: str, kind: Literal["dataset", "grid", "cube", 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. From 6b1bc4091256a237e56f2918c958add71a9b3e09 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sat, 27 Jul 2024 16:44:54 +0800 Subject: [PATCH 08/22] Sort data kind by importance --- pygmt/clib/session.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 77ef28e7f75..cd903d24efc 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -1071,7 +1071,7 @@ def put_matrix(self, dataset, matrix, pad=0): def read_data( self, infile: str, - kind: Literal["dataset", "grid", "cube", "image"], + kind: Literal["dataset", "grid", "image", "cube"], family: str | None = None, geometry: str | None = None, mode: str = "GMT_READ_NORMAL", @@ -1141,8 +1141,8 @@ def read_data( _family, _geometry, dtype = { "dataset": ("GMT_IS_DATASET", "GMT_IS_PLP", _GMT_DATASET), "grid": ("GMT_IS_GRID", "GMT_IS_SURFACE", _GMT_GRID), - "cube": ("GMT_IS_CUBE", "GMT_IS_VOLUME", _GMT_CUBE), "image": ("GMT_IS_IMAGE", "GMT_IS_SURFACE", _GMT_IMAGE), + "cube": ("GMT_IS_CUBE", "GMT_IS_VOLUME", _GMT_CUBE), }[kind] if family is None: family = _family @@ -1800,7 +1800,7 @@ def virtualfile_from_data( @contextlib.contextmanager def virtualfile_out( self, - kind: Literal["dataset", "grid", "cube", "image"] = "dataset", + kind: Literal["dataset", "grid", "image", "cube"] = "dataset", fname: str | None = None, ) -> Generator[str, None, None]: r""" @@ -1858,8 +1858,8 @@ def virtualfile_out( family, geometry = { "dataset": ("GMT_IS_DATASET", "GMT_IS_PLP"), "grid": ("GMT_IS_GRID", "GMT_IS_SURFACE"), - "cube": ("GMT_IS_CUBE", "GMT_IS_VOLUME"), "image": ("GMT_IS_IMAGE", "GMT_IS_SURFACE"), + "cube": ("GMT_IS_CUBE", "GMT_IS_VOLUME"), }[kind] direction = "GMT_OUT|GMT_IS_REFERENCE" if kind == "image" else "GMT_OUT" with self.open_virtualfile(family, geometry, direction, None) as vfile: @@ -1962,8 +1962,8 @@ def read_virtualfile( dtype = { "dataset": _GMT_DATASET, "grid": _GMT_GRID, - "cube": _GMT_CUBE, "image": _GMT_IMAGE, + "cube": _GMT_CUBE, }[kind] return ctp.cast(pointer, ctp.POINTER(dtype)) From fec544cb45fce2c5963b7e4fa6e8e5ddb96bc0e9 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 9 Aug 2024 14:11:01 +0800 Subject: [PATCH 09/22] Cache cube.nc --- pygmt/helpers/caching.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pygmt/helpers/caching.py b/pygmt/helpers/caching.py index 714f12d890e..572eca691c1 100644 --- a/pygmt/helpers/caching.py +++ b/pygmt/helpers/caching.py @@ -74,6 +74,7 @@ def cache_data(): "@Table_5_11_mean.xyz", "@capitals.gmt", "@circuit.png", + "@cube.nc", "@earth_relief_20m_holes.grd", "@fractures_06.txt", "@hotspots.txt", From 2647e41239d3525216fa2227704003ece2197b94 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Thu, 19 Sep 2024 22:17:29 +0800 Subject: [PATCH 10/22] Add a doctest --- pygmt/datatypes/cube.py | 56 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/pygmt/datatypes/cube.py b/pygmt/datatypes/cube.py index 8bac0d14659..e56f372f02f 100644 --- a/pygmt/datatypes/cube.py +++ b/pygmt/datatypes/cube.py @@ -21,6 +21,62 @@ class _GMT_CUBE(ctp.Structure): # noqa: N801 cubes. It requires a 2-D grid header and extended parameters for the 3rd dimension. header->n_bands is used for the number of layers in 3-D cubes. + + Examples + -------- + >>> import numpy as np + >>> from pygmt import which + >>> from pygmt.clib import Session + >>> cubefile = which("@cube.nc", download="c") + >>> with Session() as lib: + ... with lib.virtualfile_out(kind="cube") as voutcube: + ... lib.call_module("read", [cubefile, voutcube, "-Tu", "-Vd"]) + ... # Read the cube from the virtual file + ... cube = lib.read_virtualfile(vfname=voutcube, kind="cube").contents + ... # The cube header + ... header = cube.header.contents + ... # Access the header properties + ... print(header.n_rows, header.n_columns, header.registration) + ... print(header.wesn[:], header.inc[:]) + ... print(header.z_scale_factor, header.z_add_offset) + ... print(header.x_units, header.y_units, header.z_units) + ... 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.xy_off) + ... # Cube-specific attributes. + ... print(cube.mode, cube.z_range[:], cube.z_inc, cube.name, cube.units) + ... # The x, y, and z coordinates + ... x = cube.x[: header.n_columns] + ... y = cube.y[: header.n_rows] + ... z = cube.z[: header.n_bands] + ... # The data array (with paddings) + ... data = np.reshape( + ... cube.data[: header.n_bands * header.mx * header.my], + ... (header.my, header.mx, header.n_bands), + ... ) + ... # The data array (without paddings) + ... pad = header.pad[:] + ... data = data[pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1], :] + 11 11 0 + [0.0, 10.0, 0.0, 10.0] [1.0, 1.0] + 1.0 0.0 + b'x' b'y' b'cube' + 121 226 0 + 18 4 15 15 + [2, 2, 2, 2] + b'' 0.0 + 0 [1.0, 5.0] 0.0 b'' b'z' + >>> x + [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0] + >>> y + [10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0, 0.0] + >>> z + [1.0, 2.0, 3.0, 5.0] + >>> data.shape + (11, 11, 4) + >>> #data.min(), data.max() # The min/max are wrong. Upstream bug? + >>> #(-29.399999618530273, 169.39999389648438) """ _fields_: ClassVar = [ From 0b64d70e23ab982a571a2674f81f570625060725 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Thu, 19 Sep 2024 22:20:43 +0800 Subject: [PATCH 11/22] Fix styles --- pygmt/datatypes/cube.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygmt/datatypes/cube.py b/pygmt/datatypes/cube.py index e56f372f02f..09c48d0b62f 100644 --- a/pygmt/datatypes/cube.py +++ b/pygmt/datatypes/cube.py @@ -75,8 +75,8 @@ class _GMT_CUBE(ctp.Structure): # noqa: N801 [1.0, 2.0, 3.0, 5.0] >>> data.shape (11, 11, 4) - >>> #data.min(), data.max() # The min/max are wrong. Upstream bug? - >>> #(-29.399999618530273, 169.39999389648438) + >>> # data.min(), data.max() # The min/max are wrong. Upstream bug? + >>> # (-29.399999618530273, 169.39999389648438) """ _fields_: ClassVar = [ From 80999c614f19a9dd184d8582879102fa5e7db3fc Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 30 Sep 2024 21:38:00 +0800 Subject: [PATCH 12/22] Updates --- pygmt/datatypes/cube.py | 21 ++++++++++----------- pygmt/tests/test_clib_read_data.py | 2 +- 2 files changed, 11 insertions(+), 12 deletions(-) diff --git a/pygmt/datatypes/cube.py b/pygmt/datatypes/cube.py index 09c48d0b62f..363febd1fc0 100644 --- a/pygmt/datatypes/cube.py +++ b/pygmt/datatypes/cube.py @@ -47,14 +47,13 @@ class _GMT_CUBE(ctp.Structure): # noqa: N801 ... # Cube-specific attributes. ... print(cube.mode, cube.z_range[:], cube.z_inc, cube.name, cube.units) ... # The x, y, and z coordinates - ... x = cube.x[: header.n_columns] - ... y = cube.y[: header.n_rows] - ... z = cube.z[: header.n_bands] + ... x = np.ctypeslib.as_array(cube.x, shape=(header.n_columns,)).copy() + ... y = np.ctypeslib.as_array(cube.y, shape=(header.n_rows,)).copy() + ... z = np.ctypeslib.as_array(cube.z, shape=(header.n_bands,)).copy() ... # The data array (with paddings) - ... data = np.reshape( - ... cube.data[: header.n_bands * header.mx * header.my], - ... (header.my, header.mx, header.n_bands), - ... ) + ... data = np.ctypeslib.as_array( + ... cube.data, shape=(header.my, header.mx, header.n_bands) + ... ).copy() ... # The data array (without paddings) ... pad = header.pad[:] ... data = data[pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1], :] @@ -68,15 +67,15 @@ class _GMT_CUBE(ctp.Structure): # noqa: N801 b'' 0.0 0 [1.0, 5.0] 0.0 b'' b'z' >>> x - [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0] + array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]) >>> y - [10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0, 0.0] + array([10., 9., 8., 7., 6., 5., 4., 3., 2., 1., 0.]) >>> z - [1.0, 2.0, 3.0, 5.0] + array([1., 2., 3., 5.]) >>> data.shape (11, 11, 4) >>> # data.min(), data.max() # The min/max are wrong. Upstream bug? - >>> # (-29.399999618530273, 169.39999389648438) + >>> # (-29.4, 169.4) """ _fields_: ClassVar = [ diff --git a/pygmt/tests/test_clib_read_data.py b/pygmt/tests/test_clib_read_data.py index 00f8fba0f5b..4da7fa61b30 100644 --- a/pygmt/tests/test_clib_read_data.py +++ b/pygmt/tests/test_clib_read_data.py @@ -216,7 +216,7 @@ def test_clib_read_data_cube_actual_image(): """ with Session() as lib: with pytest.raises(GMTCLibError): - lib.read_data("@earth_day_01d_p", kind="cube", mode="GMT_CONTAINER_ONLY") + lib.read_data("@earth_day_01d", kind="cube", mode="GMT_CONTAINER_ONLY") def test_clib_read_data_fails(): From 723d9f2dbc7617e38fa5bc1ca255e928fb40ddb2 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 23 Aug 2026 12:30:59 +0800 Subject: [PATCH 13/22] Fix styling issue --- pygmt/datatypes/cube.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygmt/datatypes/cube.py b/pygmt/datatypes/cube.py index 363febd1fc0..af584a1d488 100644 --- a/pygmt/datatypes/cube.py +++ b/pygmt/datatypes/cube.py @@ -13,7 +13,7 @@ ) -class _GMT_CUBE(ctp.Structure): # noqa: N801 +class _GMT_CUBE(ctp.Structure): # ruff: ignore[invalid-class-name] """ GMT cube data structure for 3-D data. From 0117af30ee17af529ae5787c0715f3849594326e Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 23 Aug 2026 16:42:22 +0800 Subject: [PATCH 14/22] Fix for cubes --- pygmt/datatypes/cube.py | 153 ++++++++++++++++++++++++-- pygmt/datatypes/header.py | 9 +- pygmt/tests/test_clib_read_data.py | 87 +++++++++++++-- pygmt/tests/test_clib_virtualfiles.py | 21 +++- 4 files changed, 250 insertions(+), 20 deletions(-) diff --git a/pygmt/datatypes/cube.py b/pygmt/datatypes/cube.py index af584a1d488..e0d14193f97 100644 --- a/pygmt/datatypes/cube.py +++ b/pygmt/datatypes/cube.py @@ -5,10 +5,13 @@ 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_GRID_UNIT_LEN80, GMT_GRID_VARNAME_LEN80, + _parse_nameunits, gmt_grdfloat, ) @@ -20,17 +23,22 @@ class _GMT_CUBE(ctp.Structure): # ruff: ignore[invalid-class-name] The GMT_CUBE structure is a extension of the GMT_GRID structure to handle 3-D data cubes. It requires a 2-D grid header and extended parameters for the 3rd dimension. - header->n_bands is used for the number of layers in 3-D cubes. + ``header->n_bands`` is used for the number of layers in 3-D cubes. + + The ``data`` array is a stack of 2-D padded layers, i.e., layer ``k`` starts at + offset ``k * header.size``. Note that ``header.size`` is the allocated length of one + padded layer, which can exceed ``header.my * header.mx``. Examples -------- >>> import numpy as np >>> from pygmt import which >>> from pygmt.clib import Session + >>> from pygmt.datatypes import _GMT_CUBE >>> cubefile = which("@cube.nc", download="c") >>> with Session() as lib: ... with lib.virtualfile_out(kind="cube") as voutcube: - ... lib.call_module("read", [cubefile, voutcube, "-Tu", "-Vd"]) + ... lib.call_module("read", [cubefile, voutcube, "-Tu"]) ... # Read the cube from the virtual file ... cube = lib.read_virtualfile(vfname=voutcube, kind="cube").contents ... # The cube header @@ -50,13 +58,16 @@ class _GMT_CUBE(ctp.Structure): # ruff: ignore[invalid-class-name] ... x = np.ctypeslib.as_array(cube.x, shape=(header.n_columns,)).copy() ... y = np.ctypeslib.as_array(cube.y, shape=(header.n_rows,)).copy() ... z = np.ctypeslib.as_array(cube.z, shape=(header.n_bands,)).copy() - ... # The data array (with paddings) + ... # The data array (one padded layer per row) ... data = np.ctypeslib.as_array( - ... cube.data, shape=(header.my, header.mx, header.n_bands) + ... cube.data, shape=(header.n_bands, header.size) ... ).copy() - ... # The data array (without paddings) + ... # Reshape the layers to 2-D and strip the paddings ... pad = header.pad[:] - ... data = data[pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1], :] + ... data = data[:, : header.my * header.mx].reshape( + ... header.n_bands, header.my, header.mx + ... ) + ... data = data[:, pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1]] 11 11 0 [0.0, 10.0, 0.0, 10.0] [1.0, 1.0] 1.0 0.0 @@ -73,9 +84,31 @@ class _GMT_CUBE(ctp.Structure): # ruff: ignore[invalid-class-name] >>> z array([1., 2., 3., 5.]) >>> data.shape - (11, 11, 4) - >>> # data.min(), data.max() # The min/max are wrong. Upstream bug? - >>> # (-29.4, 169.4) + (4, 11, 11) + >>> print(data.min(), data.max()) + 0.0 140.0 + >>> # GMT stores rows north-first, so row 0 is y=10 and the last row is y=0. + >>> print(data[0, :, :]) + [[ 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. 100.] + [ 0. 9. 18. 27. 36. 45. 54. 63. 72. 81. 90.] + [ 0. 8. 16. 24. 32. 40. 48. 56. 64. 72. 80.] + [ 0. 7. 14. 21. 28. 35. 42. 49. 56. 63. 70.] + [ 0. 6. 12. 18. 24. 30. 36. 42. 48. 54. 60.] + [ 0. 5. 10. 15. 20. 25. 30. 35. 40. 45. 50.] + [ 0. 4. 8. 12. 16. 20. 24. 28. 32. 36. 40.] + [ 0. 3. 6. 9. 12. 15. 18. 21. 24. 27. 30.] + [ 0. 2. 4. 6. 8. 10. 12. 14. 16. 18. 20.] + [ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.] + [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]] + >>> # The northernmost row of every layer. The four layers are the same X*Y grid + >>> # scaled by 1, 1.1, 1.2 and 1.4, so this also pins the layer order. + >>> print(data[:, 0, :]) + [[ 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. 100.] + [ 0. 11. 22. 33. 44. 55. 66. 77. 88. 99. 110.] + [ 0. 12. 24. 36. 48. 60. 72. 84. 96. 108. 120.] + [ 0. 14. 28. 42. 56. 70. 84. 98. 112. 126. 140.]] + >>> print(data.max(axis=(1, 2))) + [100. 110. 120. 140.] """ _fields_: ClassVar = [ @@ -103,3 +136,105 @@ class _GMT_CUBE(ctp.Structure): # ruff: ignore[invalid-class-name] # Units in 3rd direction (complements x_units, y_units, z_units) ("units", ctp.c_char * GMT_GRID_UNIT_LEN80), ] + + def _parse_dimension(self) -> tuple[str, dict]: + """ + Get the name and attributes of the 3rd dimension. + + Unlike the x/y dimensions, the 3rd dimension is described by the cube itself + rather than by the 2-D grid header. ``self.units`` holds the dimension's + "long_name [units]" string, and ``self.z_range`` its actual range. + """ + attrs: dict = {} + long_name, units = _parse_nameunits(self.units.decode()) + if long_name: + attrs["long_name"] = long_name + if units: + attrs["units"] = units + attrs["axis"] = "Z" + attrs["actual_range"] = np.array(self.z_range[:]) + return "z", attrs + + def to_xarray(self) -> xr.DataArray: + """ + Convert a _GMT_CUBE object to a :class:`xarray.DataArray` object. + + Returns + ------- + dataarray + A 3-D :class:`xr.DataArray` object with dimensions ordered as + (3rd dimension, y, x). + + Examples + -------- + >>> from pygmt import which + >>> from pygmt.clib import Session + >>> cubefile = which("@cube.nc", download="c") + >>> with Session() as lib: + ... with lib.virtualfile_out(kind="cube") as voutcube: + ... lib.call_module("read", [cubefile, voutcube, "-Tu"]) + ... # Read the cube from the virtual file + ... cube = lib.read_virtualfile(voutcube, kind="cube") + ... # Convert to xarray.DataArray and use it later + ... da = cube.contents.to_xarray() + >>> da.name, da.dims, da.shape + ('z', ('z', 'y', 'x'), (4, 11, 11)) + >>> da.coords["z"] + Size: 32B + array([1., 2., 3., 5.]) + Coordinates: + * z (z) float64 32B 1.0 2.0 3.0 5.0 + Attributes: + long_name: z + axis: Z + actual_range: [1. 5.] + >>> # The four layers are the same X*Y grid scaled by 1, 1.1, 1.2 and 1.4. + >>> da.max(dim=("y", "x")).values + array([100., 110., 120., 140.], dtype=float32) + >>> da.gmt.registration, da.gmt.gtype + (, ) + """ + header = self.header.contents + + # The y/x dimensions come from the 2-D grid header; the 3rd one from the cube. + dims, dim_attrs = header.dims, header.dim_attrs + zdim, zdim_attrs = self._parse_dimension() + + # The coordinates, given as a tuple of the form (dims, data, attrs) + x = np.ctypeslib.as_array(self.x, shape=(header.n_columns,)).copy() + y = np.ctypeslib.as_array(self.y, shape=(header.n_rows,)).copy() + z = np.ctypeslib.as_array(self.z, shape=(header.n_bands,)).copy() + coords = [ + (zdim, z, zdim_attrs), + (dims[0], y, dim_attrs[0]), + (dims[1], x, dim_attrs[1]), + ] + + # The data array. The cube is a stack of 2-D padded layers, i.e., layer k + # starts at offset k * header.size, which can exceed header.my * header.mx. + data = np.ctypeslib.as_array( + self.data, shape=(header.n_bands, header.size) + ).copy() + pad = header.pad[:] + data = data[:, : header.my * header.mx].reshape( + header.n_bands, header.my, header.mx + ) + data = data[:, pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1]] + + # Create the xarray.DataArray object + cube = xr.DataArray( + data, coords=coords, name=header.name, attrs=header.data_attrs + ) + + # Flip the coordinates and data if necessary so that coordinates are ascending. + # `cube.sortby(list(cube.dims))` sometimes causes crashes. + # The solution comes from https://github.com/pydata/xarray/discussions/6695. + for dim in cube.dims: + if cube[dim].size > 1 and cube[dim][0] > cube[dim][1]: + cube = cube.isel({dim: slice(None, None, -1)}) + + # Set GMT accessors. + # Must put at the end, otherwise info gets lost after certain grid operations. + cube.gmt.registration = header.registration + cube.gmt.gtype = header.gtype + return cube diff --git a/pygmt/datatypes/header.py b/pygmt/datatypes/header.py index afd60c08c28..619c3e5a0ec 100644 --- a/pygmt/datatypes/header.py +++ b/pygmt/datatypes/header.py @@ -108,8 +108,10 @@ class _GMT_GRID_HEADER(ctp.Structure): # ruff: ignore[invalid-class-name] # Below are items used internally by GMT # Number of data points (n_columns * n_rows) [paddings are excluded] ("nm", ctp.c_size_t), - # Actual number of items (not bytes) required to hold this grid (mx * my), - # per band (for images) + # Actual number of items (not bytes) required to hold one layer of this grid, + # paddings included. size >= mx * my, rounded up to an even number. + # For images and cubes this is also the stride between bands/layers, i.e., + # band/layer k starts at offset k * size. ("size", ctp.c_size_t), # Bits per data value (e.g., 32 for ints/floats; 8 for bytes). # Only used for ERSI ArcInfo ASCII Exchange grids. @@ -121,7 +123,8 @@ class _GMT_GRID_HEADER(ctp.Structure): # ruff: ignore[invalid-class-name] ("complex_mode", ctp.c_uint), # Grid format ("type", ctp.c_uint), - # Number of bands [1]. Used with GMT_IMAGE containers + # Number of bands [1]. Used with GMT_IMAGE containers for the number of bands + # and with GMT_CUBE containers for the number of layers ("n_bands", ctp.c_uint), # Actual x-dimension in memory. mx = n_columns + pad[0] + pad[1] ("mx", ctp.c_uint), diff --git a/pygmt/tests/test_clib_read_data.py b/pygmt/tests/test_clib_read_data.py index e88ae83abf7..f76cf244193 100644 --- a/pygmt/tests/test_clib_read_data.py +++ b/pygmt/tests/test_clib_read_data.py @@ -5,6 +5,7 @@ from pathlib import Path import numpy as np +import numpy.testing as npt import pandas as pd import pytest import xarray as xr @@ -202,22 +203,94 @@ def test_clib_read_data_image_two_steps(expected_xrimage): xr.testing.assert_equal(xrimage, expected_xrimage) -def test_clib_read_data_cube_actual_grid(): +def test_clib_read_data_cube(): """ - Test the Session.read_data method for cube, but actually the file is a grid. + Test the Session.read_data method for cubes. """ + infile = which("@cube.nc", download="c") with Session() as lib: - with pytest.raises(GMTCLibError): - lib.read_data("@earth_relief_01d_p", kind="cube", mode="GMT_CONTAINER_ONLY") + cube = lib.read_data(infile, kind="cube").contents + header = cube.header.contents + assert header.n_rows == 11 + assert header.n_columns == 11 + assert header.n_bands == 4 # Number of layers in the cube + assert header.wesn[:] == [0.0, 10.0, 0.0, 10.0] + assert cube.z_range[:] == [1.0, 5.0] + + z = np.ctypeslib.as_array(cube.z, shape=(header.n_bands,)) + npt.assert_allclose(z, [1.0, 2.0, 3.0, 5.0]) + + # The cube data is a stack of 2-D padded layers, i.e., layer k starts at + # offset k * header.size. + data = np.ctypeslib.as_array(cube.data, shape=(header.n_bands, header.size)) + pad = header.pad[:] + data = data[:, : header.my * header.mx].reshape( + header.n_bands, header.my, header.mx + ) + data = data[:, pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1]] + assert data.shape == (4, 11, 11) + npt.assert_allclose(data.min(), 0.0) + npt.assert_allclose(data.max(), 140.0) + # The cube is X*Y scaled by 1, 1.1, 1.2 and 1.4 for the four layers, so the + # per-layer maxima also pin the layer order. + npt.assert_allclose( + [layer.max() for layer in data], [100.0, 110.0, 120.0, 140.0] + ) + # GMT stores rows north-first, so row 0 is y=10 (not the y=0 row of zeros). + npt.assert_allclose(data[0, 0], np.arange(0.0, 101.0, 10.0)) + + +def test_clib_read_data_cube_two_steps(): + """ + Test the Session.read_data method for cubes in two steps, first reading the header + and then the data. + """ + infile = which("@cube.nc", download="c") + with Session() as lib: + # Read the header first + data_ptr = lib.read_data(infile, kind="cube", mode="GMT_CONTAINER_ONLY") + cube = data_ptr.contents + header = cube.header.contents + assert header.n_rows == 11 + assert header.n_columns == 11 + assert header.n_bands == 4 + assert not cube.data # The data is not read yet + + # Read the data + lib.read_data(infile, kind="cube", mode="GMT_DATA_ONLY", data=data_ptr) + assert data_ptr.contents.data # The data is read now + + +def test_clib_read_data_cube_to_xarray(): + """ + Test converting a cube read by Session.read_data into an xarray.DataArray. + + The result is compared against the netCDF file read directly by xarray. + """ + infile = which("@cube.nc", download="c") + expected = xr.open_dataset(infile)["cube"] + with Session() as lib: + cube = lib.read_data(infile, kind="cube").contents.to_xarray() + + assert cube.dims == ("z", "y", "x") + assert cube.shape == (4, 11, 11) + # Coordinates are returned ascending, matching the CF conventions of the file. + for dim in cube.dims: + npt.assert_allclose(cube[dim], expected[dim]) + npt.assert_allclose(cube, expected) + assert cube.attrs["long_name"] == "cube" + assert cube.z.attrs["actual_range"].tolist() == [1.0, 5.0] -def test_clib_read_data_cube_actual_image(): +@pytest.mark.parametrize("infile", ["@earth_relief_01d_p", "@earth_day_01d"]) +def test_clib_read_data_cube_actual_grid_or_image(infile): """ - Test the Session.read_data method for cube, but actually the file is an image. + Test the Session.read_data method for cube, but actually the file is a grid or an + image. """ with Session() as lib: with pytest.raises(GMTCLibError): - lib.read_data("@earth_day_01d", kind="cube", mode="GMT_CONTAINER_ONLY") + lib.read_data(infile, kind="cube", mode="GMT_CONTAINER_ONLY") def test_clib_read_data_fails(): diff --git a/pygmt/tests/test_clib_virtualfiles.py b/pygmt/tests/test_clib_virtualfiles.py index 0a8b0abdfd2..c1e48ef3d08 100644 --- a/pygmt/tests/test_clib_virtualfiles.py +++ b/pygmt/tests/test_clib_virtualfiles.py @@ -1,15 +1,17 @@ """ -Test the Session.open_virtualfile method. +Test the Session methods for virtual files. """ from pathlib import Path import numpy as np +import numpy.testing as npt import pytest from pygmt import clib from pygmt.clib.session import DTYPES_NUMERIC from pygmt.exceptions import GMTCLibError, GMTValueError from pygmt.helpers import GMTTempFile +from pygmt.src import which from pygmt.tests.test_clib import mock POINTS_DATA = Path(__file__).parent / "data" / "points.txt" @@ -107,3 +109,20 @@ def test_open_virtualfile_bad_direction(): with pytest.raises(GMTValueError): with lib.open_virtualfile(*vfargs): pass + + +@pytest.mark.parametrize("kind", ["cube", None]) +def test_virtualfile_to_raster_cube(kind): + """ + Test Session.virtualfile_to_raster for cubes, with an explicit kind and with the + kind inquired from the virtual file. + """ + cubefile = which("@cube.nc", download="c") + with clib.Session() as lib: + with lib.virtualfile_out(kind="cube") as voutcube: + lib.call_module("read", [cubefile, voutcube, "-Tu"]) + cube = lib.virtualfile_to_raster(vfname=voutcube, kind=kind) + assert cube.dims == ("z", "y", "x") + assert cube.shape == (4, 11, 11) + npt.assert_allclose(cube.z, [1.0, 2.0, 3.0, 5.0]) + npt.assert_allclose(cube.max(dim=("y", "x")), [100.0, 110.0, 120.0, 140.0]) From f2b285ace1e616385c86b865d11c98f1aa05ec65 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 24 Aug 2026 18:20:18 +0800 Subject: [PATCH 15/22] Improve doctests --- pygmt/datatypes/cube.py | 32 ++++++++++++++++++++++++++------ 1 file changed, 26 insertions(+), 6 deletions(-) diff --git a/pygmt/datatypes/cube.py b/pygmt/datatypes/cube.py index e0d14193f97..06f491f4554 100644 --- a/pygmt/datatypes/cube.py +++ b/pygmt/datatypes/cube.py @@ -101,14 +101,17 @@ class _GMT_CUBE(ctp.Structure): # ruff: ignore[invalid-class-name] [ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]] >>> # The northernmost row of every layer. The four layers are the same X*Y grid - >>> # scaled by 1, 1.1, 1.2 and 1.4, so this also pins the layer order. + >>> # scaled by 1.0, 1.1, 1.2, and 1.4, respectively. >>> print(data[:, 0, :]) [[ 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. 100.] [ 0. 11. 22. 33. 44. 55. 66. 77. 88. 99. 110.] [ 0. 12. 24. 36. 48. 60. 72. 84. 96. 108. 120.] [ 0. 14. 28. 42. 56. 70. 84. 98. 112. 126. 140.]] - >>> print(data.max(axis=(1, 2))) - [100. 110. 120. 140.] + >>> # Verify that layer k equals scale[k] * outer(y, x) for every element. + >>> scale = [1.0, 1.1, 1.2, 1.4] + >>> expected = np.array([s * np.outer(y, x) for s in scale], dtype=np.float32) + >>> np.allclose(data, expected) + True """ _fields_: ClassVar = [ @@ -167,6 +170,7 @@ def to_xarray(self) -> xr.DataArray: Examples -------- + >>> import numpy as np >>> from pygmt import which >>> from pygmt.clib import Session >>> cubefile = which("@cube.nc", download="c") @@ -188,9 +192,25 @@ def to_xarray(self) -> xr.DataArray: long_name: z axis: Z actual_range: [1. 5.] - >>> # The four layers are the same X*Y grid scaled by 1, 1.1, 1.2 and 1.4. - >>> da.max(dim=("y", "x")).values - array([100., 110., 120., 140.], dtype=float32) + >>> # The four layers are the same X*Y grid scaled by 1.0, 1.1, 1.2, and 1.4. + >>> # Verify that layer k equals scale[k] * outer(y, x) for every element. + >>> scale = [1.0, 1.1, 1.2, 1.4] + >>> expected = [s * np.outer(da.y, da.x) for s in scale] + >>> np.allclose(da.values, expected) + True + >>> # Cross-check against loading the same file directly with xarray. + >>> import xarray as xr + >>> direct = xr.open_dataset(cubefile)["cube"] + >>> da.dims == direct.dims + True + >>> np.allclose(da.coords["x"], direct.coords["x"]) + True + >>> np.allclose(da.coords["y"], direct.coords["y"]) + True + >>> np.allclose(da.coords["z"], direct.coords["z"]) + True + >>> np.allclose(da.values, direct.values) + True >>> da.gmt.registration, da.gmt.gtype (, ) """ From 817846720ffe3c071d404b7deb9a3b574ca58161 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Tue, 25 Aug 2026 22:49:59 +0800 Subject: [PATCH 16/22] Improve description of n_bands --- pygmt/datatypes/header.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pygmt/datatypes/header.py b/pygmt/datatypes/header.py index 619c3e5a0ec..10c12903296 100644 --- a/pygmt/datatypes/header.py +++ b/pygmt/datatypes/header.py @@ -123,8 +123,7 @@ class _GMT_GRID_HEADER(ctp.Structure): # ruff: ignore[invalid-class-name] ("complex_mode", ctp.c_uint), # Grid format ("type", ctp.c_uint), - # Number of bands [1]. Used with GMT_IMAGE containers for the number of bands - # and with GMT_CUBE containers for the number of layers + # Number of bands for GMT_IMAGE or number of layers for GMT_CUBE [1]. ("n_bands", ctp.c_uint), # Actual x-dimension in memory. mx = n_columns + pad[0] + pad[1] ("mx", ctp.c_uint), From db98fba8ce6bd3d736bfca0557fea2f2bfa69582 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Tue, 25 Aug 2026 23:04:41 +0800 Subject: [PATCH 17/22] Revert changes in test_clib_virtualfiles.py --- pygmt/tests/test_clib_virtualfiles.py | 21 +-------------------- 1 file changed, 1 insertion(+), 20 deletions(-) diff --git a/pygmt/tests/test_clib_virtualfiles.py b/pygmt/tests/test_clib_virtualfiles.py index c1e48ef3d08..0a8b0abdfd2 100644 --- a/pygmt/tests/test_clib_virtualfiles.py +++ b/pygmt/tests/test_clib_virtualfiles.py @@ -1,17 +1,15 @@ """ -Test the Session methods for virtual files. +Test the Session.open_virtualfile method. """ from pathlib import Path import numpy as np -import numpy.testing as npt import pytest from pygmt import clib from pygmt.clib.session import DTYPES_NUMERIC from pygmt.exceptions import GMTCLibError, GMTValueError from pygmt.helpers import GMTTempFile -from pygmt.src import which from pygmt.tests.test_clib import mock POINTS_DATA = Path(__file__).parent / "data" / "points.txt" @@ -109,20 +107,3 @@ def test_open_virtualfile_bad_direction(): with pytest.raises(GMTValueError): with lib.open_virtualfile(*vfargs): pass - - -@pytest.mark.parametrize("kind", ["cube", None]) -def test_virtualfile_to_raster_cube(kind): - """ - Test Session.virtualfile_to_raster for cubes, with an explicit kind and with the - kind inquired from the virtual file. - """ - cubefile = which("@cube.nc", download="c") - with clib.Session() as lib: - with lib.virtualfile_out(kind="cube") as voutcube: - lib.call_module("read", [cubefile, voutcube, "-Tu"]) - cube = lib.virtualfile_to_raster(vfname=voutcube, kind=kind) - assert cube.dims == ("z", "y", "x") - assert cube.shape == (4, 11, 11) - npt.assert_allclose(cube.z, [1.0, 2.0, 3.0, 5.0]) - npt.assert_allclose(cube.max(dim=("y", "x")), [100.0, 110.0, 120.0, 140.0]) From 41e86a468d352aa5b03076d477779507d450d9e4 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Tue, 25 Aug 2026 23:12:46 +0800 Subject: [PATCH 18/22] Improve doctest for _GMT_CUBE --- pygmt/datatypes/cube.py | 25 ++++++++++--------------- 1 file changed, 10 insertions(+), 15 deletions(-) diff --git a/pygmt/datatypes/cube.py b/pygmt/datatypes/cube.py index 06f491f4554..cd57ec62d08 100644 --- a/pygmt/datatypes/cube.py +++ b/pygmt/datatypes/cube.py @@ -40,7 +40,7 @@ class _GMT_CUBE(ctp.Structure): # ruff: ignore[invalid-class-name] ... with lib.virtualfile_out(kind="cube") as voutcube: ... lib.call_module("read", [cubefile, voutcube, "-Tu"]) ... # Read the cube from the virtual file - ... cube = lib.read_virtualfile(vfname=voutcube, kind="cube").contents + ... cube = lib.read_virtualfile(voutcube, kind="cube").contents ... # The cube header ... header = cube.header.contents ... # Access the header properties @@ -64,6 +64,8 @@ class _GMT_CUBE(ctp.Structure): # ruff: ignore[invalid-class-name] ... ).copy() ... # Reshape the layers to 2-D and strip the paddings ... pad = header.pad[:] + ... # header.size can exceed header.my * header.mx because GMT rounds up to + ... # an even count. So we need to slice the data array to the actual size. ... data = data[:, : header.my * header.mx].reshape( ... header.n_bands, header.my, header.mx ... ) @@ -85,9 +87,7 @@ class _GMT_CUBE(ctp.Structure): # ruff: ignore[invalid-class-name] array([1., 2., 3., 5.]) >>> data.shape (4, 11, 11) - >>> print(data.min(), data.max()) - 0.0 140.0 - >>> # GMT stores rows north-first, so row 0 is y=10 and the last row is y=0. + >>> # The first layer is a grid of X * Y. >>> print(data[0, :, :]) [[ 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. 100.] [ 0. 9. 18. 27. 36. 45. 54. 63. 72. 81. 90.] @@ -100,17 +100,12 @@ class _GMT_CUBE(ctp.Structure): # ruff: ignore[invalid-class-name] [ 0. 2. 4. 6. 8. 10. 12. 14. 16. 18. 20.] [ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]] - >>> # The northernmost row of every layer. The four layers are the same X*Y grid - >>> # scaled by 1.0, 1.1, 1.2, and 1.4, respectively. - >>> print(data[:, 0, :]) - [[ 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. 100.] - [ 0. 11. 22. 33. 44. 55. 66. 77. 88. 99. 110.] - [ 0. 12. 24. 36. 48. 60. 72. 84. 96. 108. 120.] - [ 0. 14. 28. 42. 56. 70. 84. 98. 112. 126. 140.]] - >>> # Verify that layer k equals scale[k] * outer(y, x) for every element. - >>> scale = [1.0, 1.1, 1.2, 1.4] - >>> expected = np.array([s * np.outer(y, x) for s in scale], dtype=np.float32) - >>> np.allclose(data, expected) + >>> # Other layers are the same, but scaled by 1.0, 1.1, 1.2, and 1.4, respectively. + >>> np.allclose(data[1, :, :], 1.1 * data[0, :, :]) + True + >>> np.allclose(data[2, :, :], 1.2 * data[0, :, :]) + True + >>> np.allclose(data[3, :, :], 1.4 * data[0, :, :]) True """ From f3af7092c71520441875eddbf231486346f3f8c5 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Tue, 25 Aug 2026 23:38:13 +0800 Subject: [PATCH 19/22] Improve code and doctests for _GMT_CUBE.to_xarray --- pygmt/datatypes/cube.py | 57 +++++++++++++++++++++-------------------- 1 file changed, 29 insertions(+), 28 deletions(-) diff --git a/pygmt/datatypes/cube.py b/pygmt/datatypes/cube.py index cd57ec62d08..8d641fe3436 100644 --- a/pygmt/datatypes/cube.py +++ b/pygmt/datatypes/cube.py @@ -135,7 +135,7 @@ class _GMT_CUBE(ctp.Structure): # ruff: ignore[invalid-class-name] ("units", ctp.c_char * GMT_GRID_UNIT_LEN80), ] - def _parse_dimension(self) -> tuple[str, dict]: + def _parse_z_dimension(self) -> tuple[str, dict]: """ Get the name and attributes of the 3rd dimension. @@ -160,8 +160,7 @@ def to_xarray(self) -> xr.DataArray: Returns ------- dataarray - A 3-D :class:`xr.DataArray` object with dimensions ordered as - (3rd dimension, y, x). + A 3-D :class:`xr.DataArray` object with dimensions ordered as (z, y, x). Examples -------- @@ -178,42 +177,43 @@ def to_xarray(self) -> xr.DataArray: ... da = cube.contents.to_xarray() >>> da.name, da.dims, da.shape ('z', ('z', 'y', 'x'), (4, 11, 11)) + >>> da.coords["x"] + Size: 88B + array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]) + Coordinates: + * x (x) float64 88B 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 + Attributes: + long_name: x + axis: X + actual_range: [ 0. 10.] + >>> da.coords["y"] + Size: 88B + array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]) + Coordinates: + * y (y) float64 88B 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 + Attributes: + long_name: y + axis: Y + actual_range: [ 0. 10.] >>> da.coords["z"] Size: 32B array([1., 2., 3., 5.]) Coordinates: - * z (z) float64 32B 1.0 2.0 3.0 5.0 + * z (z) float64 32B 1.0 2.0 3.0 5.0 Attributes: long_name: z axis: Z actual_range: [1. 5.] - >>> # The four layers are the same X*Y grid scaled by 1.0, 1.1, 1.2, and 1.4. - >>> # Verify that layer k equals scale[k] * outer(y, x) for every element. - >>> scale = [1.0, 1.1, 1.2, 1.4] - >>> expected = [s * np.outer(da.y, da.x) for s in scale] - >>> np.allclose(da.values, expected) - True - >>> # Cross-check against loading the same file directly with xarray. - >>> import xarray as xr - >>> direct = xr.open_dataset(cubefile)["cube"] - >>> da.dims == direct.dims - True - >>> np.allclose(da.coords["x"], direct.coords["x"]) - True - >>> np.allclose(da.coords["y"], direct.coords["y"]) - True - >>> np.allclose(da.coords["z"], direct.coords["z"]) - True - >>> np.allclose(da.values, direct.values) - True >>> da.gmt.registration, da.gmt.gtype (, ) """ + # The cube header header = self.header.contents - # The y/x dimensions come from the 2-D grid header; the 3rd one from the cube. + # Get x/y dimensions and their attributes from the header. dims, dim_attrs = header.dims, header.dim_attrs - zdim, zdim_attrs = self._parse_dimension() + # Get the 3rd dimension and its attributes from the cube itself. + zdim, zdim_attrs = self._parse_z_dimension() # The coordinates, given as a tuple of the form (dims, data, attrs) x = np.ctypeslib.as_array(self.x, shape=(header.n_columns,)).copy() @@ -225,15 +225,16 @@ def to_xarray(self) -> xr.DataArray: (dims[1], x, dim_attrs[1]), ] - # The data array. The cube is a stack of 2-D padded layers, i.e., layer k - # starts at offset k * header.size, which can exceed header.my * header.mx. + # The cube is a stack of 2-D padded layers, i.e., layer k starts at offset + # k * header.size, which can exceed header.my * header.mx. data = np.ctypeslib.as_array( self.data, shape=(header.n_bands, header.size) ).copy() - pad = header.pad[:] data = data[:, : header.my * header.mx].reshape( header.n_bands, 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]] # Create the xarray.DataArray object From 7a1a03c4fe763d271e56c8a552e6cb58cfde6e7e Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Tue, 25 Aug 2026 23:49:31 +0800 Subject: [PATCH 20/22] Simplify tests --- pygmt/tests/test_clib_read_data.py | 74 ++++++++---------------------- 1 file changed, 20 insertions(+), 54 deletions(-) diff --git a/pygmt/tests/test_clib_read_data.py b/pygmt/tests/test_clib_read_data.py index f76cf244193..e0a0ba81fa3 100644 --- a/pygmt/tests/test_clib_read_data.py +++ b/pygmt/tests/test_clib_read_data.py @@ -5,7 +5,6 @@ from pathlib import Path import numpy as np -import numpy.testing as npt import pandas as pd import pytest import xarray as xr @@ -44,6 +43,14 @@ def fixture_expected_xrimage(): return None +@pytest.fixture(scope="module", name="expected_xrcube") +def fixture_expected_xrcube(): + """ + The expected xr.DataArray object for the @cube.nc file. + """ + return xr.load_dataarray(which("@cube.nc", download="c")) + + def test_clib_read_data_dataset(): """ Test the Session.read_data method for datasets. @@ -203,44 +210,18 @@ def test_clib_read_data_image_two_steps(expected_xrimage): xr.testing.assert_equal(xrimage, expected_xrimage) -def test_clib_read_data_cube(): +def test_clib_read_data_cube(expected_xrcube): """ Test the Session.read_data method for cubes. """ infile = which("@cube.nc", download="c") with Session() as lib: cube = lib.read_data(infile, kind="cube").contents - header = cube.header.contents - assert header.n_rows == 11 - assert header.n_columns == 11 - assert header.n_bands == 4 # Number of layers in the cube - assert header.wesn[:] == [0.0, 10.0, 0.0, 10.0] - assert cube.z_range[:] == [1.0, 5.0] - - z = np.ctypeslib.as_array(cube.z, shape=(header.n_bands,)) - npt.assert_allclose(z, [1.0, 2.0, 3.0, 5.0]) - - # The cube data is a stack of 2-D padded layers, i.e., layer k starts at - # offset k * header.size. - data = np.ctypeslib.as_array(cube.data, shape=(header.n_bands, header.size)) - pad = header.pad[:] - data = data[:, : header.my * header.mx].reshape( - header.n_bands, header.my, header.mx - ) - data = data[:, pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1]] - assert data.shape == (4, 11, 11) - npt.assert_allclose(data.min(), 0.0) - npt.assert_allclose(data.max(), 140.0) - # The cube is X*Y scaled by 1, 1.1, 1.2 and 1.4 for the four layers, so the - # per-layer maxima also pin the layer order. - npt.assert_allclose( - [layer.max() for layer in data], [100.0, 110.0, 120.0, 140.0] - ) - # GMT stores rows north-first, so row 0 is y=10 (not the y=0 row of zeros). - npt.assert_allclose(data[0, 0], np.arange(0.0, 101.0, 10.0)) - - -def test_clib_read_data_cube_two_steps(): + xrcube = cube.to_xarray() + xr.testing.assert_equal(xrcube, expected_xrcube) + + +def test_clib_read_data_cube_two_steps(expected_xrcube): """ Test the Session.read_data method for cubes in two steps, first reading the header and then the data. @@ -254,32 +235,17 @@ def test_clib_read_data_cube_two_steps(): assert header.n_rows == 11 assert header.n_columns == 11 assert header.n_bands == 4 + assert header.wesn[:] == [0.0, 10.0, 0.0, 10.0] + assert header.z_min == 0.0 + assert header.z_max == 140.0 assert not cube.data # The data is not read yet # Read the data lib.read_data(infile, kind="cube", mode="GMT_DATA_ONLY", data=data_ptr) - assert data_ptr.contents.data # The data is read now - - -def test_clib_read_data_cube_to_xarray(): - """ - Test converting a cube read by Session.read_data into an xarray.DataArray. - The result is compared against the netCDF file read directly by xarray. - """ - infile = which("@cube.nc", download="c") - expected = xr.open_dataset(infile)["cube"] - with Session() as lib: - cube = lib.read_data(infile, kind="cube").contents.to_xarray() - - assert cube.dims == ("z", "y", "x") - assert cube.shape == (4, 11, 11) - # Coordinates are returned ascending, matching the CF conventions of the file. - for dim in cube.dims: - npt.assert_allclose(cube[dim], expected[dim]) - npt.assert_allclose(cube, expected) - assert cube.attrs["long_name"] == "cube" - assert cube.z.attrs["actual_range"].tolist() == [1.0, 5.0] + # Full check + xrcube = data_ptr.contents.to_xarray() + xr.testing.assert_equal(xrcube, expected_xrcube) @pytest.mark.parametrize("infile", ["@earth_relief_01d_p", "@earth_day_01d"]) From 80388cf6ec4cbc2def832b5a6ac909db58d9555d Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Wed, 26 Aug 2026 00:30:18 +0800 Subject: [PATCH 21/22] Extending xarray backend --- pygmt/datatypes/cube.py | 7 ++++--- pygmt/tests/test_clib_read_data.py | 4 +++- pygmt/xarray/backend.py | 16 +++++++++++----- 3 files changed, 18 insertions(+), 9 deletions(-) diff --git a/pygmt/datatypes/cube.py b/pygmt/datatypes/cube.py index 8d641fe3436..e5b3f60dc23 100644 --- a/pygmt/datatypes/cube.py +++ b/pygmt/datatypes/cube.py @@ -176,7 +176,7 @@ def to_xarray(self) -> xr.DataArray: ... # Convert to xarray.DataArray and use it later ... da = cube.contents.to_xarray() >>> da.name, da.dims, da.shape - ('z', ('z', 'y', 'x'), (4, 11, 11)) + ('cube', ('z', 'y', 'x'), (4, 11, 11)) >>> da.coords["x"] Size: 88B array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]) @@ -237,9 +237,10 @@ def to_xarray(self) -> xr.DataArray: pad = header.pad[:] data = data[:, pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1]] - # Create the xarray.DataArray object + # Create the xarray.DataArray object. + # The cube name is stored in the header's z_units attribute. cube = xr.DataArray( - data, coords=coords, name=header.name, attrs=header.data_attrs + data, coords=coords, name=header.z_units, attrs=header.data_attrs ) # Flip the coordinates and data if necessary so that coordinates are ascending. diff --git a/pygmt/tests/test_clib_read_data.py b/pygmt/tests/test_clib_read_data.py index e0a0ba81fa3..ea189b79cf9 100644 --- a/pygmt/tests/test_clib_read_data.py +++ b/pygmt/tests/test_clib_read_data.py @@ -48,7 +48,9 @@ def fixture_expected_xrcube(): """ The expected xr.DataArray object for the @cube.nc file. """ - return xr.load_dataarray(which("@cube.nc", download="c")) + return xr.load_dataarray( + which("@cube.nc", download="c"), engine="gmt", raster_kind="cube" + ) def test_clib_read_data_dataset(): diff --git a/pygmt/xarray/backend.py b/pygmt/xarray/backend.py index cc5a0f6322c..a1ceb9582d7 100644 --- a/pygmt/xarray/backend.py +++ b/pygmt/xarray/backend.py @@ -115,7 +115,7 @@ def open_dataset( # type: ignore[override] filename_or_obj: PathLike, *, drop_variables=None, # ruff: ignore[unused-method-argument] - raster_kind: Literal["grid", "image"], + raster_kind: Literal["grid", "image", "cube"], region: Sequence[float] | str | None = None, # other backend specific keyword arguments # `chunks` and `cache` DO NOT go here, they are handled by xarray @@ -130,19 +130,25 @@ def open_dataset( # type: ignore[override] that can be read by GMT via the netCDF or GDAL C libraries. See also :gmt-docs:`reference/features.html#grid-file-format`. raster_kind - Whether to read the file as a "grid" (single-band) or "image" (multi-band). + Whether to read the file as a "grid" (single-band), "image" (multi-band), or + "cube" (stacks of 2-D grids). region The subregion of the grid or image to load, in the form of a sequence [*xmin*, *xmax*, *ymin*, *ymax*] or an ISO country code. """ - if raster_kind not in {"grid", "image"}: + if raster_kind not in {"grid", "image", "cube"}: raise GMTValueError( - raster_kind, description="raster kind", choices=["grid", "image"] + raster_kind, + description="raster kind", + choices=["grid", "image", "cube"], ) with Session() as lib: with lib.virtualfile_out(kind=raster_kind) as voutfile: - kwdict = {"R": region, "T": {"grid": "g", "image": "i"}[raster_kind]} + kwdict = { + "R": region, + "T": {"grid": "g", "image": "i", "cube": "u"}[raster_kind], + } lib.call_module( module="read", args=[filename_or_obj, voutfile, *build_arg_list(kwdict)], From e1ceda2204097b584dbeb63ac18c119669306585 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Wed, 26 Aug 2026 00:49:02 +0800 Subject: [PATCH 22/22] Fix z_unit --- pygmt/datatypes/cube.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygmt/datatypes/cube.py b/pygmt/datatypes/cube.py index e5b3f60dc23..27b65373f7a 100644 --- a/pygmt/datatypes/cube.py +++ b/pygmt/datatypes/cube.py @@ -240,7 +240,7 @@ def to_xarray(self) -> xr.DataArray: # Create the xarray.DataArray object. # The cube name is stored in the header's z_units attribute. cube = xr.DataArray( - data, coords=coords, name=header.z_units, attrs=header.data_attrs + data, coords=coords, name=header.z_units.decode(), attrs=header.data_attrs ) # Flip the coordinates and data if necessary so that coordinates are ascending.