Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 42 commits into from
Sep 29, 2024
Merged
Changes from 5 commits
Commits
Show all changes
42 commits
Select commit Hold shift + click to select a range
bcf43f0
Wrap GMT's standard data type GMT_IMAGE for images
seisman Mar 18, 2024
a052a1a
Initial implementation of to_dataarray method for _GMT_IMAGE class
weiji14 Mar 20, 2024
56a6d65
Merge branch 'main' into datatypes/gmtimage
seisman Apr 17, 2024
f71e79c
Merge branch 'main' into datatypes/gmtimage
weiji14 Jun 18, 2024
4cce4a2
Small typo fixes and add output type-hint for to_dataarray
weiji14 Jun 18, 2024
e02b650
Fix mypy error using np.array([0, 1, 2]) instead of np.arange
weiji14 Jun 18, 2024
f3d4b1f
Parse name and data_attrs from grid/image header
weiji14 Jun 18, 2024
4390136
Transpose array to (band, y, x) order and add doctest for to_dataarray
weiji14 Jun 20, 2024
5f25669
Set registration and gtype from header
weiji14 Jun 20, 2024
a3c6c14
Print basic shape and padding info in _GMT_IMAGE doctest
weiji14 Jun 20, 2024
5888e10
Only set Conventions = CF-1.7 attribute for NetCDF grid type
weiji14 Jun 20, 2024
798e658
Merge branch 'main' into datatypes/gmtimage
weiji14 Jun 20, 2024
3dbf2f2
Remove rioxarray import
weiji14 Jun 20, 2024
6b860bf
Merge branch 'main' into datatypes/gmtimage
seisman Jul 27, 2024
0bf9368
Merge branch 'main' into datatypes/gmtimage
seisman Sep 19, 2024
7d437be
Use enum for grid ids
seisman Sep 19, 2024
268e34e
Fix the band. Starting from 1
seisman Sep 19, 2024
86765e1
Refactor the tests for images
seisman Sep 19, 2024
86f3ffa
In np.reshape, a is a position-only parameter
seisman Sep 20, 2024
cc28247
Improve tests
seisman Sep 20, 2024
1e2c973
Fix one failing doctest due to xarray changes
seisman Sep 20, 2024
734dc28
The np.reshape's newshape parameter is deprecated
seisman Sep 20, 2024
919dc00
Define grid IDs using IntEnum instead of Enum
seisman Sep 20, 2024
b1eacf1
Pass the new shape as a positional parameter
seisman Sep 20, 2024
aa4fdc9
Fix failing tests
seisman Sep 20, 2024
c87a3ec
One more fix
seisman Sep 20, 2024
a20d8a2
One more fix
seisman Sep 20, 2024
926427b
Simplify a doctest
seisman Sep 20, 2024
c73328e
Improve the tests
seisman Sep 20, 2024
fb97daa
Convert ctypes array to numpy array using np.ctypeslib.as_array
seisman Sep 20, 2024
15b8d53
Fix the incorrect value due to floating number conversion in sphinter…
seisman Sep 20, 2024
8433e78
Merge branch 'ctypesarray' into datatypes/gmtimage
seisman Sep 20, 2024
3e3a6f3
Update the to_dataarray method to match the codes in GMT_GRID
seisman Sep 20, 2024
12ef40a
image data should has uint8 dtype
seisman Sep 20, 2024
f64fbb8
Further improve the tests
seisman Sep 21, 2024
4f2ae48
Merge branch 'main' into datatypes/gmtimage
seisman Sep 24, 2024
d49afed
Add a note that currently only 3-band images are supported
seisman Sep 24, 2024
f70bec0
Merge branch 'main' into datatypes/gmtimage
seisman Sep 28, 2024
2fd13fb
Remove the old GMTGridID enums from pygmt/datatypes/header.py
seisman Sep 28, 2024
9972ba1
A minor fix
seisman Sep 28, 2024
62f0ce0
Set CF-1.7-specific attributes for netCDF formats only
seisman Sep 29, 2024
f715aee
Improve the logic for determine grid/image gtype based on upstream GM…
seisman Sep 29, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 14 additions & 8 deletions pygmt/clib/session.py
Original file line number Diff line number Diff line change
@@ -25,7 +25,7 @@
vectors_to_arrays,
)
from pygmt.clib.loading import load_libgmt
from pygmt.datatypes import _GMT_DATASET, _GMT_GRID
from pygmt.datatypes import _GMT_DATASET, _GMT_GRID, _GMT_IMAGE
from pygmt.exceptions import (
GMTCLibError,
GMTCLibNoSessionError,
@@ -1697,7 +1697,9 @@ def virtualfile_from_data(

@contextlib.contextmanager
def virtualfile_out(
self, kind: Literal["dataset", "grid"] = "dataset", fname: str | None = None
self,
kind: Literal["dataset", "grid", "image"] = "dataset",
fname: str | None = None,
):
r"""
Create a virtual file or an actual file for storing output data.
@@ -1710,8 +1712,8 @@ def virtualfile_out(
Parameters
----------
kind
The data kind of the virtual file to create. Valid values are ``"dataset"``
and ``"grid"``. Ignored if ``fname`` is specified.
The data kind of the virtual file to create. Valid values are ``"dataset"``,
``"grid"``, and ``"image"``. Ignored if ``fname`` is specified.
fname
The name of the actual file to write the output data. No virtual file will
be created.
@@ -1754,8 +1756,11 @@ def virtualfile_out(
family, geometry = {
"dataset": ("GMT_IS_DATASET", "GMT_IS_PLP"),
"grid": ("GMT_IS_GRID", "GMT_IS_SURFACE"),
"image": ("GMT_IS_IMAGE", "GMT_IS_SURFACE"),
}[kind]
with self.open_virtualfile(family, geometry, "GMT_OUT", None) as vfile:
with self.open_virtualfile(
family, geometry, "GMT_OUT|GMT_IS_REFERENCE", None
seisman marked this conversation as resolved.
Show resolved Hide resolved
seisman marked this conversation as resolved.
Show resolved Hide resolved
) as vfile:
yield vfile

def inquire_virtualfile(self, vfname: str) -> int:
@@ -1801,7 +1806,8 @@ def read_virtualfile(
Name of the virtual file to read.
kind
Cast the data into a GMT data container. Valid values are ``"dataset"``,
``"grid"`` and ``None``. If ``None``, will return a ctypes void pointer.
``"grid"``, ``"image"`` and ``None``. If ``None``, will return a ctypes void
pointer.

Examples
--------
@@ -1849,9 +1855,9 @@ def read_virtualfile(
# _GMT_DATASET).
if kind is None: # Return the ctypes void pointer
return pointer
if kind in ["image", "cube"]:
if kind == "cube":
raise NotImplementedError(f"kind={kind} is not supported yet.")
dtype = {"dataset": _GMT_DATASET, "grid": _GMT_GRID}[kind]
dtype = {"dataset": _GMT_DATASET, "grid": _GMT_GRID, "image": _GMT_IMAGE}[kind]
return ctp.cast(pointer, ctp.POINTER(dtype))

def virtualfile_to_dataset(
1 change: 1 addition & 0 deletions pygmt/datatypes/__init__.py
Original file line number Diff line number Diff line change
@@ -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
102 changes: 102 additions & 0 deletions pygmt/datatypes/image.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,102 @@
"""
seisman marked this conversation as resolved.
Show resolved Hide resolved
Wrapper for the GMT_IMAGE data type.
"""

import ctypes as ctp
from typing import ClassVar

import numpy as np
import xarray as xr
from pygmt.datatypes.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]
>>> da = xr.DataArray(
... data=data,
... dims=["y", "x", "band"],
... coords={"y": y, "x": x, "band": [1, 2, 3]},
... )
>>> da = da.transpose("band", "y", "x")
>>> da = da.sortby(list(data.dims))
>>> da.plot.imshow()
"""

_fields_: ClassVar = [
# Data type, e.g. GMT_FLOAT
("type", ctp.c_int),
# Array with color lookup values
("colormap", ctp.POINTER(ctp.c_int)),
# Number of colors in a paletted image
("n_indexed_colors", ctp.c_int),
# Pointer to full GMT header for the image
("header", ctp.POINTER(_GMT_GRID_HEADER)),
# Pointer to actual image
("data", ctp.POINTER(ctp.c_ubyte)),
# Pointer to an optional transparency layer stored in a separate variable
("alpha", ctp.POINTER(ctp.c_ubyte)),
# Color interpolation
("color_interp", ctp.c_char_p),
# Pointer to the x-coordinate vector
("x", ctp.POINTER(ctp.c_double)),
# Pointer to the y-coordinate vector
("y", ctp.POINTER(ctp.c_double)),
# Book-keeping variables "hidden" from the API
("hidden", ctp.c_void_p),
]

def to_dataarray(self) -> xr.DataArray:
"""
Convert a _GMT_IMAGE object to an :class:`xarray.DataArray` object.

Returns
-------
dataarray
A :class:`xarray.DataArray` object.
"""

# Get image header
header: _GMT_GRID_HEADER = self.header.contents

Check warning on line 83 in pygmt/datatypes/image.py

Codecov / codecov/patch

pygmt/datatypes/image.py#L83

Added line #L83 was not covered by tests

# Get DataArray without padding
pad = header.pad[:]
data: np.ndarray = np.reshape(

Check warning on line 87 in pygmt/datatypes/image.py

Codecov / codecov/patch

pygmt/datatypes/image.py#L86-L87

Added lines #L86 - L87 were not covered by tests
a=self.data[: header.n_bands * header.mx * header.my],
seisman marked this conversation as resolved.
Show resolved Hide resolved
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] = {

Check warning on line 93 in pygmt/datatypes/image.py

Codecov / codecov/patch

pygmt/datatypes/image.py#L93

Added line #L93 was not covered by tests
"x": self.x[: header.n_columns],
"y": self.y[: header.n_rows],
"band": np.arange(stop=3, dtype=np.uint8),
weiji14 marked this conversation as resolved.
Show resolved Hide resolved
}

# Create the xarray.DataArray object
image = xr.DataArray(data=data, coords=coords, dims=("y", "x", "band"))

Check warning on line 100 in pygmt/datatypes/image.py

Codecov / codecov/patch

pygmt/datatypes/image.py#L100

Added line #L100 was not covered by tests

return image

Check warning on line 102 in pygmt/datatypes/image.py

Codecov / codecov/patch

pygmt/datatypes/image.py#L102

Added line #L102 was not covered by tests
Loading