WIP: Implement virtualfile_from_image

[weiji14]

Description of proposed changes

Enable passing in 3-band images to GMT via a virtualfile mechanism instead of using a temporary GeoTIFF file (thus superseding #2590). This means rioxarray won't need to be installed for plotting 3-band images with grdimage anymore!

TODO:

• Initial implementation
• Deprecate tempfile_from_image
• Add/update unit tests

Notes:

• Implementation based on virtualfile_from_grid originally created in Refactor the gmt.clib package #210
• May decide to wrap this into the work at DRAFT: Passing GMT_GRID as virtual file #3099 (have a general virtualfile_from_xr method that works with any no. of bands?)
• GMT.jl implementation at https://github.com/GenericMappingTools/GMT.jl/blob/v1.19.0/src/gmt_main.jl#L458-L527

Supersedes #2590, xref #1555.

Reminders

• Run make format and make check to make sure the code follows the style guide.
• Add tests for new features or tests that would have caught the bug that you're fixing.
• Add new public functions/methods/classes to doc/api/index.rst.
• Write detailed docstrings for all functions/methods.
• If wrapping a new module, open a 'Wrap new GMT module' issue and submit reasonably-sized PRs.
• If adding new functionality, add an example to docstrings or tutorials.
• Use underscores (not hyphens) in names of Python files and directories.

Slash Commands

You can write slash commands (/command) in the first line of a comment to perform
specific operations. Supported slash command is:

• /format: automatically format and lint the code

[seisman]

• Implementation based on virtualfile_from_grid originally created in Refactor the gmt.clib package #210

Most tests pass in CI, which means the virtualfile_from_image function almost works. Since the new virtualfile_from_image is almost the same as the existing virtualfile_from_grid function, why not just modify the virtualfile_from_grid function to make it support both grids and images?

• May decide to wrap this into the work at DRAFT: Passing GMT_GRID as virtual file #3099 (have a general virtualfile_from_xr method that works with any no. of bands?)

PR #3099 will likely take longer to implement due to upstream bugs. So better to extend the virtualfile_from_grid first.

[weiji14]

• Implementation based on virtualfile_from_grid originally created in Refactor the gmt.clib package #210

Most tests pass in CI, which means the virtualfile_from_image function almost works. Since the new virtualfile_from_image is almost the same as the existing virtualfile_from_grid function, why not just modify the virtualfile_from_grid function to make it support both grids and images?

Actually if you look at the Python 3.12 / NumPy 2.1 tests (where the test_grdimage_image.py tests are being run), you'll see some segfaults. Let's maybe have a separate virtualfile_from_image method first as we improve the implementation so that it is easier to examine things, and decide whether to combine with virtualfile_from_grid later.

• May decide to wrap this into the work at DRAFT: Passing GMT_GRID as virtual file #3099 (have a general virtualfile_from_xr method that works with any no. of bands?)

PR #3099 will likely take longer to implement due to upstream bugs. So better to extend the virtualfile_from_grid first.

👍

[weiji14]

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?

[seisman]

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

[weiji14]

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)

[seisman]

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

[seisman]

• 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