From bcf43f0f9e8bc45137129bd48662386180d4d7b6 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 18 Mar 2024 22:34:50 +0800 Subject: [PATCH 01/34] Wrap GMT's standard data type GMT_IMAGE for images --- pygmt/clib/session.py | 22 ++++++++----- pygmt/datatypes/__init__.py | 1 + pygmt/datatypes/image.py | 66 +++++++++++++++++++++++++++++++++++++ 3 files changed, 81 insertions(+), 8 deletions(-) create mode 100644 pygmt/datatypes/image.py diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 6e6e495e244..bfea38a5bf1 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -23,7 +23,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, @@ -1615,7 +1615,9 @@ def virtualfile_in( # noqa: PLR0912 @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. @@ -1628,8 +1630,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. @@ -1672,12 +1674,15 @@ 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 read_virtualfile( - self, vfname: str, kind: Literal["dataset", "grid", None] = None + self, vfname: str, kind: Literal["dataset", "grid", "image", None] = None ): """ Read data from a virtual file and optionally cast into a GMT data container. @@ -1688,7 +1693,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 -------- @@ -1735,7 +1741,7 @@ def read_virtualfile( # _GMT_DATASET). if kind is None: # Return the ctypes void pointer return pointer - 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( diff --git a/pygmt/datatypes/__init__.py b/pygmt/datatypes/__init__.py index 237a050a9f7..3489dd19d10 100644 --- a/pygmt/datatypes/__init__.py +++ b/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 diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py new file mode 100644 index 00000000000..74af128dc3c --- /dev/null +++ b/pygmt/datatypes/image.py @@ -0,0 +1,66 @@ +""" +Wrapper for the GMT_IMAGE data type. +""" + +import ctypes as ctp +from typing import ClassVar + +from pygmt.datatypes.grid 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] + >>> data = xr.DataArray( + ... data, dims=["y", "x", "band"], coords={"y": y, "x": x, "band": [1, 2, 3]} + ... ) + >>> data = data.transpose("band", "y", "x") + >>> data = data.sortby(list(data.dims)) + >>> data.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), + ] From a052a1af1e7a3118dcf7e6d7e3d28bf236b2e4d0 Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Wed, 20 Mar 2024 20:38:57 +1300 Subject: [PATCH 02/34] Initial implementation of to_dataarray method for _GMT_IMAGE class Minimal xarray.DataArray output with data and coordinates, no metadata yet. --- pygmt/datatypes/image.py | 46 +++++++++++++++++++++++++++++++++++----- 1 file changed, 41 insertions(+), 5 deletions(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index 74af128dc3c..fcd052ee695 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -5,6 +5,8 @@ import ctypes as ctp from typing import ClassVar +import numpy as np +import xarray as xr from pygmt.datatypes.grid import _GMT_GRID_HEADER @@ -34,12 +36,14 @@ class _GMT_IMAGE(ctp.Structure): # noqa: N801 ... 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] - >>> data = xr.DataArray( - ... data, dims=["y", "x", "band"], coords={"y": y, "x": x, "band": [1, 2, 3]} + >>> da = xr.DataArray( + ... data=data, + ... dims=["y", "x", "band"], + ... coords={"y": y, "x": x, "band": [1, 2, 3]}, ... ) - >>> data = data.transpose("band", "y", "x") - >>> data = data.sortby(list(data.dims)) - >>> data.plot.imshow() + >>> da = da.transpose("band", "y", "x") + >>> da = da.sortby(list(data.dims)) + >>> da.plot.imshow() """ _fields_: ClassVar = [ @@ -64,3 +68,35 @@ class _GMT_IMAGE(ctp.Structure): # noqa: N801 # Book-keeping variables "hidden" from the API ("hidden", ctp.c_void_p), ] + + def to_dataarray(self): + """ + Convert a _GMT_GRID object to an :class:`xarray.DataArray` object. + + Returns + ------- + dataarray + A :class:`xarray.DataArray` object. + """ + + # Get grid 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] = { + "x": self.x[: header.n_columns], + "y": self.y[: header.n_rows], + "band": np.arange(stop=3, dtype=np.uint8), + } + + # Create the xarray.DataArray object + image = xr.DataArray(data=data, coords=coords, dims=("y", "x", "band")) + + return image From 4cce4a2c4fa404caebb21581a496eb6f0326bbd1 Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Tue, 18 Jun 2024 21:29:08 +1200 Subject: [PATCH 03/34] Small typo fixes and add output type-hint for to_dataarray --- pygmt/datatypes/image.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index fcd052ee695..d1e59ba9bf2 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -69,9 +69,9 @@ class _GMT_IMAGE(ctp.Structure): # noqa: N801 ("hidden", ctp.c_void_p), ] - def to_dataarray(self): + def to_dataarray(self) -> xr.DataArray: """ - Convert a _GMT_GRID object to an :class:`xarray.DataArray` object. + Convert a _GMT_IMAGE object to an :class:`xarray.DataArray` object. Returns ------- @@ -79,7 +79,7 @@ def to_dataarray(self): A :class:`xarray.DataArray` object. """ - # Get grid header + # Get image header header: _GMT_GRID_HEADER = self.header.contents # Get DataArray without padding From e02b65016a74524436e14294a9b641b2cfcdee32 Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Tue, 18 Jun 2024 21:47:56 +1200 Subject: [PATCH 04/34] Fix mypy error using np.array([0, 1, 2]) instead of np.arange --- pygmt/datatypes/image.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index d1e59ba9bf2..3cfd37e5cc8 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -90,10 +90,10 @@ def to_dataarray(self) -> xr.DataArray: )[pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1], :] # Get x and y coordinates - coords: dict[str, list] = { + coords: dict[str, list | np.ndarray] = { "x": self.x[: header.n_columns], "y": self.y[: header.n_rows], - "band": np.arange(stop=3, dtype=np.uint8), + "band": np.array([0, 1, 2], dtype=np.uint8), } # Create the xarray.DataArray object From f3d4b1fec79ee9f791088a0aba0208be1fbb4a5f Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Tue, 18 Jun 2024 22:08:32 +1200 Subject: [PATCH 05/34] Parse name and data_attrs from grid/image header Extra metadata from the _GMT_GRID_HEADER struct. --- pygmt/datatypes/image.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index 3cfd37e5cc8..6d425480890 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -7,7 +7,7 @@ import numpy as np import xarray as xr -from pygmt.datatypes.grid import _GMT_GRID_HEADER +from pygmt.datatypes.header import _GMT_GRID_HEADER class _GMT_IMAGE(ctp.Structure): # noqa: N801 @@ -97,6 +97,12 @@ def to_dataarray(self) -> xr.DataArray: } # Create the xarray.DataArray object - image = xr.DataArray(data=data, coords=coords, dims=("y", "x", "band")) + image = xr.DataArray( + data=data, + coords=coords, + dims=("y", "x", "band"), + name=header.name, + attrs=header.data_attrs, + ) return image From 43901367f93027d745046385dda909a41d80a781 Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Thu, 20 Jun 2024 21:36:10 +1200 Subject: [PATCH 06/34] Transpose array to (band, y, x) order and add doctest for to_dataarray Reorder the dimensions to follow Channel, Height, Width (CHW) convention. Also added doctest checking output DataArray object and the image's x and y coordinates. --- pygmt/datatypes/image.py | 68 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 67 insertions(+), 1 deletion(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index 6d425480890..c86017945c9 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -77,6 +77,72 @@ def to_dataarray(self) -> xr.DataArray: ------- dataarray A :class:`xarray.DataArray` object. + + Examples + -------- + >>> from pygmt.clib import Session + >>> with Session() as lib: + ... with lib.virtualfile_out(kind="image") as voutimg: + ... lib.call_module("read", ["@earth_day_01d", voutimg, "-Ti"]) + ... # Read the image from the virtual file + ... image = lib.read_virtualfile(voutimg, kind="image") + ... # Convert to xarray.DataArray and use it later + ... da = image.contents.to_dataarray() + >>> da # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS + Size: 2MB + array([[[ 10, 10, 10, ..., 10, 10, 10], + [ 10, 10, 10, ..., 10, 10, 10], + [ 10, 10, 10, ..., 10, 10, 10], + ..., + [192, 193, 193, ..., 193, 192, 191], + [204, 206, 206, ..., 205, 206, 204], + [208, 210, 210, ..., 210, 210, 208]], + + [[ 10, 10, 10, ..., 10, 10, 10], + [ 10, 10, 10, ..., 10, 10, 10], + [ 10, 10, 10, ..., 10, 10, 10], + ..., + [186, 187, 188, ..., 187, 186, 185], + [196, 198, 198, ..., 197, 197, 196], + [199, 201, 201, ..., 201, 202, 199]], + + [[ 51, 51, 51, ..., 51, 51, 51], + [ 51, 51, 51, ..., 51, 51, 51], + [ 51, 51, 51, ..., 51, 51, 51], + ..., + [177, 179, 179, ..., 178, 177, 177], + [185, 187, 187, ..., 187, 186, 185], + [189, 191, 191, ..., 191, 191, 189]]]) + Coordinates: + * x (x) float64 3kB -179.5 -178.5 -177.5 -176.5 ... 177.5 178.5 179.5 + * y (y) float64 1kB 89.5 88.5 87.5 86.5 ... -86.5 -87.5 -88.5 -89.5 + * band (band) uint8 3B 0 1 2 + Attributes: + title: + history: + description: + long_name: z + actual_range: [ 1.79769313e+308 -1.79769313e+308] + + >>> da.coords["x"] # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS + Size: 3kB + array([-179.5, -178.5, -177.5, ..., 177.5, 178.5, 179.5]) + Coordinates: + * x (x) float64 3kB -179.5 -178.5 -177.5 -176.5 ... 177.5 178.5 179.5 + + >>> da.coords["y"] # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS + Size: 1kB + array([ 89.5, 88.5, 87.5, 86.5, 85.5, 84.5, 83.5, 82.5, 81.5, 80.5, + 79.5, 78.5, 77.5, 76.5, 75.5, 74.5, 73.5, 72.5, 71.5, 70.5, + 69.5, 68.5, 67.5, 66.5, 65.5, 64.5, 63.5, 62.5, 61.5, 60.5, + ... + -0.5, -1.5, -2.5, -3.5, -4.5, -5.5, -6.5, -7.5, -8.5, -9.5, + ... + -60.5, -61.5, -62.5, -63.5, -64.5, -65.5, -66.5, -67.5, -68.5, -69.5, + -70.5, -71.5, -72.5, -73.5, -74.5, -75.5, -76.5, -77.5, -78.5, -79.5, + -80.5, -81.5, -82.5, -83.5, -84.5, -85.5, -86.5, -87.5, -88.5, -89.5]) + Coordinates: + * y (y) float64 1kB 89.5 88.5 87.5 86.5 ... -86.5 -87.5 -88.5 -89.5 """ # Get image header @@ -103,6 +169,6 @@ def to_dataarray(self) -> xr.DataArray: dims=("y", "x", "band"), name=header.name, attrs=header.data_attrs, - ) + ).transpose("band", "y", "x") return image From 5f25669583133e727a289c33eb853712236ebebe Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Thu, 20 Jun 2024 21:38:16 +1200 Subject: [PATCH 07/34] Set registration and gtype from header Get the registration and gtype info from the grid header and apply it to the GMT accessor attributes. --- pygmt/datatypes/image.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index c86017945c9..7db5e15b49c 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -143,6 +143,9 @@ def to_dataarray(self) -> xr.DataArray: -80.5, -81.5, -82.5, -83.5, -84.5, -85.5, -86.5, -87.5, -88.5, -89.5]) Coordinates: * y (y) float64 1kB 89.5 88.5 87.5 86.5 ... -86.5 -87.5 -88.5 -89.5 + + >>> da.gmt.registration, da.gmt.gtype + (1, 0) """ # Get image header @@ -171,4 +174,8 @@ def to_dataarray(self) -> xr.DataArray: attrs=header.data_attrs, ).transpose("band", "y", "x") + # Set GMT accessors. + # Must put at the end, otherwise info gets lost after certain image operations. + image.gmt.registration = header.registration + image.gmt.gtype = header.gtype return image From a3c6c14098d4e2f50152ac955beae942f048a3c0 Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Thu, 20 Jun 2024 21:45:45 +1200 Subject: [PATCH 08/34] Print basic shape and padding info in _GMT_IMAGE doctest Trying to match some of the doctests in _GMT_GRID. --- pygmt/datatypes/image.py | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index 7db5e15b49c..b99e9694885 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -24,26 +24,28 @@ class _GMT_IMAGE(ctp.Structure): # noqa: N801 >>> 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) + ... # Read the image from the virtual file + ... image = lib.read_virtualfile(vfname=voutimg, kind="image").contents + ... # The image header + ... header = image.header.contents + ... # Access the header properties + ... print(image.type, header.n_bands, header.n_rows, header.n_columns) ... print(header.pad[:]) + ... # The x and y coordinates + ... x = image.x[: header.n_columns] + ... y = image.y[: header.n_rows] + ... # The data array (with paddings) ... data = np.reshape( - ... ds.data[: header.n_bands * header.mx * header.my], + ... image.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], :] - ... 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() + ... print(data.shape) + 1 3 180 360 + [2, 2, 2, 2] + (180, 360, 3) """ _fields_: ClassVar = [ From 5888e105a25063889d413ed65a6b694152458348 Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Thu, 20 Jun 2024 21:50:05 +1200 Subject: [PATCH 09/34] Only set Conventions = CF-1.7 attribute for NetCDF grid type Remove hardcoded attribute that was only meant for NetCDF files, so that GeoTIFF files won't have it. --- pygmt/datatypes/header.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pygmt/datatypes/header.py b/pygmt/datatypes/header.py index 04e10ac0c72..ab109521131 100644 --- a/pygmt/datatypes/header.py +++ b/pygmt/datatypes/header.py @@ -203,7 +203,8 @@ def data_attrs(self) -> dict[str, Any]: Attributes for the data variable from the grid header. """ attrs: dict[str, Any] = {} - attrs["Conventions"] = "CF-1.7" + if self.type == 18: # Grid file format: ns = GMT netCDF format + attrs["Conventions"] = "CF-1.7" attrs["title"] = self.title.decode() attrs["history"] = self.command.decode() attrs["description"] = self.remark.decode() From 3dbf2f24e277784d232c89bbcb0b37e36300bf37 Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Thu, 20 Jun 2024 22:01:16 +1200 Subject: [PATCH 10/34] Remove rioxarray import --- pygmt/datatypes/image.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index b99e9694885..ded7c818ae3 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -19,7 +19,6 @@ class _GMT_IMAGE(ctp.Structure): # noqa: N801 >>> 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: From 7d437be0675d7a986bafa937c8fe19557d7bcce6 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 00:34:58 +0800 Subject: [PATCH 11/34] Use enum for grid ids --- pygmt/datatypes/header.py | 43 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 42 insertions(+), 1 deletion(-) diff --git a/pygmt/datatypes/header.py b/pygmt/datatypes/header.py index ab109521131..db1428036bc 100644 --- a/pygmt/datatypes/header.py +++ b/pygmt/datatypes/header.py @@ -3,6 +3,7 @@ """ import ctypes as ctp +from enum import Enum from typing import Any, ClassVar import numpy as np @@ -23,6 +24,40 @@ gmt_grdfloat = ctp.c_float +class GMTGridID(Enum): + """ + Enum for the GMT grid format ID. + + Defined in gmt_grdio.h. + """ + + GMT_GRD_UNKNOWN_FMT = 0 # if grid format cannot be auto-detected + GMT_GRID_IS_BF = 1 # GMT native, C-binary format (32-bit float) + GMT_GRID_IS_BS = 2 # GMT native, C-binary format (16-bit integer) + GMT_GRID_IS_RB = 3 # SUN rasterfile format (8-bit standard) + GMT_GRID_IS_BB = 4 # GMT native, C-binary format (8-bit integer) + GMT_GRID_IS_BM = 5 # GMT native, C-binary format (bit-mask) + GMT_GRID_IS_SF = 6 # Golden Software Surfer format 6 (32-bit float) + GMT_GRID_IS_CB = 7 # GMT netCDF format (8-bit integer) + GMT_GRID_IS_CS = 8 # GMT netCDF format (16-bit integer) + GMT_GRID_IS_CI = 9 # GMT netCDF format (32-bit integer) + GMT_GRID_IS_CF = 10 # GMT netCDF format (32-bit float) + GMT_GRID_IS_CD = 11 # GMT netCDF format (64-bit float) + GMT_GRID_IS_RF = 12 # GEODAS grid format GRD98 (NGDC) + GMT_GRID_IS_BI = 13 # GMT native, C-binary format (32-bit integer) + GMT_GRID_IS_BD = 14 # GMT native, C-binary format (64-bit float) + GMT_GRID_IS_NB = 15 # GMT netCDF format (8-bit integer) + GMT_GRID_IS_NS = 16 # GMT netCDF format (16-bit integer) + GMT_GRID_IS_NI = 17 # GMT netCDF format (32-bit integer) + GMT_GRID_IS_NF = 18 # GMT netCDF format (32-bit float) + GMT_GRID_IS_ND = 19 # GMT netCDF format (64-bit float) + GMT_GRID_IS_SD = 20 # Golden Software Surfer format 7 (64-bit float, read-only) + GMT_GRID_IS_AF = 21 # Atlantic Geoscience Center format AGC (32-bit float) + GMT_GRID_IS_GD = 22 # Import through GDAL + GMT_GRID_IS_EI = 23 # ESRI Arc/Info ASCII Grid Interchange format (ASCII integer) + GMT_GRID_IS_EF = 24 # ESRI Arc/Info ASCII Grid Interchange format (ASCII float) + + def _parse_nameunits(nameunits: str) -> tuple[str, str | None]: """ Get the long_name and units attributes from x_units/y_units/z_units in the grid @@ -203,7 +238,13 @@ def data_attrs(self) -> dict[str, Any]: Attributes for the data variable from the grid header. """ attrs: dict[str, Any] = {} - if self.type == 18: # Grid file format: ns = GMT netCDF format + if self.type in { + GMTGridID.GMT_GRID_IS_NB, + GMTGridID.GMT_GRID_IS_NS, + GMTGridID.GMT_GRID_IS_NI, + GMTGridID.GMT_GRID_IS_NF, + GMTGridID.GMT_GRID_IS_ND, + }: attrs["Conventions"] = "CF-1.7" attrs["title"] = self.title.decode() attrs["history"] = self.command.decode() From 268e34e61086073f5f47ab5b406ccc7adce110ad Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 00:53:24 +0800 Subject: [PATCH 12/34] Fix the band. Starting from 1 --- pygmt/datatypes/image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index c36960f2a6a..4451eb83666 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -188,7 +188,7 @@ def to_dataarray(self) -> xr.DataArray: 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), + "band": np.array([1, 2, 3], dtype=np.uint8), } # Create the xarray.DataArray object From 86765e1492e1a878aa871119e897515396972942 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 00:59:30 +0800 Subject: [PATCH 13/34] Refactor the tests for images --- pygmt/tests/test_clib_read_data.py | 36 ++++++++++++++++-------------- 1 file changed, 19 insertions(+), 17 deletions(-) diff --git a/pygmt/tests/test_clib_read_data.py b/pygmt/tests/test_clib_read_data.py index b2f92403df0..67f297fab2c 100644 --- a/pygmt/tests/test_clib_read_data.py +++ b/pygmt/tests/test_clib_read_data.py @@ -14,7 +14,7 @@ from pygmt.src import which try: - import rioxarray # noqa: F401 + import rioxarray _HAS_RIOXARRAY = True except ImportError: @@ -29,6 +29,16 @@ def fixture_expected_xrgrid(): return load_dataarray(which("@static_earth_relief.nc")) +@pytest.fixture(scope="module", name="expected_xrimage") +def fixture_expected_xrimage(): + """ + The expected xr.DataArray object for the @earth_day_01d_p file. + """ + with rioxarray.open_rasterio(which("@earth_day_01d_p")) as da: + dataarray = da.load().drop_vars("spatial_ref") + return dataarray + + def test_clib_read_data_dataset(): """ Test the Session.read_data method for datasets. @@ -102,7 +112,7 @@ def test_clib_read_data_grid_two_steps(expected_xrgrid): xr.testing.assert_equal(xrgrid, expected_xrgrid) -def test_clib_read_data_grid_actual_image(): +def test_clib_read_data_grid_actual_image(expected_xrimage): """ Test the Session.read_data method for grid, but actually the file is an image. """ @@ -120,34 +130,25 @@ def test_clib_read_data_grid_actual_image(): if _HAS_RIOXARRAY: # Full check if rioxarray is installed. xrimage = image.to_dataarray() - expected_xrimage = xr.open_dataarray( - which("@earth_day_01d_p"), engine="rasterio" - ) assert expected_xrimage.band.size == 3 # 3-band image. xr.testing.assert_equal( xrimage, - expected_xrimage.isel(band=0) - .drop_vars(["band", "spatial_ref"]) - .sortby("y"), + expected_xrimage.isel(band=0).drop_vars(["band"]).sortby("y"), ) # Note: Simplify the tests for images after GMT_IMAGE.to_dataarray() is implemented. -def test_clib_read_data_image(): +def test_clib_read_data_image(expected_xrimage): """ Test the Session.read_data method for images. """ with Session() as lib: image = lib.read_data("@earth_day_01d_p", kind="image").contents - header = image.header.contents - assert header.n_rows == 180 - assert header.n_columns == 360 - assert header.n_bands == 3 - assert header.wesn[:] == [-180.0, 180.0, -90.0, 90.0] - assert image.data + xrimage = image.to_dataarray() + xr.testing.assert_equal(xrimage, expected_xrimage) -def test_clib_read_data_image_two_steps(): +def test_clib_read_data_image_two_steps(expected_xrimage): """ Test the Session.read_data method for images in two steps, first reading the header and then the data. @@ -166,7 +167,8 @@ def test_clib_read_data_image_two_steps(): # Read the data lib.read_data(infile, kind="image", mode="GMT_DATA_ONLY", data=data_ptr) - assert image.data + xrimage = data_ptr.contents.to_dataarray() + xr.testing.assert_equal(xrimage, expected_xrimage) def test_clib_read_data_fails(): From 86f3ffa2e2e97a628dcc597197eae1f3b61427cf Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 10:56:51 +0800 Subject: [PATCH 14/34] In np.reshape, a is a position-only parameter --- pygmt/datatypes/image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index 4451eb83666..d1f5592d1fa 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -180,7 +180,7 @@ def to_dataarray(self) -> xr.DataArray: # Get DataArray without padding pad = header.pad[:] data: np.ndarray = np.reshape( - a=self.data[: header.n_bands * header.mx * header.my], + 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], :] From cc282470e7c379daad9b910d85fe1b1952183b29 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 12:17:03 +0800 Subject: [PATCH 15/34] Improve tests --- pygmt/tests/test_clib_read_data.py | 31 +++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/pygmt/tests/test_clib_read_data.py b/pygmt/tests/test_clib_read_data.py index 67f297fab2c..5f3ad68fd83 100644 --- a/pygmt/tests/test_clib_read_data.py +++ b/pygmt/tests/test_clib_read_data.py @@ -108,6 +108,8 @@ def test_clib_read_data_grid_two_steps(expected_xrgrid): # Read the data lib.read_data(infile, kind="grid", mode="GMT_DATA_ONLY", data=data_ptr) + + # Full check xrgrid = data_ptr.contents.to_dataarray() xr.testing.assert_equal(xrgrid, expected_xrgrid) @@ -117,10 +119,7 @@ def test_clib_read_data_grid_actual_image(expected_xrimage): Test the Session.read_data method for grid, but actually the file is an image. """ with Session() as lib: - data_ptr = lib.read_data( - "@earth_day_01d_p", kind="grid", mode="GMT_CONTAINER_AND_DATA" - ) - image = data_ptr.contents + image = lib.read_data("@earth_day_01d_p", kind="grid").contents header = image.header.contents assert header.n_rows == 180 assert header.n_columns == 360 @@ -128,24 +127,34 @@ def test_clib_read_data_grid_actual_image(expected_xrimage): # Explicitly check n_bands. Only one band is read for 3-band images. assert header.n_bands == 1 + xrimage = image.to_dataarray() if _HAS_RIOXARRAY: # Full check if rioxarray is installed. - xrimage = image.to_dataarray() assert expected_xrimage.band.size == 3 # 3-band image. xr.testing.assert_equal( xrimage, expected_xrimage.isel(band=0).drop_vars(["band"]).sortby("y"), ) + else: + assert xrimage.shape == (180, 360) -# Note: Simplify the tests for images after GMT_IMAGE.to_dataarray() is implemented. def test_clib_read_data_image(expected_xrimage): """ Test the Session.read_data method for images. """ with Session() as lib: image = lib.read_data("@earth_day_01d_p", kind="image").contents + header = image.header.contents + assert header.n_rows == 180 + assert header.n_columns == 360 + assert header.n_bands == 3 + assert header.wesn[:] == [-180.0, 180.0, -90.0, 90.0] + xrimage = image.to_dataarray() - xr.testing.assert_equal(xrimage, expected_xrimage) + if _HAS_RIOXARRAY: # Full check if rioxarray is installed. + xr.testing.assert_equal(xrimage, expected_xrimage) + else: + assert xrimage.shape == (3, 180, 360) def test_clib_read_data_image_two_steps(expected_xrimage): @@ -167,8 +176,12 @@ def test_clib_read_data_image_two_steps(expected_xrimage): # Read the data lib.read_data(infile, kind="image", mode="GMT_DATA_ONLY", data=data_ptr) - xrimage = data_ptr.contents.to_dataarray() - xr.testing.assert_equal(xrimage, expected_xrimage) + + xrimage = image.to_dataarray() + if _HAS_RIOXARRAY: # Full check if rioxarray is installed. + xr.testing.assert_equal(xrimage, expected_xrimage) + else: + assert xrimage.shape == (3, 180, 360) def test_clib_read_data_fails(): From 1e2c973976fac0e73d967931bd52cf06cedc8f93 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 14:46:51 +0800 Subject: [PATCH 16/34] Fix one failing doctest due to xarray changes --- pygmt/datatypes/image.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index d1f5592d1fa..b6a2ef5d04e 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -140,9 +140,9 @@ def to_dataarray(self) -> xr.DataArray: [185, 187, 187, ..., 187, 186, 185], [189, 191, 191, ..., 191, 191, 189]]]) Coordinates: - * x (x) float64 3kB -179.5 -178.5 -177.5 -176.5 ... 177.5 178.5 179.5 - * y (y) float64 1kB 89.5 88.5 87.5 86.5 ... -86.5 -87.5 -88.5 -89.5 - * band (band) uint8 3B 0 1 2 + * x (x) float64... -179.5 -178.5 -177.5 -176.5 ... 177.5 178.5 179.5 + * y (y) float64... 89.5 88.5 87.5 86.5 ... -86.5 -87.5 -88.5 -89.5 + * band (band) uint8... 1 2 3 Attributes: title: history: From 734dc28185fcccd856ba93e23e92dab40b274669 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 14:47:32 +0800 Subject: [PATCH 17/34] The np.reshape's newshape parameter is deprecated --- pygmt/datatypes/image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index b6a2ef5d04e..4bbdc352968 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -181,7 +181,7 @@ def to_dataarray(self) -> xr.DataArray: pad = header.pad[:] data: np.ndarray = np.reshape( self.data[: header.n_bands * header.mx * header.my], - newshape=(header.my, header.mx, header.n_bands), + shape=(header.my, header.mx, header.n_bands), )[pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1], :] # Get x and y coordinates From 919dc0065955da884e9815210d241da88f1e5259 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 14:57:57 +0800 Subject: [PATCH 18/34] Define grid IDs using IntEnum instead of Enum --- pygmt/datatypes/header.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygmt/datatypes/header.py b/pygmt/datatypes/header.py index db1428036bc..6f3f5386fb3 100644 --- a/pygmt/datatypes/header.py +++ b/pygmt/datatypes/header.py @@ -3,7 +3,7 @@ """ import ctypes as ctp -from enum import Enum +from enum import IntEnum from typing import Any, ClassVar import numpy as np @@ -24,7 +24,7 @@ gmt_grdfloat = ctp.c_float -class GMTGridID(Enum): +class GMTGridID(IntEnum): """ Enum for the GMT grid format ID. From b1eacf145b746f6fbfb944eb5d61ca4513ea6d23 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 15:04:01 +0800 Subject: [PATCH 19/34] Pass the new shape as a positional parameter --- pygmt/datatypes/image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index 4bbdc352968..25c3f2677f3 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -181,7 +181,7 @@ def to_dataarray(self) -> xr.DataArray: pad = header.pad[:] data: np.ndarray = np.reshape( self.data[: header.n_bands * header.mx * header.my], - shape=(header.my, header.mx, header.n_bands), + (header.my, header.mx, header.n_bands), )[pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1], :] # Get x and y coordinates From aa4fdc94d4c2756d6a00f32e4a0c8ec3daa6b207 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 15:17:10 +0800 Subject: [PATCH 20/34] Fix failing tests --- pygmt/datatypes/image.py | 2 +- pygmt/tests/test_clib_read_data.py | 8 +++++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index 25c3f2677f3..22d07779ec5 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -115,7 +115,7 @@ def to_dataarray(self) -> xr.DataArray: ... # Convert to xarray.DataArray and use it later ... da = image.contents.to_dataarray() >>> da # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS - Size: 2MB + ... array([[[ 10, 10, 10, ..., 10, 10, 10], [ 10, 10, 10, ..., 10, 10, 10], [ 10, 10, 10, ..., 10, 10, 10], diff --git a/pygmt/tests/test_clib_read_data.py b/pygmt/tests/test_clib_read_data.py index 5f3ad68fd83..bd91f4774c7 100644 --- a/pygmt/tests/test_clib_read_data.py +++ b/pygmt/tests/test_clib_read_data.py @@ -34,9 +34,11 @@ def fixture_expected_xrimage(): """ The expected xr.DataArray object for the @earth_day_01d_p file. """ - with rioxarray.open_rasterio(which("@earth_day_01d_p")) as da: - dataarray = da.load().drop_vars("spatial_ref") - return dataarray + if _HAS_RIOXARRAY: + with rioxarray.open_rasterio(which("@earth_day_01d_p")) as da: + dataarray = da.load().drop_vars("spatial_ref") + return dataarray + return None def test_clib_read_data_dataset(): From c87a3ec50aa06e8ba2a19fb3691155372bc1edc9 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 15:20:38 +0800 Subject: [PATCH 21/34] One more fix --- pygmt/datatypes/image.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index 22d07779ec5..6eb590f2c8b 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -151,13 +151,13 @@ def to_dataarray(self) -> xr.DataArray: actual_range: [ 1.79769313e+308 -1.79769313e+308] >>> da.coords["x"] # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS - Size: 3kB + ... array([-179.5, -178.5, -177.5, ..., 177.5, 178.5, 179.5]) Coordinates: * x (x) float64 3kB -179.5 -178.5 -177.5 -176.5 ... 177.5 178.5 179.5 >>> da.coords["y"] # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS - Size: 1kB + ... array([ 89.5, 88.5, 87.5, 86.5, 85.5, 84.5, 83.5, 82.5, 81.5, 80.5, 79.5, 78.5, 77.5, 76.5, 75.5, 74.5, 73.5, 72.5, 71.5, 70.5, 69.5, 68.5, 67.5, 66.5, 65.5, 64.5, 63.5, 62.5, 61.5, 60.5, From a20d8a273276f438b1b70071859a8b26b19554bb Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 15:23:43 +0800 Subject: [PATCH 22/34] One more fix --- pygmt/datatypes/image.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index 6eb590f2c8b..8cdc0d13d32 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -154,7 +154,7 @@ def to_dataarray(self) -> xr.DataArray: ... array([-179.5, -178.5, -177.5, ..., 177.5, 178.5, 179.5]) Coordinates: - * x (x) float64 3kB -179.5 -178.5 -177.5 -176.5 ... 177.5 178.5 179.5 + * x (x) float64... -179.5 -178.5 -177.5 -176.5 ... 177.5 178.5 179.5 >>> da.coords["y"] # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS ... @@ -168,7 +168,7 @@ def to_dataarray(self) -> xr.DataArray: -70.5, -71.5, -72.5, -73.5, -74.5, -75.5, -76.5, -77.5, -78.5, -79.5, -80.5, -81.5, -82.5, -83.5, -84.5, -85.5, -86.5, -87.5, -88.5, -89.5]) Coordinates: - * y (y) float64 1kB 89.5 88.5 87.5 86.5 ... -86.5 -87.5 -88.5 -89.5 + * y (y) float64... 89.5 88.5 87.5 86.5 ... -86.5 -87.5 -88.5 -89.5 >>> da.gmt.registration, da.gmt.gtype (1, 0) From 926427bb28f183f624856625b1860adf8f9263de Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 15:37:54 +0800 Subject: [PATCH 23/34] Simplify a doctest --- pygmt/datatypes/image.py | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index 8cdc0d13d32..9c5f2085e1b 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -158,22 +158,13 @@ def to_dataarray(self) -> xr.DataArray: >>> da.coords["y"] # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS ... - array([ 89.5, 88.5, 87.5, 86.5, 85.5, 84.5, 83.5, 82.5, 81.5, 80.5, - 79.5, 78.5, 77.5, 76.5, 75.5, 74.5, 73.5, 72.5, 71.5, 70.5, - 69.5, 68.5, 67.5, 66.5, 65.5, 64.5, 63.5, 62.5, 61.5, 60.5, - ... - -0.5, -1.5, -2.5, -3.5, -4.5, -5.5, -6.5, -7.5, -8.5, -9.5, - ... - -60.5, -61.5, -62.5, -63.5, -64.5, -65.5, -66.5, -67.5, -68.5, -69.5, - -70.5, -71.5, -72.5, -73.5, -74.5, -75.5, -76.5, -77.5, -78.5, -79.5, - -80.5, -81.5, -82.5, -83.5, -84.5, -85.5, -86.5, -87.5, -88.5, -89.5]) + array([ 89.5, 88.5, 87.5, 86.5, ..., -87.5, -88.5, -89.5]) Coordinates: * y (y) float64... 89.5 88.5 87.5 86.5 ... -86.5 -87.5 -88.5 -89.5 >>> da.gmt.registration, da.gmt.gtype (1, 0) """ - # Get image header header: _GMT_GRID_HEADER = self.header.contents From c73328e5669884ed46895fbc59a54fc84dcf0c53 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 16:07:19 +0800 Subject: [PATCH 24/34] Improve the tests --- pygmt/tests/test_clib_read_data.py | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/pygmt/tests/test_clib_read_data.py b/pygmt/tests/test_clib_read_data.py index bd91f4774c7..9311fdc8f62 100644 --- a/pygmt/tests/test_clib_read_data.py +++ b/pygmt/tests/test_clib_read_data.py @@ -4,6 +4,7 @@ from pathlib import Path +import numpy as np import pandas as pd import pytest import xarray as xr @@ -138,6 +139,13 @@ def test_clib_read_data_grid_actual_image(expected_xrimage): ) else: assert xrimage.shape == (180, 360) + assert xrimage.coords["x"].data.min() == -179.5 + assert xrimage.coords["x"].data.max() == 179.5 + assert xrimage.coords["y"].data.min() == -89.5 + assert xrimage.coords["y"].data.max() == 89.5 + assert xrimage.data.min() == 10.0 + assert xrimage.data.max() == 255.0 + assert xrimage.data.dtype == np.float64 def test_clib_read_data_image(expected_xrimage): @@ -149,14 +157,21 @@ def test_clib_read_data_image(expected_xrimage): header = image.header.contents assert header.n_rows == 180 assert header.n_columns == 360 - assert header.n_bands == 3 assert header.wesn[:] == [-180.0, 180.0, -90.0, 90.0] + assert header.n_bands == 3 xrimage = image.to_dataarray() if _HAS_RIOXARRAY: # Full check if rioxarray is installed. xr.testing.assert_equal(xrimage, expected_xrimage) else: assert xrimage.shape == (3, 180, 360) + assert xrimage.coords["x"].data.min() == -179.5 + assert xrimage.coords["x"].data.max() == 179.5 + assert xrimage.coords["y"].data.min() == -89.5 + assert xrimage.coords["y"].data.max() == 89.5 + assert xrimage.data.min() == 10 + assert xrimage.data.max() == 255 + assert xrimage.data.dtype == np.int64 def test_clib_read_data_image_two_steps(expected_xrimage): @@ -184,6 +199,13 @@ def test_clib_read_data_image_two_steps(expected_xrimage): xr.testing.assert_equal(xrimage, expected_xrimage) else: assert xrimage.shape == (3, 180, 360) + assert xrimage.coords["x"].data.min() == -179.5 + assert xrimage.coords["x"].data.max() == 179.5 + assert xrimage.coords["y"].data.min() == -89.5 + assert xrimage.coords["y"].data.max() == 89.5 + assert xrimage.data.min() == 10 + assert xrimage.data.max() == 255 + assert xrimage.data.dtype == np.int64 def test_clib_read_data_fails(): From fb97daa9f032295a94736770d64c480463ca4812 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 19:58:52 +0800 Subject: [PATCH 25/34] Convert ctypes array to numpy array using np.ctypeslib.as_array --- pygmt/datatypes/grid.py | 62 +++++++++++++++++++++------------------- pygmt/datatypes/image.py | 19 ++++++------ 2 files changed, 41 insertions(+), 40 deletions(-) diff --git a/pygmt/datatypes/grid.py b/pygmt/datatypes/grid.py index b420db6d7f1..4be2599afd5 100644 --- a/pygmt/datatypes/grid.py +++ b/pygmt/datatypes/grid.py @@ -40,16 +40,15 @@ class _GMT_GRID(ctp.Structure): # noqa: N801 ... 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]) + ... x = np.ctypeslib.as_array(grid.x, shape=(header.n_columns,)).copy() + ... y = np.ctypeslib.as_array(grid.y, shape=(header.n_rows,)).copy() ... # The data array (with paddings) - ... data = np.reshape( - ... grid.data[: header.mx * header.my], (header.my, header.mx) - ... ) + ... data = np.ctypeslib.as_array( + ... grid.data, shape=(header.my, header.mx) + ... ).copy() ... # 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 @@ -61,22 +60,26 @@ class _GMT_GRID(ctp.Structure): # noqa: N801 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. ]] + >>> x + array([-54.5, -53.5, -52.5, -51.5, -50.5, -49.5, -48.5, -47.5]) + >>> y + array([-10.5, -11.5, -12.5, -13.5, ..., -20.5, -21.5, -22.5, -23.5]) + >>> data + array([[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. ]], + dtype=float32) """ _fields_: ClassVar = [ @@ -126,7 +129,8 @@ def to_dataarray(self) -> xr.DataArray: [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. ]]) + [347.5, 331.5, 309. , 282. , 190. , 208. , 299.5, 348. ]], + dtype=float32) 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 @@ -169,16 +173,14 @@ def to_dataarray(self) -> xr.DataArray: # 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]), - ] + x = np.ctypeslib.as_array(self.x, shape=(header.n_columns,)).copy() + y = np.ctypeslib.as_array(self.y, shape=(header.n_rows,)).copy() + coords = [(dims[0], y, dim_attrs[0]), (dims[1], x, dim_attrs[1])] # The data array without paddings + data = np.ctypeslib.as_array(self.data, shape=(header.my, header.mx)).copy() 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] - ] + data = data[pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1]] # Create the xarray.DataArray object grid = xr.DataArray( diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index 672c640c67d..aaccc75cdf5 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -38,13 +38,12 @@ class _GMT_IMAGE(ctp.Structure): # noqa: N801 ... # Image-specific attributes. ... print(image.type, image.n_indexed_colors) ... # The x and y coordinates - ... x = image.x[: header.n_columns] - ... y = image.y[: header.n_rows] + ... x = np.ctypeslib.as_array(image.x, shape=(header.n_columns,)).copy() + ... y = np.ctypeslib.as_array(image.y, shape=(header.n_rows,)).copy() ... # The data array (with paddings) - ... data = np.reshape( - ... image.data[: header.n_bands * header.mx * header.my], - ... (header.my, header.mx, header.n_bands), - ... ) + ... data = np.ctypeslib.as_array( + ... image.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], :] @@ -60,10 +59,10 @@ class _GMT_IMAGE(ctp.Structure): # noqa: N801 [2, 2, 2, 2] b'BRPa' 0.5 1 0 - >>> x - [-179.5, -178.5, ..., 178.5, 179.5] - >>> y - [89.5, 88.5, ..., -88.5, -89.5] + >>> x # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS + array([-179.5, -178.5, ..., 178.5, 179.5]) + >>> y # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS + array([ 89.5, 88.5, ..., -88.5, -89.5]) >>> data.shape (180, 360, 3) >>> data.min(), data.max() From 15b8d530dc1c8de1e01317bbe303ecc83c6cfeeb Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 22:25:41 +0800 Subject: [PATCH 26/34] Fix the incorrect value due to floating number conversion in sphinterpolate --- pygmt/tests/test_sphinterpolate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygmt/tests/test_sphinterpolate.py b/pygmt/tests/test_sphinterpolate.py index 5e5485dc2cc..d7e24eba780 100644 --- a/pygmt/tests/test_sphinterpolate.py +++ b/pygmt/tests/test_sphinterpolate.py @@ -41,4 +41,4 @@ def test_sphinterpolate_no_outgrid(mars): npt.assert_allclose(temp_grid.max(), 14628.144) npt.assert_allclose(temp_grid.min(), -6908.1987) npt.assert_allclose(temp_grid.median(), 118.96849) - npt.assert_allclose(temp_grid.mean(), 272.60578) + npt.assert_allclose(temp_grid.mean(), 272.60593) From 3e3a6f3a56c21272898754296be8ca1427f038cb Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 23:57:48 +0800 Subject: [PATCH 27/34] Update the to_dataarray method to match the codes in GMT_GRID --- pygmt/datatypes/image.py | 51 ++++++++++++++++++++++++---------------- 1 file changed, 31 insertions(+), 20 deletions(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index f03751414cd..545ed2ea349 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -65,6 +65,8 @@ class _GMT_IMAGE(ctp.Structure): # noqa: N801 array([-179.5, -178.5, ..., 178.5, 179.5]) >>> y # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS array([ 89.5, 88.5, ..., -88.5, -89.5]) + >>> data.dtype + dtype('uint8') >>> data.shape (180, 360, 3) >>> data.min(), data.max() @@ -108,7 +110,7 @@ def to_dataarray(self) -> xr.DataArray: >>> from pygmt.clib import Session >>> with Session() as lib: ... with lib.virtualfile_out(kind="image") as voutimg: - ... lib.call_module("read", ["@earth_day_01d", voutimg, "-Ti"]) + ... lib.call_module("read", ["@earth_day_01d_p", voutimg, "-Ti"]) ... # Read the image from the virtual file ... image = lib.read_virtualfile(voutimg, kind="image") ... # Convert to xarray.DataArray and use it later @@ -137,10 +139,10 @@ def to_dataarray(self) -> xr.DataArray: ..., [177, 179, 179, ..., 178, 177, 177], [185, 187, 187, ..., 187, 186, 185], - [189, 191, 191, ..., 191, 191, 189]]]) + [189, 191, 191, ..., 191, 191, 189]]], dtype=uint8) Coordinates: - * x (x) float64... -179.5 -178.5 -177.5 -176.5 ... 177.5 178.5 179.5 * y (y) float64... 89.5 88.5 87.5 86.5 ... -86.5 -87.5 -88.5 -89.5 + * x (x) float64... -179.5 -178.5 -177.5 -176.5 ... 177.5 178.5 179.5 * band (band) uint8... 1 2 3 Attributes: title: @@ -154,41 +156,50 @@ def to_dataarray(self) -> xr.DataArray: array([-179.5, -178.5, -177.5, ..., 177.5, 178.5, 179.5]) Coordinates: * x (x) float64... -179.5 -178.5 -177.5 -176.5 ... 177.5 178.5 179.5 - + Attributes: + long_name: x + axis: X + actual_range: [-180. 180.] >>> da.coords["y"] # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS ... array([ 89.5, 88.5, 87.5, 86.5, ..., -87.5, -88.5, -89.5]) Coordinates: * y (y) float64... 89.5 88.5 87.5 86.5 ... -86.5 -87.5 -88.5 -89.5 - + Attributes: + long_name: y + axis: Y + actual_range: [-90. 90.] >>> da.gmt.registration, da.gmt.gtype (1, 0) """ - # Get image header - header: _GMT_GRID_HEADER = self.header.contents + # The image 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) + x = np.ctypeslib.as_array(self.x, shape=(header.n_columns,)).copy() + y = np.ctypeslib.as_array(self.y, shape=(header.n_rows,)).copy() + coords = [ + (dims[0], y, dim_attrs[0]), + (dims[1], x, dim_attrs[1]), + ("band", np.array([1, 2, 3], dtype=np.uint8), None), + ] # Get DataArray without padding + data = np.ctypeslib.as_array( + self.data, shape=(header.my, header.mx, header.n_bands) + ).copy() pad = header.pad[:] - data: np.ndarray = np.reshape( - self.data[: header.n_bands * header.mx * header.my], - (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([1, 2, 3], dtype=np.uint8), - } + data = data[pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1], :] # Create the xarray.DataArray object image = xr.DataArray( data=data, coords=coords, - dims=("y", "x", "band"), name=header.name, attrs=header.data_attrs, - ).transpose("band", "y", "x") + ).transpose("band", dims[0], dims[1]) # Set GMT accessors. # Must put at the end, otherwise info gets lost after certain image operations. From 12ef40ab566ec325847b63e2648a8ba1318ac546 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sat, 21 Sep 2024 00:31:15 +0800 Subject: [PATCH 28/34] image data should has uint8 dtype --- pygmt/tests/test_clib_read_data.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygmt/tests/test_clib_read_data.py b/pygmt/tests/test_clib_read_data.py index 9311fdc8f62..bf89a1d8086 100644 --- a/pygmt/tests/test_clib_read_data.py +++ b/pygmt/tests/test_clib_read_data.py @@ -171,7 +171,7 @@ def test_clib_read_data_image(expected_xrimage): assert xrimage.coords["y"].data.max() == 89.5 assert xrimage.data.min() == 10 assert xrimage.data.max() == 255 - assert xrimage.data.dtype == np.int64 + assert xrimage.data.dtype == np.uint8 def test_clib_read_data_image_two_steps(expected_xrimage): From f64fbb84825fb068a6c01194f81fa885709a4277 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sat, 21 Sep 2024 11:45:46 +0800 Subject: [PATCH 29/34] Further improve the tests --- pygmt/datatypes/header.py | 2 +- pygmt/tests/test_clib_read_data.py | 67 +++++++++++++----------------- 2 files changed, 31 insertions(+), 38 deletions(-) diff --git a/pygmt/datatypes/header.py b/pygmt/datatypes/header.py index 6f3f5386fb3..448367b7edb 100644 --- a/pygmt/datatypes/header.py +++ b/pygmt/datatypes/header.py @@ -244,7 +244,7 @@ def data_attrs(self) -> dict[str, Any]: GMTGridID.GMT_GRID_IS_NI, GMTGridID.GMT_GRID_IS_NF, GMTGridID.GMT_GRID_IS_ND, - }: + }: # netCDF format attrs["Conventions"] = "CF-1.7" attrs["title"] = self.title.decode() attrs["history"] = self.command.decode() diff --git a/pygmt/tests/test_clib_read_data.py b/pygmt/tests/test_clib_read_data.py index bf89a1d8086..fd34fbced9f 100644 --- a/pygmt/tests/test_clib_read_data.py +++ b/pygmt/tests/test_clib_read_data.py @@ -123,29 +123,27 @@ def test_clib_read_data_grid_actual_image(expected_xrimage): """ with Session() as lib: image = lib.read_data("@earth_day_01d_p", kind="grid").contents - header = image.header.contents - assert header.n_rows == 180 - assert header.n_columns == 360 - assert header.wesn[:] == [-180.0, 180.0, -90.0, 90.0] # Explicitly check n_bands. Only one band is read for 3-band images. - assert header.n_bands == 1 + assert image.header.contents.n_bands == 1 xrimage = image.to_dataarray() + assert xrimage.shape == (180, 360) + assert xrimage.coords["x"].data.min() == -179.5 + assert xrimage.coords["x"].data.max() == 179.5 + assert xrimage.coords["y"].data.min() == -89.5 + assert xrimage.coords["y"].data.max() == 89.5 + assert xrimage.data.min() == 10.0 + assert xrimage.data.max() == 255.0 + # Data are stored as uint8 in images but are converted to float32 when reading + # into a GMT_GRID container. + assert xrimage.data.dtype == np.float32 + if _HAS_RIOXARRAY: # Full check if rioxarray is installed. assert expected_xrimage.band.size == 3 # 3-band image. xr.testing.assert_equal( xrimage, expected_xrimage.isel(band=0).drop_vars(["band"]).sortby("y"), ) - else: - assert xrimage.shape == (180, 360) - assert xrimage.coords["x"].data.min() == -179.5 - assert xrimage.coords["x"].data.max() == 179.5 - assert xrimage.coords["y"].data.min() == -89.5 - assert xrimage.coords["y"].data.max() == 89.5 - assert xrimage.data.min() == 10.0 - assert xrimage.data.max() == 255.0 - assert xrimage.data.dtype == np.float64 def test_clib_read_data_image(expected_xrimage): @@ -154,24 +152,19 @@ def test_clib_read_data_image(expected_xrimage): """ with Session() as lib: image = lib.read_data("@earth_day_01d_p", kind="image").contents - header = image.header.contents - assert header.n_rows == 180 - assert header.n_columns == 360 - assert header.wesn[:] == [-180.0, 180.0, -90.0, 90.0] - assert header.n_bands == 3 xrimage = image.to_dataarray() + assert xrimage.shape == (3, 180, 360) + assert xrimage.coords["x"].data.min() == -179.5 + assert xrimage.coords["x"].data.max() == 179.5 + assert xrimage.coords["y"].data.min() == -89.5 + assert xrimage.coords["y"].data.max() == 89.5 + assert xrimage.data.min() == 10 + assert xrimage.data.max() == 255 + assert xrimage.data.dtype == np.uint8 + if _HAS_RIOXARRAY: # Full check if rioxarray is installed. xr.testing.assert_equal(xrimage, expected_xrimage) - else: - assert xrimage.shape == (3, 180, 360) - assert xrimage.coords["x"].data.min() == -179.5 - assert xrimage.coords["x"].data.max() == 179.5 - assert xrimage.coords["y"].data.min() == -89.5 - assert xrimage.coords["y"].data.max() == 89.5 - assert xrimage.data.min() == 10 - assert xrimage.data.max() == 255 - assert xrimage.data.dtype == np.uint8 def test_clib_read_data_image_two_steps(expected_xrimage): @@ -195,17 +188,17 @@ def test_clib_read_data_image_two_steps(expected_xrimage): lib.read_data(infile, kind="image", mode="GMT_DATA_ONLY", data=data_ptr) xrimage = image.to_dataarray() + assert xrimage.shape == (3, 180, 360) + assert xrimage.coords["x"].data.min() == -179.5 + assert xrimage.coords["x"].data.max() == 179.5 + assert xrimage.coords["y"].data.min() == -89.5 + assert xrimage.coords["y"].data.max() == 89.5 + assert xrimage.data.min() == 10 + assert xrimage.data.max() == 255 + assert xrimage.data.dtype == np.uint8 + if _HAS_RIOXARRAY: # Full check if rioxarray is installed. xr.testing.assert_equal(xrimage, expected_xrimage) - else: - assert xrimage.shape == (3, 180, 360) - assert xrimage.coords["x"].data.min() == -179.5 - assert xrimage.coords["x"].data.max() == 179.5 - assert xrimage.coords["y"].data.min() == -89.5 - assert xrimage.coords["y"].data.max() == 89.5 - assert xrimage.data.min() == 10 - assert xrimage.data.max() == 255 - assert xrimage.data.dtype == np.int64 def test_clib_read_data_fails(): From d49afed50d4fc52992c640436c882f78870aa099 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Tue, 24 Sep 2024 10:16:53 +0800 Subject: [PATCH 30/34] Add a note that currently only 3-band images are supported --- pygmt/datatypes/image.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index 545ed2ea349..86ea7d691bd 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -175,6 +175,14 @@ def to_dataarray(self) -> xr.DataArray: # The image header header = self.header.contents + if header.n_bands != 3: + msg = ( + f"The raster image has {header.n_bands}-band(s). " + "Currently only 3-band images are supported. " + "Please consider submitting a feature request to us. " + ) + raise NotImplementedError(msg) + # 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) From 2fd13fb7b8efd9b8438440b84199244a71397e05 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sat, 28 Sep 2024 21:25:44 +0800 Subject: [PATCH 31/34] Remove the old GMTGridID enums from pygmt/datatypes/header.py --- pygmt/datatypes/header.py | 35 ----------------------------------- 1 file changed, 35 deletions(-) diff --git a/pygmt/datatypes/header.py b/pygmt/datatypes/header.py index a097df19eeb..d7806e50a24 100644 --- a/pygmt/datatypes/header.py +++ b/pygmt/datatypes/header.py @@ -3,7 +3,6 @@ """ import ctypes as ctp -from enum import IntEnum from typing import Any, ClassVar import numpy as np @@ -25,40 +24,6 @@ gmt_grdfloat = ctp.c_float -class GMTGridID(IntEnum): - """ - Enum for the GMT grid format ID. - - Defined in gmt_grdio.h. - """ - - GMT_GRD_UNKNOWN_FMT = 0 # if grid format cannot be auto-detected - GMT_GRID_IS_BF = 1 # GMT native, C-binary format (32-bit float) - GMT_GRID_IS_BS = 2 # GMT native, C-binary format (16-bit integer) - GMT_GRID_IS_RB = 3 # SUN rasterfile format (8-bit standard) - GMT_GRID_IS_BB = 4 # GMT native, C-binary format (8-bit integer) - GMT_GRID_IS_BM = 5 # GMT native, C-binary format (bit-mask) - GMT_GRID_IS_SF = 6 # Golden Software Surfer format 6 (32-bit float) - GMT_GRID_IS_CB = 7 # GMT netCDF format (8-bit integer) - GMT_GRID_IS_CS = 8 # GMT netCDF format (16-bit integer) - GMT_GRID_IS_CI = 9 # GMT netCDF format (32-bit integer) - GMT_GRID_IS_CF = 10 # GMT netCDF format (32-bit float) - GMT_GRID_IS_CD = 11 # GMT netCDF format (64-bit float) - GMT_GRID_IS_RF = 12 # GEODAS grid format GRD98 (NGDC) - GMT_GRID_IS_BI = 13 # GMT native, C-binary format (32-bit integer) - GMT_GRID_IS_BD = 14 # GMT native, C-binary format (64-bit float) - GMT_GRID_IS_NB = 15 # GMT netCDF format (8-bit integer) - GMT_GRID_IS_NS = 16 # GMT netCDF format (16-bit integer) - GMT_GRID_IS_NI = 17 # GMT netCDF format (32-bit integer) - GMT_GRID_IS_NF = 18 # GMT netCDF format (32-bit float) - GMT_GRID_IS_ND = 19 # GMT netCDF format (64-bit float) - GMT_GRID_IS_SD = 20 # Golden Software Surfer format 7 (64-bit float, read-only) - GMT_GRID_IS_AF = 21 # Atlantic Geoscience Center format AGC (32-bit float) - GMT_GRID_IS_GD = 22 # Import through GDAL - GMT_GRID_IS_EI = 23 # ESRI Arc/Info ASCII Grid Interchange format (ASCII integer) - GMT_GRID_IS_EF = 24 # ESRI Arc/Info ASCII Grid Interchange format (ASCII float) - - def _parse_nameunits(nameunits: str) -> tuple[str, str | None]: """ Get the long_name and units attributes from x_units/y_units/z_units in the grid From 9972ba17446fd780b18a74008aff545b520815e0 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sat, 28 Sep 2024 22:34:37 +0800 Subject: [PATCH 32/34] A minor fix --- pygmt/datatypes/image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index 86ea7d691bd..037db5c2a6e 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -177,7 +177,7 @@ def to_dataarray(self) -> xr.DataArray: if header.n_bands != 3: msg = ( - f"The raster image has {header.n_bands}-band(s). " + f"The raster image has {header.n_bands} band(s). " "Currently only 3-band images are supported. " "Please consider submitting a feature request to us. " ) From 62f0ce056cb6bdd32521905529246d67faae2b8d Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 29 Sep 2024 12:26:46 +0800 Subject: [PATCH 33/34] Set CF-1.7-specific attributes for netCDF formats only --- pygmt/datatypes/grid.py | 2 +- pygmt/datatypes/header.py | 10 +++++----- pygmt/datatypes/image.py | 4 ---- 3 files changed, 6 insertions(+), 10 deletions(-) diff --git a/pygmt/datatypes/grid.py b/pygmt/datatypes/grid.py index 4be2599afd5..18ad679aeef 100644 --- a/pygmt/datatypes/grid.py +++ b/pygmt/datatypes/grid.py @@ -139,8 +139,8 @@ def to_dataarray(self) -> xr.DataArray: 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.] + long_name: elevation (m) >>> da.coords["lon"] # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS ... array([-54.5, -53.5, -52.5, -51.5, -50.5, -49.5, -48.5, -47.5]) diff --git a/pygmt/datatypes/header.py b/pygmt/datatypes/header.py index d7806e50a24..1cef6d38b38 100644 --- a/pygmt/datatypes/header.py +++ b/pygmt/datatypes/header.py @@ -210,17 +210,17 @@ def data_attrs(self) -> dict[str, Any]: GridFormat.NI, GridFormat.NF, GridFormat.ND, - }: # Only set the 'Conventions' attribute for netCDF. + }: # Set attributes specific to CF-1.7 conventions attrs["Conventions"] = "CF-1.7" - attrs["title"] = self.title.decode() - attrs["history"] = self.command.decode() - attrs["description"] = self.remark.decode() + attrs["title"] = self.title.decode() + attrs["history"] = self.command.decode() + attrs["description"] = self.remark.decode() + attrs["actual_range"] = np.array([self.z_min, self.z_max]) long_name, units = _parse_nameunits(self.z_units.decode()) if long_name: attrs["long_name"] = long_name if units: attrs["units"] = units - attrs["actual_range"] = np.array([self.z_min, self.z_max]) return attrs @property diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index 037db5c2a6e..773c3cdb78a 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -145,11 +145,7 @@ def to_dataarray(self) -> xr.DataArray: * x (x) float64... -179.5 -178.5 -177.5 -176.5 ... 177.5 178.5 179.5 * band (band) uint8... 1 2 3 Attributes: - title: - history: - description: long_name: z - actual_range: [ 1.79769313e+308 -1.79769313e+308] >>> da.coords["x"] # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS ... From f715aee1d1a5497d4018c2db32d83d77661f88d8 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 29 Sep 2024 15:27:58 +0800 Subject: [PATCH 34/34] Improve the logic for determine grid/image gtype based on upstream GMT codes --- pygmt/datatypes/header.py | 12 +++++++++++- pygmt/datatypes/image.py | 2 +- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/pygmt/datatypes/header.py b/pygmt/datatypes/header.py index 1cef6d38b38..197eaa26e45 100644 --- a/pygmt/datatypes/header.py +++ b/pygmt/datatypes/header.py @@ -250,6 +250,16 @@ def gtype(self) -> int: "lon"/"lat" or have units "degrees_east"/"degrees_north", then the grid is assumed to be geographic. """ + gtype = 0 # Cartesian by default + dims = self.dims - gtype = 1 if dims[0] == "lat" and dims[1] == "lon" else 0 + if dims[0] == "lat" and dims[1] == "lon": + # Check dimensions for grids that following CF-conventions + gtype = 1 + elif self.ProjRefPROJ4 is not None: + # Check ProjRefPROJ4 for images imported via GDAL. + # The logic comes from GMT's `gmtlib_read_image_info` function. + projref = self.ProjRefPROJ4.decode() + if "longlat" in projref or "latlong" in projref: + gtype = 1 return gtype diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index 773c3cdb78a..743cc003681 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -166,7 +166,7 @@ def to_dataarray(self) -> xr.DataArray: axis: Y actual_range: [-90. 90.] >>> da.gmt.registration, da.gmt.gtype - (1, 0) + (1, 1) """ # The image header header = self.header.contents