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.

Sign up for GitHub

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

WIP: Implement virtualfile_from_image #3468

Draft
wants to merge 5 commits into
base: main
Choose a base branch
from
Draft

WIP: Implement virtualfile_from_image #3468

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
331f9aa
Initial implementation of virtualfile_from_image
weiji14 Sep 29, 2024
8c63ae6
Update pygmt/clib/session.py
weiji14 Sep 30, 2024
5cc4cc7
Merge branch 'main' into vfile/from_image
weiji14 Sep 30, 2024
7828153
Typo
weiji14 Sep 30, 2024
149d2ae
Merge branch 'main' into vfile/from_image
weiji14 Nov 6, 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
5 changes: 3 additions & 2 deletions pygmt/clib/conversion.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,9 +84,10 @@ def dataarray_to_matrix(grid):
>>> print(inc)
[2.0, 2.0]
"""
if len(grid.dims) != 2:
if len(grid.dims) not in {2, 3}:
raise GMTInvalidInput(
f"Invalid number of grid dimensions '{len(grid.dims)}'. Must be 2."
f"Invalid number of grid/image dimensions '{len(grid.dims)}'. "
"Must be 2 for grid, or 3 for image."
)
# Extract region and inc from the grid
region = []
Expand Down
77 changes: 67 additions & 10 deletions pygmt/clib/session.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,12 +29,7 @@
from pygmt.clib.loading import get_gmt_version, load_libgmt
from pygmt.datatypes import _GMT_DATASET, _GMT_GRID, _GMT_IMAGE
from pygmt.exceptions import GMTCLibError, GMTCLibNoSessionError, GMTInvalidInput
from pygmt.helpers import (
_validate_data_input,
data_kind,
tempfile_from_geojson,
tempfile_from_image,
)
from pygmt.helpers import _validate_data_input, data_kind, tempfile_from_geojson

FAMILIES = [
"GMT_IS_DATASET", # Entity is a data table
Expand Down Expand Up @@ -868,6 +863,10 @@
... gmttype = ses._check_dtype_and_dim(data, ndim=2)
... gmttype == ses["GMT_FLOAT"]
True
>>> data = np.ones((5, 3, 2), dtype="uint8")
>>> with Session() as ses:
... gmttype = ses._check_dtype_and_dim(data, ndim=3)
... gmttype == ses["GMT_UCHAR"]
weiji14 marked this conversation as resolved.
Show resolved Hide resolved
"""
# Check that the array has the given number of dimensions
if array.ndim != ndim:
Expand Down Expand Up @@ -1006,9 +1005,9 @@
f"Failed to put strings of type {strings.dtype} into dataset"
)

def put_matrix(self, dataset, matrix, pad=0):
def put_matrix(self, dataset, matrix, pad=0, ndim=2):
"""
Attach a numpy 2-D array to a GMT dataset.
Attach a numpy n-D (2-D or 3-D) array to a GMT dataset.

Use this function to attach numpy array data to a GMT dataset and pass
it to GMT modules. Wraps ``GMT_Put_Matrix``.
Expand Down Expand Up @@ -1048,7 +1047,7 @@
restype=ctp.c_int,
)

gmt_type = self._check_dtype_and_dim(matrix, ndim=2)
gmt_type = self._check_dtype_and_dim(matrix, ndim=ndim)
matrix_pointer = matrix.ctypes.data_as(ctp.c_void_p)
status = c_put_matrix(
self.session_pointer, dataset, gmt_type, pad, matrix_pointer
Expand Down Expand Up @@ -1610,6 +1609,64 @@
with self.open_virtualfile(*args) as vfile:
yield vfile

@contextlib.contextmanager
def virtualfile_from_image(self, image: xr.DataArray):
"""
Store a image in a virtual file.

Use the virtual file name to pass in the data in your image to a GMT module.
Images must be :class:`xarray.DataArray` instances.

Context manager (use in a ``with`` block). Yields the virtual file name that you
can pass as an argument to a GMT module call. Closes the virtual file upon exit
of the ``with`` block.

The virtual file will contain the image as a ``GMT_MATRIX`` with extra metadata.

Use this instead of creating a data container and virtual file by hand with
:meth:`pygmt.clib.Session.create_data`, :meth:`pygmt.clib.Session.put_matrix`,
and :meth:`pygmt.clib.Session.open_virtualfile`.

The image data matrix must be C contiguous in memory. If it is not (e.g., it is
a slice of a larger array), the array will be copied to make sure it is.

Parameters
----------
image : :class:`xarray.DataArray`
The image that will be included in the virtual file.

Yields
------
fname : str
The name of virtual file. Pass this as a file name argument to a GMT module.

