From bcf43f0f9e8bc45137129bd48662386180d4d7b6 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 18 Mar 2024 22:34:50 +0800 Subject: [PATCH 01/73] 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/73] 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 59d523c2f5cccd30448e6fc1288cb11b39fa4e76 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Tue, 16 Apr 2024 15:05:12 +0800 Subject: [PATCH 03/73] pygmt.grdcut: Support both grid and image output --- pygmt/clib/session.py | 24 +++++++++++++++++++----- pygmt/datatypes/__init__.py | 1 + pygmt/datatypes/image.py | 9 +++++++++ pygmt/src/grdcut.py | 36 +++++++++++++++++++++++------------- 4 files changed, 52 insertions(+), 18 deletions(-) create mode 100644 pygmt/datatypes/image.py diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 8a8b52df8e5..cd4637093e0 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -25,7 +25,7 @@ vectors_to_arrays, ) from pygmt.clib.loading import load_libgmt -from pygmt.datatypes import _GMT_DATASET, _GMT_GRID +from pygmt.datatypes import _GMT_DATASET, _GMT_GRID, _GMT_IMAGE from pygmt.exceptions import ( GMTCLibError, GMTCLibNoSessionError, @@ -33,6 +33,7 @@ GMTVersionError, ) from pygmt.helpers import ( + GMTTempFile, data_kind, fmt_docstring, tempfile_from_geojson, @@ -1649,7 +1650,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. @@ -1706,8 +1709,11 @@ def virtualfile_out( family, geometry = { "dataset": ("GMT_IS_DATASET", "GMT_IS_PLP"), "grid": ("GMT_IS_GRID", "GMT_IS_SURFACE"), + "image": ("GMT_IS_IMAGE", "GMT_IS_SURFACE"), }[kind] - with self.open_virtualfile(family, geometry, "GMT_OUT", None) as vfile: + with self.open_virtualfile( + family, geometry, "GMT_OUT|GMT_IS_REFERENCE", None + ) as vfile: yield vfile def inquire_virtualfile(self, vfname: str) -> int: @@ -1801,9 +1807,13 @@ def read_virtualfile( # _GMT_DATASET). if kind is None: # Return the ctypes void pointer return pointer - if kind in ["image", "cube"]: + if kind == "cube": raise NotImplementedError(f"kind={kind} is not supported yet.") - dtype = {"dataset": _GMT_DATASET, "grid": _GMT_GRID}[kind] + dtype = { + "dataset": _GMT_DATASET, + "grid": _GMT_GRID, + "image": _GMT_IMAGE, + }[kind] return ctp.cast(pointer, ctp.POINTER(dtype)) def virtualfile_to_dataset( @@ -2013,6 +2023,10 @@ def virtualfile_to_raster( self["GMT_IS_IMAGE"]: "image", self["GMT_IS_CUBE"]: "cube", }[family] + if kind == "image": + with GMTTempFile(suffix=".tif") as tmpfile: + self.call_module("write", f"{vfname} {tmpfile.name} -Ti") + return xr.load_dataarray(tmpfile.name) return self.read_virtualfile(vfname, kind=kind).contents.to_dataarray() def extract_region(self): 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..f67aefc1fe2 --- /dev/null +++ b/pygmt/datatypes/image.py @@ -0,0 +1,9 @@ +""" +Wrapper for the GMT_GRID data type. +""" + +import ctypes as ctp + + +class _GMT_IMAGE(ctp.Structure): # noqa: N801 + pass diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index 9ffcda213d6..279e4cee838 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -4,20 +4,19 @@ from pygmt.clib import Session from pygmt.helpers import ( - GMTTempFile, build_arg_string, + data_kind, fmt_docstring, kwargs_to_strings, use_alias, ) -from pygmt.io import load_dataarray +from pygmt.src.which import which __doctest_skip__ = ["grdcut"] @fmt_docstring @use_alias( - G="outgrid", R="region", J="projection", N="extend", @@ -27,7 +26,7 @@ f="coltypes", ) @kwargs_to_strings(R="sequence") -def grdcut(grid, **kwargs): +def grdcut(grid, outgrid: str | None = None, **kwargs): r""" Extract subregion from a grid. @@ -99,13 +98,24 @@ def grdcut(grid, **kwargs): >>> # 12° E to 15° E and a latitude range of 21° N to 24° N >>> new_grid = pygmt.grdcut(grid=grid, region=[12, 15, 21, 24]) """ - with GMTTempFile(suffix=".nc") as tmpfile: - with Session() as lib: - with lib.virtualfile_in(check_kind="raster", data=grid) as vingrd: - if (outgrid := kwargs.get("G")) is None: - kwargs["G"] = outgrid = tmpfile.name # output to tmpfile - lib.call_module( - module="grdcut", args=build_arg_string(kwargs, infile=vingrd) - ) + inkind = data_kind(grid) + outkind = "image" if inkind == "image" else "grid" + if inkind == "file": + realpath = which(str(grid)) + if isinstance(realpath, list): + realpath = realpath[0] + if realpath.endswith(".tif"): + outkind = "image" - return load_dataarray(outgrid) if outgrid == tmpfile.name else None + with Session() as lib: + with ( + lib.virtualfile_in(check_kind="raster", data=grid) as vingrd, + lib.virtualfile_out(kind=outkind, fname=outgrid) as voutgrd, + ): + kwargs["G"] = voutgrd + lib.call_module( + module="grdcut", args=build_arg_string(kwargs, infile=vingrd) + ) + return lib.virtualfile_to_raster( + outgrid=outgrid, kind=outkind, vfname=voutgrd + ) From cea337472399ed17d8af7d472832f8f793b9f2b1 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 19 Apr 2024 10:33:29 +0800 Subject: [PATCH 04/73] Fix --- pygmt/src/grdcut.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index 3ff83849a93..49e71c6e548 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -113,9 +113,7 @@ def grdcut(grid, outgrid: str | None = None, **kwargs): lib.virtualfile_out(kind=outkind, fname=outgrid) as voutgrd, ): kwargs["G"] = voutgrd - lib.call_module( - module="grdcut", args=build_arg_string(kwargs, infile=vingrd) - ) + lib.call_module(module="grdcut", args=build_arg_list(kwargs, infile=vingrd)) return lib.virtualfile_to_raster( outgrid=outgrid, kind=outkind, vfname=voutgrd ) From 80d9837dea710a48048ad62977a8f53220798dfc Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sat, 20 Apr 2024 01:00:06 +0800 Subject: [PATCH 05/73] Refactor --- pygmt/src/grdcut.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index 49e71c6e548..3656397b3b6 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -99,13 +99,15 @@ def grdcut(grid, outgrid: str | None = None, **kwargs): >>> new_grid = pygmt.grdcut(grid=grid, region=[12, 15, 21, 24]) """ inkind = data_kind(grid) - outkind = "image" if inkind == "image" else "grid" - if inkind == "file": - realpath = which(str(grid)) - if isinstance(realpath, list): - realpath = realpath[0] - if realpath.endswith(".tif"): - outkind = "image" + match inkind: + case "image" | "grid": + outkind = inkind + case "file": + realpath = which(str(grid)) + if isinstance(realpath, list): + realpath = realpath[0] + if realpath.endswith(".tif"): + outkind = "image" with Session() as lib: with ( From 22fba56ca34961d3db78921ced7e237a87c46305 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sat, 20 Apr 2024 01:01:59 +0800 Subject: [PATCH 06/73] fix --- pygmt/src/grdcut.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index 3656397b3b6..6125629a522 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -103,7 +103,7 @@ def grdcut(grid, outgrid: str | None = None, **kwargs): case "image" | "grid": outkind = inkind case "file": - realpath = which(str(grid)) + realpath = which(grid, download="a") if isinstance(realpath, list): realpath = realpath[0] if realpath.endswith(".tif"): 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 07/73] 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 08/73] 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 09/73] 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 10/73] 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 11/73] 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 12/73] 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 13/73] 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 14/73] 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 3a24ebd3cc4659d005cdc1c326ff97af6816c73f Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Thu, 20 Jun 2024 21:52:23 +0800 Subject: [PATCH 15/73] Apply suggestions from code review Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com> --- pygmt/datatypes/image.py | 2 +- pygmt/src/grdcut.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pygmt/datatypes/image.py b/pygmt/datatypes/image.py index f67aefc1fe2..c4ea0c6d3b3 100644 --- a/pygmt/datatypes/image.py +++ b/pygmt/datatypes/image.py @@ -1,5 +1,5 @@ """ -Wrapper for the GMT_GRID data type. +Wrapper for the GMT_IMAGE data type. """ import ctypes as ctp diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index 6125629a522..577f92f5f04 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -28,7 +28,7 @@ @kwargs_to_strings(R="sequence") def grdcut(grid, outgrid: str | None = None, **kwargs): r""" - Extract subregion from a grid. + Extract subregion from a grid or image. Produce a new ``outgrid`` file which is a subregion of ``grid``. The subregion is specified with ``region``; the specified range must not exceed From 5e390d423cb25bbdccef50f94ab110259c0c2cc4 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Thu, 20 Jun 2024 21:55:25 +0800 Subject: [PATCH 16/73] Address reviewer's comments --- pygmt/clib/session.py | 4 ++-- pygmt/src/grdcut.py | 3 +-- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index e2fc0cf5850..dbaa0e607fa 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -1713,8 +1713,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. diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index 577f92f5f04..9a63936137b 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -106,8 +106,7 @@ def grdcut(grid, outgrid: str | None = None, **kwargs): realpath = which(grid, download="a") if isinstance(realpath, list): realpath = realpath[0] - if realpath.endswith(".tif"): - outkind = "image" + outkind = "image" if realpath.endswith(".tif") else "grid" with Session() as lib: with ( From 003383da831be4c82c89dc4187266d0abf0555ec Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 21 Jun 2024 14:36:27 +0800 Subject: [PATCH 17/73] Fix GMT_OUT --- pygmt/clib/session.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index dbaa0e607fa..062d8ca337d 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -1759,9 +1759,8 @@ def virtualfile_out( "grid": ("GMT_IS_GRID", "GMT_IS_SURFACE"), "image": ("GMT_IS_IMAGE", "GMT_IS_SURFACE"), }[kind] - with self.open_virtualfile( - family, geometry, "GMT_OUT|GMT_IS_REFERENCE", None - ) as vfile: + direction = "GMT_OUT|GMT_IS_REFERENCE" if kind == "image" else "GMT_OUT" + with self.open_virtualfile(family, geometry, direction, None) as vfile: yield vfile def inquire_virtualfile(self, vfname: str) -> int: From 377941a2a1f00d943e6aea78c10b5643643708e7 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 7 Jul 2024 16:59:11 +0800 Subject: [PATCH 18/73] Revert changes for _GMT_IMAGE --- pygmt/clib/session.py | 10 +++------- pygmt/datatypes/__init__.py | 1 - pygmt/datatypes/image.py | 9 --------- 3 files changed, 3 insertions(+), 17 deletions(-) delete mode 100644 pygmt/datatypes/image.py diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 062d8ca337d..e965516741b 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -25,7 +25,7 @@ vectors_to_arrays, ) from pygmt.clib.loading import load_libgmt -from pygmt.datatypes import _GMT_DATASET, _GMT_GRID, _GMT_IMAGE +from pygmt.datatypes import _GMT_DATASET, _GMT_GRID from pygmt.exceptions import ( GMTCLibError, GMTCLibNoSessionError, @@ -1854,13 +1854,9 @@ def read_virtualfile( # _GMT_DATASET). if kind is None: # Return the ctypes void pointer return pointer - if kind == "cube": + if kind in ["image", "cube"]: raise NotImplementedError(f"kind={kind} is not supported yet.") - dtype = { - "dataset": _GMT_DATASET, - "grid": _GMT_GRID, - "image": _GMT_IMAGE, - }[kind] + dtype = {"dataset": _GMT_DATASET, "grid": _GMT_GRID}[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 3489dd19d10..237a050a9f7 100644 --- a/pygmt/datatypes/__init__.py +++ b/pygmt/datatypes/__init__.py @@ -4,4 +4,3 @@ 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 deleted file mode 100644 index c4ea0c6d3b3..00000000000 --- a/pygmt/datatypes/image.py +++ /dev/null @@ -1,9 +0,0 @@ -""" -Wrapper for the GMT_IMAGE data type. -""" - -import ctypes as ctp - - -class _GMT_IMAGE(ctp.Structure): # noqa: N801 - pass From 20617f54c209b07650e3df60ddb26c8b93d8a35f Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 7 Jul 2024 17:03:04 +0800 Subject: [PATCH 19/73] Use rioxarray.open_rasterio for loading images --- pygmt/clib/session.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index e965516741b..e0a895174b1 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -1759,6 +1759,7 @@ def virtualfile_out( "grid": ("GMT_IS_GRID", "GMT_IS_SURFACE"), "image": ("GMT_IS_IMAGE", "GMT_IS_SURFACE"), }[kind] + # For unkown reasons, 'GMT_OUT' crashes for 'image' kind. direction = "GMT_OUT|GMT_IS_REFERENCE" if kind == "image" else "GMT_OUT" with self.open_virtualfile(family, geometry, direction, None) as vfile: yield vfile @@ -2071,10 +2072,14 @@ def virtualfile_to_raster( self["GMT_IS_IMAGE"]: "image", self["GMT_IS_CUBE"]: "cube", }[family] - if kind == "image": + + if kind == "image": # Use temporary file for images + import rioxarray + with GMTTempFile(suffix=".tif") as tmpfile: self.call_module("write", f"{vfname} {tmpfile.name} -Ti") - return xr.load_dataarray(tmpfile.name) + with rioxarray.open_rasterio(tmpfile.name) as dataarray: + return dataarray.load() return self.read_virtualfile(vfname, kind=kind).contents.to_dataarray() def extract_region(self): From a99871862ef66b14dee2d620046ece245c728a81 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 7 Jul 2024 17:05:22 +0800 Subject: [PATCH 20/73] Check if rioxarray is installed --- pygmt/clib/session.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index e0a895174b1..0bc111d3aec 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -39,6 +39,13 @@ tempfile_from_image, ) +try: + import rioxarray + + _HAS_RIOXARRAY = True +except ImportError: + _HAS_RIOXARRAY = False + FAMILIES = [ "GMT_IS_DATASET", # Entity is a data table "GMT_IS_GRID", # Entity is a grid @@ -2074,8 +2081,12 @@ def virtualfile_to_raster( }[family] if kind == "image": # Use temporary file for images - import rioxarray - + if not _HAS_RIOXARRAY: + raise ImportError( + "Package `rioxarray` is required to be installed to load images. " + "Please use `python -m pip install rioxarray` or " + "`mamba install -c conda-forge rioxarray` to install the package." + ) with GMTTempFile(suffix=".tif") as tmpfile: self.call_module("write", f"{vfname} {tmpfile.name} -Ti") with rioxarray.open_rasterio(tmpfile.name) as dataarray: From 86cab444ae6b32c0b7e5142b7a73e4430fa4a6fc Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 7 Jul 2024 17:27:44 +0800 Subject: [PATCH 21/73] Improve grdcut --- pygmt/src/grdcut.py | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index 9a63936137b..e9356308a3c 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -98,15 +98,24 @@ def grdcut(grid, outgrid: str | None = None, **kwargs): >>> # 12° E to 15° E and a latitude range of 21° N to 24° N >>> new_grid = pygmt.grdcut(grid=grid, region=[12, 15, 21, 24]) """ + # Determine the output data kind based on the input data kind. inkind = data_kind(grid) match inkind: case "image" | "grid": outkind = inkind case "file": - realpath = which(grid, download="a") - if isinstance(realpath, list): - realpath = realpath[0] - outkind = "image" if realpath.endswith(".tif") else "grid" + realpath = str(grid) + if realpath.startswith("@"): # Is a remote file without suffix + realpath = which(grid, download="a") + if isinstance(realpath, list): + realpath = realpath[0] + + if realpath.endswith( + ".jpg", ".jpeg", ".png", ".tif", ".tiff", ".bmp", ".webp" + ): + outkind = "image" + else: # Fall back to grid. + outkind = "grid" with Session() as lib: with ( From 6031bab31ba0c013a4713d10d3c73c3254441718 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 7 Jul 2024 18:33:19 +0800 Subject: [PATCH 22/73] Fix typos in grdcut --- pygmt/src/grdcut.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index e9356308a3c..41e5ac5596d 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -111,7 +111,7 @@ def grdcut(grid, outgrid: str | None = None, **kwargs): realpath = realpath[0] if realpath.endswith( - ".jpg", ".jpeg", ".png", ".tif", ".tiff", ".bmp", ".webp" + (".jpg", ".jpeg", ".png", ".tif", ".tiff", ".bmp", ".webp") ): outkind = "image" else: # Fall back to grid. From eb0af2dd9d8ecf0fc5ee8bdcb1a2cd775070d463 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 7 Jul 2024 18:33:32 +0800 Subject: [PATCH 23/73] Add tests for grdcut images --- pygmt/tests/test_grdcut.py | 63 +++++++++++++++++++++++++++++++++++++- 1 file changed, 62 insertions(+), 1 deletion(-) diff --git a/pygmt/tests/test_grdcut.py b/pygmt/tests/test_grdcut.py index dbf5dd21f49..0241679e5da 100644 --- a/pygmt/tests/test_grdcut.py +++ b/pygmt/tests/test_grdcut.py @@ -5,11 +5,18 @@ import numpy as np import pytest import xarray as xr -from pygmt import grdcut, load_dataarray +from pygmt import grdcut, load_dataarray, which from pygmt.exceptions import GMTInvalidInput from pygmt.helpers import GMTTempFile from pygmt.helpers.testing import load_static_earth_relief +try: + import rioxarray + + _HAS_RIOXARRAY = True +except ImportError: + _HAS_RIOXARRAY = False + @pytest.fixture(scope="module", name="grid") def fixture_grid(): @@ -43,6 +50,41 @@ def fixture_expected_grid(): ) +@pytest.fixture(scope="module", name="expected_image") +def fixture_expected_image(): + """ + Load the expected grdcut image result. + """ + return xr.DataArray( + data=np.array( + [ + [[90, 93, 95, 90], [91, 90, 91, 91], [91, 90, 89, 90]], + [[87, 88, 88, 89], [88, 87, 86, 85], [90, 90, 89, 88]], + [[48, 49, 49, 45], [49, 48, 47, 45], [48, 47, 48, 46]], + ], + dtype=np.uint8, + ), + coords={ + "band": [1, 2, 3], + "x": [-52.5, -51.5, -50.5, -49.5], + "y": [-17.5, -18.5, -19.5], + "spatial_ref": 0, + }, + dims=["band", "y", "x"], + attrs={ + "scale_factor": 1.0, + "add_offset": 0.0, + "_FillValue": 0, + "STATISTICS_MAXIMUM": 95, + "STATISTICS_MEAN": 90.91666666666667, + "STATISTICS_MINIMUM": 89, + "STATISTICS_STDDEV": 1.5523280008498, + "STATISTICS_VALID_PERCENT": 100, + "AREA_OR_POINT": "Area", + }, + ) + + def test_grdcut_dataarray_in_file_out(grid, expected_grid, region): """ Test grdcut on an input DataArray, and output to a grid file. @@ -70,3 +112,22 @@ def test_grdcut_fails(): """ with pytest.raises(GMTInvalidInput): grdcut(np.arange(10).reshape((5, 2))) + + +def test_grdcut_image_file(region, expected_image): + """ + Test grdcut on an input image file. + """ + result = grdcut("@earth_day_01d", region=region) + xr.testing.assert_allclose(a=result, b=expected_image) + + +@pytest.mark.skipif(not _HAS_RIOXARRAY, reason="rioxarray is not installed") +def test_grdcut_image_dataarray(region, expected_image): + """ + Test grdcut on an input xarray.DataArray object. + """ + path = which("@earth_day_01d", download="a") + with rioxarray.open_rasterio(path) as raster: + result = grdcut(raster, region=region) + xr.testing.assert_allclose(a=result, b=expected_image) From 7f6ca7de4dbda799d173096e00bcccf37d5ddc46 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 7 Jul 2024 18:39:13 +0800 Subject: [PATCH 24/73] Fix one failing test --- pygmt/tests/test_grdcut.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pygmt/tests/test_grdcut.py b/pygmt/tests/test_grdcut.py index 0241679e5da..a757e16d969 100644 --- a/pygmt/tests/test_grdcut.py +++ b/pygmt/tests/test_grdcut.py @@ -114,6 +114,7 @@ def test_grdcut_fails(): grdcut(np.arange(10).reshape((5, 2))) +@pytest.mark.skipif(not _HAS_RIOXARRAY, reason="rioxarray is not installed") def test_grdcut_image_file(region, expected_image): """ Test grdcut on an input image file. From 21b194ad7056845470fea216e48ed60e0eb5b967 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 7 Jul 2024 18:48:56 +0800 Subject: [PATCH 25/73] Fix open_rasterio --- pygmt/clib/session.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 0bc111d3aec..35d6ecb960f 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -2089,8 +2089,7 @@ def virtualfile_to_raster( ) with GMTTempFile(suffix=".tif") as tmpfile: self.call_module("write", f"{vfname} {tmpfile.name} -Ti") - with rioxarray.open_rasterio(tmpfile.name) as dataarray: - return dataarray.load() + return rioxarray.open_rasterio(tmpfile.name) # type: ignore[return-value] return self.read_virtualfile(vfname, kind=kind).contents.to_dataarray() def extract_region(self): From e7eaf5cbd8ed67e145c28ffe17239fc2ff4c8179 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 7 Jul 2024 19:09:33 +0800 Subject: [PATCH 26/73] Fix open_rasterio --- pygmt/tests/test_grdcut.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygmt/tests/test_grdcut.py b/pygmt/tests/test_grdcut.py index a757e16d969..73c7437ae76 100644 --- a/pygmt/tests/test_grdcut.py +++ b/pygmt/tests/test_grdcut.py @@ -129,6 +129,6 @@ def test_grdcut_image_dataarray(region, expected_image): Test grdcut on an input xarray.DataArray object. """ path = which("@earth_day_01d", download="a") - with rioxarray.open_rasterio(path) as raster: - result = grdcut(raster, region=region) + raster = rioxarray.open_rasterio(path) + result = grdcut(raster, region=region) xr.testing.assert_allclose(a=result, b=expected_image) From e3c8569577ac96cacc70dd79296ff68d464870b2 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 7 Jul 2024 19:25:56 +0800 Subject: [PATCH 27/73] Make sure the image is loaded --- pygmt/clib/session.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 35d6ecb960f..8d8b0d61ea9 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -2089,7 +2089,8 @@ def virtualfile_to_raster( ) with GMTTempFile(suffix=".tif") as tmpfile: self.call_module("write", f"{vfname} {tmpfile.name} -Ti") - return rioxarray.open_rasterio(tmpfile.name) # type: ignore[return-value] + dataarray = rioxarray.open_rasterio(tmpfile.name) # type: ignore[return-value] + return dataarray.load() return self.read_virtualfile(vfname, kind=kind).contents.to_dataarray() def extract_region(self): From 1c8312cb3550690c6be0dfc940008572652ca818 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 8 Jul 2024 07:02:02 +0800 Subject: [PATCH 28/73] Update pygmt/clib/session.py Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com> --- pygmt/clib/session.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 8d8b0d61ea9..46a09a78b60 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -2089,8 +2089,8 @@ def virtualfile_to_raster( ) with GMTTempFile(suffix=".tif") as tmpfile: self.call_module("write", f"{vfname} {tmpfile.name} -Ti") - dataarray = rioxarray.open_rasterio(tmpfile.name) # type: ignore[return-value] - return dataarray.load() + dataarray = rioxarray.open_rasterio(tmpfile.name).load() + return dataarray return self.read_virtualfile(vfname, kind=kind).contents.to_dataarray() def extract_region(self): From 391343032f3229996ae0acb3c4ff927698b15307 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 8 Jul 2024 09:20:05 +0800 Subject: [PATCH 29/73] Use rioxarray.open_rasterio in a context manager --- pygmt/clib/session.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 46a09a78b60..7dd824983dc 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -2089,7 +2089,8 @@ def virtualfile_to_raster( ) with GMTTempFile(suffix=".tif") as tmpfile: self.call_module("write", f"{vfname} {tmpfile.name} -Ti") - dataarray = rioxarray.open_rasterio(tmpfile.name).load() + with rioxarray.open_rasterio(tmpfile.name) as da: + dataarray = da.load() return dataarray return self.read_virtualfile(vfname, kind=kind).contents.to_dataarray() From ab77187d8bd7fdc2b265dc3da3dc279a17619c6e Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 8 Jul 2024 10:43:11 +0800 Subject: [PATCH 30/73] Fix mypy errors --- pygmt/clib/session.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 7dd824983dc..981d3f92c0b 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -2089,9 +2089,9 @@ def virtualfile_to_raster( ) with GMTTempFile(suffix=".tif") as tmpfile: self.call_module("write", f"{vfname} {tmpfile.name} -Ti") - with rioxarray.open_rasterio(tmpfile.name) as da: - dataarray = da.load() - return dataarray + with rioxarray.open_rasterio(tmpfile.name) as da: # type: ignore[union-attr] + dataarray = da.load() # type: ignore[union-attr] + return dataarray # type: ignore[return-value] return self.read_virtualfile(vfname, kind=kind).contents.to_dataarray() def extract_region(self): From 6f3e4746db66e2c68425d3355f2eaad095559ea6 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 8 Jul 2024 10:45:18 +0800 Subject: [PATCH 31/73] Move grdcut image tests to a separate test file --- pygmt/tests/test_grdcut.py | 29 +-------------------------- pygmt/tests/test_grdcut_image.py | 34 ++++++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+), 28 deletions(-) create mode 100644 pygmt/tests/test_grdcut_image.py diff --git a/pygmt/tests/test_grdcut.py b/pygmt/tests/test_grdcut.py index 73c7437ae76..f0ad5e3c178 100644 --- a/pygmt/tests/test_grdcut.py +++ b/pygmt/tests/test_grdcut.py @@ -5,18 +5,11 @@ import numpy as np import pytest import xarray as xr -from pygmt import grdcut, load_dataarray, which +from pygmt import grdcut, load_dataarray from pygmt.exceptions import GMTInvalidInput from pygmt.helpers import GMTTempFile from pygmt.helpers.testing import load_static_earth_relief -try: - import rioxarray - - _HAS_RIOXARRAY = True -except ImportError: - _HAS_RIOXARRAY = False - @pytest.fixture(scope="module", name="grid") def fixture_grid(): @@ -112,23 +105,3 @@ def test_grdcut_fails(): """ with pytest.raises(GMTInvalidInput): grdcut(np.arange(10).reshape((5, 2))) - - -@pytest.mark.skipif(not _HAS_RIOXARRAY, reason="rioxarray is not installed") -def test_grdcut_image_file(region, expected_image): - """ - Test grdcut on an input image file. - """ - result = grdcut("@earth_day_01d", region=region) - xr.testing.assert_allclose(a=result, b=expected_image) - - -@pytest.mark.skipif(not _HAS_RIOXARRAY, reason="rioxarray is not installed") -def test_grdcut_image_dataarray(region, expected_image): - """ - Test grdcut on an input xarray.DataArray object. - """ - path = which("@earth_day_01d", download="a") - raster = rioxarray.open_rasterio(path) - result = grdcut(raster, region=region) - xr.testing.assert_allclose(a=result, b=expected_image) diff --git a/pygmt/tests/test_grdcut_image.py b/pygmt/tests/test_grdcut_image.py new file mode 100644 index 00000000000..3fc76ffa589 --- /dev/null +++ b/pygmt/tests/test_grdcut_image.py @@ -0,0 +1,34 @@ +""" +Test pygmt.grdcut on images. +""" + +import pytest +import xarray as xr +from pygmt import grdcut, which + +try: + import rioxarray + + _HAS_RIOXARRAY = True +except ImportError: + _HAS_RIOXARRAY = False + + +@pytest.mark.skipif(not _HAS_RIOXARRAY, reason="rioxarray is not installed") +def test_grdcut_image_file(region, expected_image): + """ + Test grdcut on an input image file. + """ + result = grdcut("@earth_day_01d", region=region) + xr.testing.assert_allclose(a=result, b=expected_image) + + +@pytest.mark.skipif(not _HAS_RIOXARRAY, reason="rioxarray is not installed") +def test_grdcut_image_dataarray(region, expected_image): + """ + Test grdcut on an input xarray.DataArray object. + """ + path = which("@earth_day_01d", download="a") + raster = rioxarray.open_rasterio(path) + result = grdcut(raster, region=region) + xr.testing.assert_allclose(a=result, b=expected_image) From 5b07dd99360515e685d781d18d4ef47c35a68304 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 8 Jul 2024 10:55:17 +0800 Subject: [PATCH 32/73] Fix copy & paste errors --- pygmt/tests/test_grdcut.py | 35 ------------------------- pygmt/tests/test_grdcut_image.py | 44 ++++++++++++++++++++++++++++++++ 2 files changed, 44 insertions(+), 35 deletions(-) diff --git a/pygmt/tests/test_grdcut.py b/pygmt/tests/test_grdcut.py index f0ad5e3c178..dbf5dd21f49 100644 --- a/pygmt/tests/test_grdcut.py +++ b/pygmt/tests/test_grdcut.py @@ -43,41 +43,6 @@ def fixture_expected_grid(): ) -@pytest.fixture(scope="module", name="expected_image") -def fixture_expected_image(): - """ - Load the expected grdcut image result. - """ - return xr.DataArray( - data=np.array( - [ - [[90, 93, 95, 90], [91, 90, 91, 91], [91, 90, 89, 90]], - [[87, 88, 88, 89], [88, 87, 86, 85], [90, 90, 89, 88]], - [[48, 49, 49, 45], [49, 48, 47, 45], [48, 47, 48, 46]], - ], - dtype=np.uint8, - ), - coords={ - "band": [1, 2, 3], - "x": [-52.5, -51.5, -50.5, -49.5], - "y": [-17.5, -18.5, -19.5], - "spatial_ref": 0, - }, - dims=["band", "y", "x"], - attrs={ - "scale_factor": 1.0, - "add_offset": 0.0, - "_FillValue": 0, - "STATISTICS_MAXIMUM": 95, - "STATISTICS_MEAN": 90.91666666666667, - "STATISTICS_MINIMUM": 89, - "STATISTICS_STDDEV": 1.5523280008498, - "STATISTICS_VALID_PERCENT": 100, - "AREA_OR_POINT": "Area", - }, - ) - - def test_grdcut_dataarray_in_file_out(grid, expected_grid, region): """ Test grdcut on an input DataArray, and output to a grid file. diff --git a/pygmt/tests/test_grdcut_image.py b/pygmt/tests/test_grdcut_image.py index 3fc76ffa589..31b3c3b7ae5 100644 --- a/pygmt/tests/test_grdcut_image.py +++ b/pygmt/tests/test_grdcut_image.py @@ -2,6 +2,7 @@ Test pygmt.grdcut on images. """ +import numpy as np import pytest import xarray as xr from pygmt import grdcut, which @@ -14,6 +15,49 @@ _HAS_RIOXARRAY = False +@pytest.fixture(scope="module", name="region") +def fixture_region(): + """ + Set the data region. + """ + return [-53, -49, -20, -17] + + +@pytest.fixture(scope="module", name="expected_image") +def fixture_expected_image(): + """ + Load the expected grdcut image result. + """ + return xr.DataArray( + data=np.array( + [ + [[90, 93, 95, 90], [91, 90, 91, 91], [91, 90, 89, 90]], + [[87, 88, 88, 89], [88, 87, 86, 85], [90, 90, 89, 88]], + [[48, 49, 49, 45], [49, 48, 47, 45], [48, 47, 48, 46]], + ], + dtype=np.uint8, + ), + coords={ + "band": [1, 2, 3], + "x": [-52.5, -51.5, -50.5, -49.5], + "y": [-17.5, -18.5, -19.5], + "spatial_ref": 0, + }, + dims=["band", "y", "x"], + attrs={ + "scale_factor": 1.0, + "add_offset": 0.0, + "_FillValue": 0, + "STATISTICS_MAXIMUM": 95, + "STATISTICS_MEAN": 90.91666666666667, + "STATISTICS_MINIMUM": 89, + "STATISTICS_STDDEV": 1.5523280008498, + "STATISTICS_VALID_PERCENT": 100, + "AREA_OR_POINT": "Area", + }, + ) + + @pytest.mark.skipif(not _HAS_RIOXARRAY, reason="rioxarray is not installed") def test_grdcut_image_file(region, expected_image): """ From 31272abe223c4debe3f28f6112a4ebd315c9b101 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 8 Jul 2024 12:26:14 +0800 Subject: [PATCH 33/73] Run codspeed benchmark for test_grdcut_image_dataarray Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com> --- pygmt/tests/test_grdcut_image.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pygmt/tests/test_grdcut_image.py b/pygmt/tests/test_grdcut_image.py index 31b3c3b7ae5..bd6dbfba41b 100644 --- a/pygmt/tests/test_grdcut_image.py +++ b/pygmt/tests/test_grdcut_image.py @@ -67,6 +67,7 @@ def test_grdcut_image_file(region, expected_image): xr.testing.assert_allclose(a=result, b=expected_image) +@pytest.mark.benchmark @pytest.mark.skipif(not _HAS_RIOXARRAY, reason="rioxarray is not installed") def test_grdcut_image_dataarray(region, expected_image): """ From 279595b7a498bc3ad4f0e73ef1b8b16aa2977dc1 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 5 Aug 2024 16:24:58 +0800 Subject: [PATCH 34/73] Add the raster_kind function to determine the raster kind --- pygmt/clib/session.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 82f155d44df..9d70b6e8d95 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -2231,3 +2231,29 @@ def extract_region(self) -> np.ndarray: if status != 0: raise GMTCLibError("Failed to extract region from current figure.") return region + + +def raster_kind(raster: str): + """ + Determine the raster kind. + + >>> raster_kind("@earth_relief_01d") + 'grid' + >>> raster_kind("@static_earth_relief.nc") + 'grid' + >>> raster_kind("@earth_day_01d") + 'image' + >>> raster_kind("@hotspots.txt") + """ + with Session() as lib: + try: + img = lib.read_data(infile=raster, kind="image", mode="GMT_CONTAINER_ONLY") + return "image" if img.contents.header.contents.n_bands == 3 else "grid" + except GMTCLibError: + pass + try: + _ = lib.read_data(infile=raster, kind="grid", mode="GMT_CONTAINER_ONLY") + return "grid" + except GMTCLibError: + pass + return None From 7def4b51d330e8cff341823db8d8a9202197617f Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 5 Aug 2024 16:54:45 +0800 Subject: [PATCH 35/73] Simplify the grdcut function --- pygmt/src/grdcut.py | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index 41e5ac5596d..b12acf34d77 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -3,6 +3,7 @@ """ from pygmt.clib import Session +from pygmt.clib.session import raster_kind from pygmt.helpers import ( build_arg_list, data_kind, @@ -10,7 +11,6 @@ kwargs_to_strings, use_alias, ) -from pygmt.src.which import which __doctest_skip__ = ["grdcut"] @@ -104,18 +104,7 @@ def grdcut(grid, outgrid: str | None = None, **kwargs): case "image" | "grid": outkind = inkind case "file": - realpath = str(grid) - if realpath.startswith("@"): # Is a remote file without suffix - realpath = which(grid, download="a") - if isinstance(realpath, list): - realpath = realpath[0] - - if realpath.endswith( - (".jpg", ".jpeg", ".png", ".tif", ".tiff", ".bmp", ".webp") - ): - outkind = "image" - else: # Fall back to grid. - outkind = "grid" + outkind = raster_kind(grid) with Session() as lib: with ( From 7d437be0675d7a986bafa937c8fe19557d7bcce6 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 00:34:58 +0800 Subject: [PATCH 36/73] 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 37/73] 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 38/73] 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 39/73] 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 40/73] 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 41/73] 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 42/73] 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 43/73] 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 44/73] 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 45/73] 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 46/73] 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 47/73] 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 48/73] 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 49/73] 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 bf9275c1534e4ad2a659f6d7c887d97a0b8727e2 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 16:43:28 +0800 Subject: [PATCH 50/73] Remove the workaround for images --- pygmt/clib/session.py | 21 --------------------- 1 file changed, 21 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 97a56bbf52a..36935807081 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -30,20 +30,12 @@ from pygmt.datatypes import _GMT_DATASET, _GMT_GRID, _GMT_IMAGE from pygmt.exceptions import GMTCLibError, GMTCLibNoSessionError, GMTInvalidInput from pygmt.helpers import ( - GMTTempFile, _validate_data_input, data_kind, tempfile_from_geojson, tempfile_from_image, ) -try: - import rioxarray - - _HAS_RIOXARRAY = True -except ImportError: - _HAS_RIOXARRAY = False - FAMILIES = [ "GMT_IS_DATASET", # Entity is a data table "GMT_IS_GRID", # Entity is a grid @@ -2258,19 +2250,6 @@ def virtualfile_to_raster( self["GMT_IS_IMAGE"]: "image", self["GMT_IS_CUBE"]: "cube", }[family] - - if kind == "image": # Use temporary file for images - if not _HAS_RIOXARRAY: - raise ImportError( - "Package `rioxarray` is required to be installed to load images. " - "Please use `python -m pip install rioxarray` or " - "`mamba install -c conda-forge rioxarray` to install the package." - ) - with GMTTempFile(suffix=".tif") as tmpfile: - self.call_module("write", f"{vfname} {tmpfile.name} -Ti") - with rioxarray.open_rasterio(tmpfile.name) as da: # type: ignore[union-attr] - dataarray = da.load() # type: ignore[union-attr] - return dataarray # type: ignore[return-value] return self.read_virtualfile(vfname, kind=kind).contents.to_dataarray() def extract_region(self) -> np.ndarray: From fb97daa9f032295a94736770d64c480463ca4812 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 20 Sep 2024 19:58:52 +0800 Subject: [PATCH 51/73] 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 52/73] 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 53/73] 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 54/73] 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 55/73] 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 56/73] 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 a97d0b3d779fb964f6a58a8179ce44ac847e536b Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sat, 28 Sep 2024 15:28:13 +0800 Subject: [PATCH 57/73] Apply suggestions from code review Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com> --- pygmt/clib/session.py | 12 ++++++------ pygmt/src/grdcut.py | 6 +++--- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 36935807081..f11084a9cff 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -1942,7 +1942,7 @@ def virtualfile_out( "grid": ("GMT_IS_GRID", "GMT_IS_SURFACE"), "image": ("GMT_IS_IMAGE", "GMT_IS_SURFACE"), }[kind] - # For unkown reasons, 'GMT_OUT' crashes for 'image' kind. + # For unknown reasons, 'GMT_OUT' crashes for 'image' kind. direction = "GMT_OUT|GMT_IS_REFERENCE" if kind == "image" else "GMT_OUT" with self.open_virtualfile(family, geometry, direction, None) as vfile: yield vfile @@ -2314,17 +2314,17 @@ def extract_region(self) -> np.ndarray: return region -def raster_kind(raster: str): +def _raster_kind(raster: str) -> Literal["grid", "image"] | None: """ Determine the raster kind. - >>> raster_kind("@earth_relief_01d") + >>> _raster_kind("@earth_relief_01d") 'grid' - >>> raster_kind("@static_earth_relief.nc") + >>> _raster_kind("@static_earth_relief.nc") 'grid' - >>> raster_kind("@earth_day_01d") + >>> _raster_kind("@earth_day_01d") 'image' - >>> raster_kind("@hotspots.txt") + >>> _raster_kind("@hotspots.txt") """ with Session() as lib: try: diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index e8ddb931c7d..5b87788c6e7 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -4,7 +4,7 @@ import xarray as xr from pygmt.clib import Session -from pygmt.clib.session import raster_kind +from pygmt.clib.session import _raster_kind from pygmt.helpers import ( build_arg_list, data_kind, @@ -105,7 +105,7 @@ def grdcut(grid, outgrid: str | None = None, **kwargs) -> xr.DataArray | None: case "image" | "grid": outkind = inkind case "file": - outkind = raster_kind(grid) + outkind = _raster_kind(grid) with Session() as lib: with ( @@ -115,5 +115,5 @@ def grdcut(grid, outgrid: str | None = None, **kwargs) -> xr.DataArray | None: kwargs["G"] = voutgrd lib.call_module(module="grdcut", args=build_arg_list(kwargs, infile=vingrd)) return lib.virtualfile_to_raster( - outgrid=outgrid, kind=outkind, vfname=voutgrd + vfname=voutgrd, kind=outkind, outgrid=outgrid ) From 2fd13fb7b8efd9b8438440b84199244a71397e05 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sat, 28 Sep 2024 21:25:44 +0800 Subject: [PATCH 58/73] 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 59/73] 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 9ec00be8f6e8574cb6c28d14adf98fa4de6bb9de Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 29 Sep 2024 21:34:24 +0800 Subject: [PATCH 60/73] Let _raster_kind return grid by default --- pygmt/clib/session.py | 10 ++++++++-- pygmt/src/grdcut.py | 4 ++++ 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 7e32aaf0fbf..a15f3fa711c 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -2314,10 +2314,12 @@ def extract_region(self) -> np.ndarray: return region -def _raster_kind(raster: str) -> Literal["grid", "image"] | None: +def _raster_kind(raster: str) -> Literal["grid", "image"]: """ Determine the raster kind. + Examples + -------- >>> _raster_kind("@earth_relief_01d") 'grid' >>> _raster_kind("@static_earth_relief.nc") @@ -2325,7 +2327,11 @@ def _raster_kind(raster: str) -> Literal["grid", "image"] | None: >>> _raster_kind("@earth_day_01d") 'image' >>> _raster_kind("@hotspots.txt") + 'grid' """ + # The logic here is because: an image can be read into a grid container, but a grid + # can't be read into an image container. So, try to read the file as an image first. + # If fails, try to read it as a grid. with Session() as lib: try: img = lib.read_data(infile=raster, kind="image", mode="GMT_CONTAINER_ONLY") @@ -2337,4 +2343,4 @@ def _raster_kind(raster: str) -> Literal["grid", "image"] | None: return "grid" except GMTCLibError: pass - return None + return "grid" # Fallback to "grid" and let GMT determine the type. diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index 5b87788c6e7..7f33b6f2784 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -5,6 +5,7 @@ import xarray as xr from pygmt.clib import Session from pygmt.clib.session import _raster_kind +from pygmt.exceptions import GMTInvalidInput from pygmt.helpers import ( build_arg_list, data_kind, @@ -106,6 +107,9 @@ def grdcut(grid, outgrid: str | None = None, **kwargs) -> xr.DataArray | None: outkind = inkind case "file": outkind = _raster_kind(grid) + case "_": + msg = f"Unsupported data type {type(grid)}." + raise GMTInvalidInput(msg) with Session() as lib: with ( From f3a2f8e6f7ff17d8e7c1db85e800058977f10dce Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 29 Sep 2024 21:37:41 +0800 Subject: [PATCH 61/73] Simplify the grdcut image tests --- pygmt/tests/test_grdcut_image.py | 16 +++------------- 1 file changed, 3 insertions(+), 13 deletions(-) diff --git a/pygmt/tests/test_grdcut_image.py b/pygmt/tests/test_grdcut_image.py index bd6dbfba41b..e25e2264cfb 100644 --- a/pygmt/tests/test_grdcut_image.py +++ b/pygmt/tests/test_grdcut_image.py @@ -41,39 +41,29 @@ def fixture_expected_image(): "band": [1, 2, 3], "x": [-52.5, -51.5, -50.5, -49.5], "y": [-17.5, -18.5, -19.5], - "spatial_ref": 0, }, dims=["band", "y", "x"], attrs={ "scale_factor": 1.0, "add_offset": 0.0, - "_FillValue": 0, - "STATISTICS_MAXIMUM": 95, - "STATISTICS_MEAN": 90.91666666666667, - "STATISTICS_MINIMUM": 89, - "STATISTICS_STDDEV": 1.5523280008498, - "STATISTICS_VALID_PERCENT": 100, - "AREA_OR_POINT": "Area", }, ) -@pytest.mark.skipif(not _HAS_RIOXARRAY, reason="rioxarray is not installed") +@pytest.mark.benchmark def test_grdcut_image_file(region, expected_image): """ Test grdcut on an input image file. """ - result = grdcut("@earth_day_01d", region=region) + result = grdcut("@earth_day_01d_p", region=region) xr.testing.assert_allclose(a=result, b=expected_image) -@pytest.mark.benchmark @pytest.mark.skipif(not _HAS_RIOXARRAY, reason="rioxarray is not installed") def test_grdcut_image_dataarray(region, expected_image): """ Test grdcut on an input xarray.DataArray object. """ - path = which("@earth_day_01d", download="a") - raster = rioxarray.open_rasterio(path) + raster = rioxarray.open_rasterio(which("@earth_day_01d", download="a")).load() result = grdcut(raster, region=region) xr.testing.assert_allclose(a=result, b=expected_image) From 3c12e2b786b376cd46759694a70f1e5816143ca4 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 29 Sep 2024 21:48:10 +0800 Subject: [PATCH 62/73] Add one more test for file in & file out --- pygmt/tests/test_grdcut_image.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/pygmt/tests/test_grdcut_image.py b/pygmt/tests/test_grdcut_image.py index e25e2264cfb..2e34a5a6465 100644 --- a/pygmt/tests/test_grdcut_image.py +++ b/pygmt/tests/test_grdcut_image.py @@ -2,10 +2,13 @@ Test pygmt.grdcut on images. """ +from pathlib import Path + import numpy as np import pytest import xarray as xr from pygmt import grdcut, which +from pygmt.helpers import GMTTempFile try: import rioxarray @@ -67,3 +70,16 @@ def test_grdcut_image_dataarray(region, expected_image): raster = rioxarray.open_rasterio(which("@earth_day_01d", download="a")).load() result = grdcut(raster, region=region) xr.testing.assert_allclose(a=result, b=expected_image) + + +def test_grdcut_image_file_in_file_out(region, expected_image): + """ + Test grdcut on an input image file and outputs to another image file. + """ + with GMTTempFile(suffix=".tif") as tmp: + result = grdcut("@earth_day_01d_p", region=region, outgrid=tmp.filename) + assert result is None + assert Path(tmp.filename).stat().st_size > 0 + if _HAS_RIOXARRAY: + raster = rioxarray.open_rasterio(which("@earth_day_01d", download="a")).load() + xr.testing.assert_allclose(a=raster, b=expected_image) From f852b0d1a8e2d1ed8802f359de6a3e625377ef43 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Sun, 29 Sep 2024 22:12:02 +0800 Subject: [PATCH 63/73] Fix typos --- pygmt/tests/test_grdcut_image.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygmt/tests/test_grdcut_image.py b/pygmt/tests/test_grdcut_image.py index 2e34a5a6465..b5bc433af9e 100644 --- a/pygmt/tests/test_grdcut_image.py +++ b/pygmt/tests/test_grdcut_image.py @@ -77,9 +77,9 @@ def test_grdcut_image_file_in_file_out(region, expected_image): Test grdcut on an input image file and outputs to another image file. """ with GMTTempFile(suffix=".tif") as tmp: - result = grdcut("@earth_day_01d_p", region=region, outgrid=tmp.filename) + result = grdcut("@earth_day_01d_p", region=region, outgrid=tmp.name) assert result is None - assert Path(tmp.filename).stat().st_size > 0 + assert Path(tmp.name).stat().st_size > 0 if _HAS_RIOXARRAY: raster = rioxarray.open_rasterio(which("@earth_day_01d", download="a")).load() xr.testing.assert_allclose(a=raster, b=expected_image) From bb1a0b0032770c44091105d461cee2febe6cc07c Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 30 Sep 2024 13:08:24 +0800 Subject: [PATCH 64/73] Use the new load_blue_marble function --- pygmt/tests/test_grdcut_image.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/pygmt/tests/test_grdcut_image.py b/pygmt/tests/test_grdcut_image.py index b5bc433af9e..fe26f3d1bc0 100644 --- a/pygmt/tests/test_grdcut_image.py +++ b/pygmt/tests/test_grdcut_image.py @@ -7,7 +7,8 @@ import numpy as np import pytest import xarray as xr -from pygmt import grdcut, which +from pygmt import grdcut +from pygmt.datasets import load_blue_marble from pygmt.helpers import GMTTempFile try: @@ -67,7 +68,7 @@ def test_grdcut_image_dataarray(region, expected_image): """ Test grdcut on an input xarray.DataArray object. """ - raster = rioxarray.open_rasterio(which("@earth_day_01d", download="a")).load() + raster = load_blue_marble() result = grdcut(raster, region=region) xr.testing.assert_allclose(a=result, b=expected_image) @@ -80,6 +81,6 @@ def test_grdcut_image_file_in_file_out(region, expected_image): result = grdcut("@earth_day_01d_p", region=region, outgrid=tmp.name) assert result is None assert Path(tmp.name).stat().st_size > 0 - if _HAS_RIOXARRAY: - raster = rioxarray.open_rasterio(which("@earth_day_01d", download="a")).load() - xr.testing.assert_allclose(a=raster, b=expected_image) + if _HAS_RIOXARRAY: + raster = rioxarray.open_rasterio(tmp.name).load() + xr.testing.assert_allclose(a=raster, b=expected_image) From 584b5af83cc0e3999b67482858f612b43e39defa Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 30 Sep 2024 13:20:16 +0800 Subject: [PATCH 65/73] Drop the spatial_ref coord --- pygmt/tests/test_grdcut_image.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pygmt/tests/test_grdcut_image.py b/pygmt/tests/test_grdcut_image.py index fe26f3d1bc0..1a0ccaadc0c 100644 --- a/pygmt/tests/test_grdcut_image.py +++ b/pygmt/tests/test_grdcut_image.py @@ -59,7 +59,7 @@ def test_grdcut_image_file(region, expected_image): """ Test grdcut on an input image file. """ - result = grdcut("@earth_day_01d_p", region=region) + result = grdcut("@earth_day_01d", region=region) xr.testing.assert_allclose(a=result, b=expected_image) @@ -78,9 +78,9 @@ def test_grdcut_image_file_in_file_out(region, expected_image): Test grdcut on an input image file and outputs to another image file. """ with GMTTempFile(suffix=".tif") as tmp: - result = grdcut("@earth_day_01d_p", region=region, outgrid=tmp.name) + result = grdcut("@earth_day_01d", region=region, outgrid=tmp.name) assert result is None assert Path(tmp.name).stat().st_size > 0 if _HAS_RIOXARRAY: - raster = rioxarray.open_rasterio(tmp.name).load() + raster = rioxarray.open_rasterio(tmp.name).load().drop_vars("spatial_ref") xr.testing.assert_allclose(a=raster, b=expected_image) From fac51f1cc5b3d64595b9762477027bc6a14faf97 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 25 Nov 2024 19:17:33 +0800 Subject: [PATCH 66/73] Update _raster_kind --- pygmt/clib/session.py | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 2b983115548..8294bacc571 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -206,6 +206,14 @@ def info(self) -> dict[str, str]: } return self._info + def __init__(self, silent: bool = False): + """ + Initialize a GMT session. + + Does nothing but set the Session variables. + """ + self._silent = silent + def __enter__(self): """ Create a GMT API session. @@ -382,6 +390,8 @@ def print_func(file_pointer, message): # noqa: ARG001 We'll capture the messages and print them to stderr so that they will show up on the Jupyter notebook. """ + if self._silent: + return 0 # Have to use try..except due to upstream GMT bug in GMT <= 6.5.0. # See https://github.com/GenericMappingTools/pygmt/issues/3205. try: @@ -2401,10 +2411,13 @@ def _raster_kind(raster: str) -> Literal["grid", "image"]: # The logic here is because: an image can be read into a grid container, but a grid # can't be read into an image container. So, try to read the file as an image first. # If fails, try to read it as a grid. - with Session() as lib: + with Session(silent=True) as lib: + from osgeo import gdal + + gdal.PushErrorHandler("CPLQuietErrorHandler") try: - img = lib.read_data(infile=raster, kind="image", mode="GMT_CONTAINER_ONLY") - return "image" if img.contents.header.contents.n_bands == 3 else "grid" + _ = lib.read_data(infile=raster, kind="image", mode="GMT_CONTAINER_ONLY") + return "image" except GMTCLibError: pass try: @@ -2412,4 +2425,5 @@ def _raster_kind(raster: str) -> Literal["grid", "image"]: return "grid" except GMTCLibError: pass - return "grid" # Fallback to "grid" and let GMT determine the type. + gdal.PopErrorHandler() + return "grid" # Fallback to "grid". From 517ee5225ec8f9c7f8b226d4e70c21a4a2a783cb Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 9 Dec 2024 18:51:27 +0800 Subject: [PATCH 67/73] Revert "Update _raster_kind" This reverts commit fac51f1cc5b3d64595b9762477027bc6a14faf97. --- pygmt/clib/session.py | 22 ++++------------------ 1 file changed, 4 insertions(+), 18 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 8294bacc571..2b983115548 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -206,14 +206,6 @@ def info(self) -> dict[str, str]: } return self._info - def __init__(self, silent: bool = False): - """ - Initialize a GMT session. - - Does nothing but set the Session variables. - """ - self._silent = silent - def __enter__(self): """ Create a GMT API session. @@ -390,8 +382,6 @@ def print_func(file_pointer, message): # noqa: ARG001 We'll capture the messages and print them to stderr so that they will show up on the Jupyter notebook. """ - if self._silent: - return 0 # Have to use try..except due to upstream GMT bug in GMT <= 6.5.0. # See https://github.com/GenericMappingTools/pygmt/issues/3205. try: @@ -2411,13 +2401,10 @@ def _raster_kind(raster: str) -> Literal["grid", "image"]: # The logic here is because: an image can be read into a grid container, but a grid # can't be read into an image container. So, try to read the file as an image first. # If fails, try to read it as a grid. - with Session(silent=True) as lib: - from osgeo import gdal - - gdal.PushErrorHandler("CPLQuietErrorHandler") + with Session() as lib: try: - _ = lib.read_data(infile=raster, kind="image", mode="GMT_CONTAINER_ONLY") - return "image" + img = lib.read_data(infile=raster, kind="image", mode="GMT_CONTAINER_ONLY") + return "image" if img.contents.header.contents.n_bands == 3 else "grid" except GMTCLibError: pass try: @@ -2425,5 +2412,4 @@ def _raster_kind(raster: str) -> Literal["grid", "image"]: return "grid" except GMTCLibError: pass - gdal.PopErrorHandler() - return "grid" # Fallback to "grid". + return "grid" # Fallback to "grid" and let GMT determine the type. From 7ca58fb0196a9d72c7e17b24c1cac6df23a1cb4f Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 9 Dec 2024 18:52:45 +0800 Subject: [PATCH 68/73] Remove the _raster_kind function --- pygmt/clib/session.py | 33 --------------------------------- 1 file changed, 33 deletions(-) diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index ac600fb0991..ee37c55d59f 100644 --- a/pygmt/clib/session.py +++ b/pygmt/clib/session.py @@ -2017,7 +2017,6 @@ def virtualfile_out( "grid": ("GMT_IS_GRID", "GMT_IS_SURFACE"), "image": ("GMT_IS_IMAGE", "GMT_IS_SURFACE"), }[kind] - # For unknown reasons, 'GMT_OUT' crashes for 'image' kind. direction = "GMT_OUT|GMT_IS_REFERENCE" if kind == "image" else "GMT_OUT" with self.open_virtualfile(family, geometry, direction, None) as vfile: yield vfile @@ -2389,35 +2388,3 @@ def extract_region(self) -> np.ndarray: msg = "Failed to extract region from current figure." raise GMTCLibError(msg) return region - - -def _raster_kind(raster: str) -> Literal["grid", "image"]: - """ - Determine the raster kind. - - Examples - -------- - >>> _raster_kind("@earth_relief_01d") - 'grid' - >>> _raster_kind("@static_earth_relief.nc") - 'grid' - >>> _raster_kind("@earth_day_01d") - 'image' - >>> _raster_kind("@hotspots.txt") - 'grid' - """ - # The logic here is because: an image can be read into a grid container, but a grid - # can't be read into an image container. So, try to read the file as an image first. - # If fails, try to read it as a grid. - with Session() as lib: - try: - img = lib.read_data(infile=raster, kind="image", mode="GMT_CONTAINER_ONLY") - return "image" if img.contents.header.contents.n_bands == 3 else "grid" - except GMTCLibError: - pass - try: - _ = lib.read_data(infile=raster, kind="grid", mode="GMT_CONTAINER_ONLY") - return "grid" - except GMTCLibError: - pass - return "grid" # Fallback to "grid" and let GMT determine the type. From 66580d48bd73920260247dfdce4af490d9aeff8a Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 9 Dec 2024 18:57:05 +0800 Subject: [PATCH 69/73] Add the 'kind' parameter --- pygmt/src/grdcut.py | 13 ++++++++++--- pygmt/tests/test_grdcut_image.py | 4 ++-- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index 21996b14777..ba2e84617e9 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -2,9 +2,10 @@ grdcut - Extract subregion from a grid. """ +from typing import Literal + import xarray as xr from pygmt.clib import Session -from pygmt.clib.session import _raster_kind from pygmt.exceptions import GMTInvalidInput from pygmt.helpers import ( build_arg_list, @@ -28,7 +29,9 @@ f="coltypes", ) @kwargs_to_strings(R="sequence") -def grdcut(grid, outgrid: str | None = None, **kwargs) -> xr.DataArray | None: +def grdcut( + grid, kind: Literal["grid", "image"] = "grid", outgrid: str | None = None, **kwargs +) -> xr.DataArray | None: r""" Extract subregion from a grid or image. @@ -48,6 +51,10 @@ def grdcut(grid, outgrid: str | None = None, **kwargs) -> xr.DataArray | None: Parameters ---------- {grid} + kind + The raster data kind. Valid values are ``grid`` and ``image``. When the input + ``grid`` is a file name, it's hard to determine if the file is a grid or an + image, so we need to specify the kind explicitly. The default is ``grid``. {outgrid} {projection} {region} @@ -106,7 +113,7 @@ def grdcut(grid, outgrid: str | None = None, **kwargs) -> xr.DataArray | None: case "image" | "grid": outkind = inkind case "file": - outkind = _raster_kind(grid) + outkind = kind case "_": msg = f"Unsupported data type {type(grid)}." raise GMTInvalidInput(msg) diff --git a/pygmt/tests/test_grdcut_image.py b/pygmt/tests/test_grdcut_image.py index 1a0ccaadc0c..a01c5227a4d 100644 --- a/pygmt/tests/test_grdcut_image.py +++ b/pygmt/tests/test_grdcut_image.py @@ -59,7 +59,7 @@ def test_grdcut_image_file(region, expected_image): """ Test grdcut on an input image file. """ - result = grdcut("@earth_day_01d", region=region) + result = grdcut("@earth_day_01d", region=region, kind="image") xr.testing.assert_allclose(a=result, b=expected_image) @@ -69,7 +69,7 @@ def test_grdcut_image_dataarray(region, expected_image): Test grdcut on an input xarray.DataArray object. """ raster = load_blue_marble() - result = grdcut(raster, region=region) + result = grdcut(raster, region=region, kind="image") xr.testing.assert_allclose(a=result, b=expected_image) From 194ee90599b160a8fa4a626bd26f07bf55dfd815 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 9 Dec 2024 19:09:26 +0800 Subject: [PATCH 70/73] Minor update --- pygmt/src/grdcut.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index ba2e84617e9..04ae714d3ca 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -108,8 +108,7 @@ def grdcut( >>> new_grid = pygmt.grdcut(grid=grid, region=[12, 15, 21, 24]) """ # Determine the output data kind based on the input data kind. - inkind = data_kind(grid) - match inkind: + match inkind := data_kind(grid): case "image" | "grid": outkind = inkind case "file": From 473428ced7e0474aa8dffad21545a9e85e3569b5 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 9 Dec 2024 20:49:22 +0800 Subject: [PATCH 71/73] Avoid keep the file open --- pygmt/tests/test_grdcut_image.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pygmt/tests/test_grdcut_image.py b/pygmt/tests/test_grdcut_image.py index a01c5227a4d..0fa90499050 100644 --- a/pygmt/tests/test_grdcut_image.py +++ b/pygmt/tests/test_grdcut_image.py @@ -82,5 +82,6 @@ def test_grdcut_image_file_in_file_out(region, expected_image): assert result is None assert Path(tmp.name).stat().st_size > 0 if _HAS_RIOXARRAY: - raster = rioxarray.open_rasterio(tmp.name).load().drop_vars("spatial_ref") - xr.testing.assert_allclose(a=raster, b=expected_image) + with rioxarray.open_rasterio(tmp.name) as raster: + image = raster.load().drop_vars("spatial_ref") + xr.testing.assert_allclose(a=image, b=expected_image) From 84b2ed45bf35238759259c52010946193bb3b6af Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 9 Dec 2024 22:33:39 +0800 Subject: [PATCH 72/73] Remove one unnecessary pytest.skipif marker --- pygmt/tests/test_grdcut_image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygmt/tests/test_grdcut_image.py b/pygmt/tests/test_grdcut_image.py index 0fa90499050..fc0417154c3 100644 --- a/pygmt/tests/test_grdcut_image.py +++ b/pygmt/tests/test_grdcut_image.py @@ -63,7 +63,7 @@ def test_grdcut_image_file(region, expected_image): xr.testing.assert_allclose(a=result, b=expected_image) -@pytest.mark.skipif(not _HAS_RIOXARRAY, reason="rioxarray is not installed") +@pytest.mark.benchmark def test_grdcut_image_dataarray(region, expected_image): """ Test grdcut on an input xarray.DataArray object. From 796809680c4bd51ef4c9bee229511cedf031a240 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Mon, 9 Dec 2024 22:46:52 +0800 Subject: [PATCH 73/73] Add it back because we still need rioxarray --- pygmt/tests/test_grdcut_image.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pygmt/tests/test_grdcut_image.py b/pygmt/tests/test_grdcut_image.py index fc0417154c3..585a7b1b4a5 100644 --- a/pygmt/tests/test_grdcut_image.py +++ b/pygmt/tests/test_grdcut_image.py @@ -64,6 +64,7 @@ def test_grdcut_image_file(region, expected_image): @pytest.mark.benchmark +@pytest.mark.skipif(not _HAS_RIOXARRAY, reason="rioxarray is not installed") def test_grdcut_image_dataarray(region, expected_image): """ Test grdcut on an input xarray.DataArray object.