"""
_gtype = {0: "GMT_GRID_IS_CARTESIAN", 1: "GMT_GRID_IS_GEO"}[image.gmt.gtype]
_reg = {0: "GMT_GRID_NODE_REG", 1: "GMT_GRID_PIXEL_REG"}[image.gmt.registration]

Check warning on line 1645 in pygmt/clib/session.py

View check run for this annotation

Codecov / codecov/patch

pygmt/clib/session.py#L1644-L1645 

Added lines #L1644 - L1645 were not covered by tests

# Conversion to a C-contiguous array needs to be done here and not in put_matrix
# because we need to maintain a reference to the copy while it is being used by
# the C API. Otherwise, the array would be garbage collected and the memory
# freed. Creating it in this context manager guarantees that the copy will be
# around until the virtual file is closed. The conversion is implicit in
# dataarray_to_matrix.
matrix, region, inc = dataarray_to_matrix(image)

Check warning on line 1653 in pygmt/clib/session.py

View check run for this annotation

Codecov / codecov/patch

pygmt/clib/session.py#L1653 

Added line #L1653 was not covered by tests
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A little unsure about the dimension order when flattening the 3-D image into a 2-D matrix (hard to decipher from https://docs.generic-mapping-tools.org/6.5/devdocs/api.html#gmt-images). Does GMT expect band sequential (BSQ) order, or something else?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe this issue (GenericMappingTools/gmt#3299) is helpful.

Copy link
Member Author

@weiji14 weiji14 Oct 1, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah yes, that is very helpful, so Paul said this:

  • The exact arrangement of r,g,b,a values in the image array I->data depends on the chosen setting via API_IMAGE_LAYOUT. The default in GMT is TRPa [RGBARGBARGBA...]. Modules can thus set the format most suitable for their use.

  • We support (T|B)(R|C)(B|P). First letter says image is top down or bottom up, second letter says we go by rows or column, and third letter says the r,g,b is interleaved by bands (all R first, followed by all G, then B, then optionally A) or pixel (rgb).

The 'P' in 'TRP' sounds like band interleaved by pixel (BIP) if I'm interpreting this correctly. So something like:

1. (R->G->B)->(R->G->B)->(R->G->B)
2. (R->G->B)->(R->G->B)->(R->G->B)
3. (R->G->B)->(R->G->B)->(R->G->B)

where we go first along row 1 from left-to-right with RGB values of each pixel, and then row 2, row 3, etc.

That said, the format appears to be determined by API_IMAGE_LAYOUT, and I tried to get it using the following code:

import pygmt

with pygmt.clib.Session() as lib:
    print(lib.info["image layout"])  # API_IMAGE_LAYOUT

but it comes up empty even though I'm using GMT 6.5.0? There's a note here that mentions API_IMAGE_LAYOUT is not defined if GMT is not compiled with GDAL:

pygmt/pygmt/clib/session.py

Lines 202 to 207 in f8db417

# For GMT<6.4.0, API_IMAGE_LAYOUT is not defined if GMT is not
# compiled with GDAL. Since GMT 6.4.0, GDAL is a required GMT
# dependency. The code block can be refactored after we bump
# the minimum required GMT version to 6.4.0.
with contextlib.suppress(GMTCLibError):
self._info["image layout"] = self.get_default("API_IMAGE_LAYOUT")

Unsure what is going on, maybe I should test this out on the branch in #3450. Edit: See #3450 (comment)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

From the GMT source codes, it seems the value is set after reading an image via GDAL. More specifically, these lines:

https://github.com/GenericMappingTools/gmt/blob/54dbea3217e45c982d83721ebd2740c84870c637/src/gmt_gdalread.c#L822.

But I still can't get its value even after reading an image. Need more try.

As for reading an image, the layout information can be obtained from header.mem_layout, see

pygmt/pygmt/datatypes/image.py

Line 62 in f8db417

b'BRPa' 0.5

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  • The exact arrangement of r,g,b,a values in the image array I->data depends on the chosen setting via API_IMAGE_LAYOUT. The default in GMT is TRPa [RGBARGBARGBA...].

I'm unsure if this is correct. From the source codes, it seems the default image layout is TRBa.

https://github.com/GenericMappingTools/gmt/blob/426abad288ef4b87c2cecf4c9633fd849866e45d/src/gmt_constants.h#L322


family = "GMT_IS_IMAGE|GMT_VIA_MATRIX"
geometry = "GMT_IS_SURFACE"
gmt_image = self.create_data(

Check warning on line 1657 in pygmt/clib/session.py

View check run for this annotation

Codecov / codecov/patch

pygmt/clib/session.py#L1655-L1657 

Added lines #L1655 - L1657 were not covered by tests
family,
geometry,
mode=f"GMT_CONTAINER_ONLY|{_gtype}",
ranges=region[0:4], # (xmin, xmax, ymin, ymax) only, leave out (zmin, zmax)
inc=inc[0:2], # (x-inc, y-inc) only, leave out z-inc
registration=_reg,
)
self.put_matrix(gmt_image, matrix, ndim=3)
args = (family, geometry, "GMT_IN|GMT_IS_REFERENCE", gmt_image)
with self.open_virtualfile(*args) as vfile:
yield vfile

Check warning on line 1668 in pygmt/clib/session.py

View check run for this annotation

Codecov / codecov/patch

pygmt/clib/session.py#L1665-L1668 

Added lines #L1665 - L1668 were not covered by tests

@contextlib.contextmanager
def virtualfile_from_stringio(self, stringio: io.StringIO):
r"""
Expand Down Expand Up @@ -1796,7 +1853,7 @@
"arg": contextlib.nullcontext,
"geojson": tempfile_from_geojson,
"grid": self.virtualfile_from_grid,
"image": tempfile_from_image,
"image": self.virtualfile_from_image,
"stringio": self.virtualfile_from_stringio,
# Note: virtualfile_from_matrix is not used because a matrix can be
# converted to vectors instead, and using vectors allows for better
Expand Down
Loading