pygmt.grdcut: Refactor to store output in virtualfiles for grids #3115

seisman · 2024-03-17T15:24:17Z

This PR tries to get rid of temporary files from pygmt.grdcut (#2730).

The GMT_GRID wrapper (PR #2398) can store output 2-D grid in virtual files, but doesn't work for output images. PR #3128 will wrap GMT_IMAGE to store output images in virtual files, but it's still WIP.

This PR makes necessary changes to make pygmt.grdcut work for both grids and images. For images, we still use temporary files but it should be refactored again after PR #3128.

Minimal xarray.DataArray output with data and coordinates, no metadata yet.

pygmt/src/grdcut.py

pygmt/clib/session.py

pygmt/src/grdcut.py

Special case `earth_day` to be loaded as GMT_IMAGE rather than GMT_GRID. Need to use `GMT_OUT|GMT_IS_REFERENCE` in virtualfile_out to avoid segfault, xref #3115.

Extra metadata from the _GMT_GRID_HEADER struct.

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.

Get the registration and gtype info from the grid header and apply it to the GMT accessor attributes.

Trying to match some of the doctests in _GMT_GRID.

…polate

pygmt/clib/session.py

pygmt/src/grdcut.py

pygmt/clib/session.py

pygmt/src/grdcut.py

weiji14 · 2024-09-28T04:53:31Z

pygmt/clib/session.py

+        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


Just double-checking that there's not a function in GMT C already to determine grid/image type? I know there's GMT_Inquire_VirtualFile for virtualfiles, but is there not one for just files? I saw this GMT_Get_Info, but unsure if it'll be usable here.

https://github.com/GenericMappingTools/gmt/blob/7809736ba32d87a4a96b15444419eb176c6a35f3/src/gmt_grdio.c#L3377

gmt_raster_type is the function that GMT uses to determine the raster type, but it's not a public API function.

Oo, magic byte parsing 😆 If that function becomes public someday, maybe we'll use it.

Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com>

weiji14 · 2024-09-29T04:17:45Z

pygmt/clib/session.py

@@ -2311,3 +2312,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) -> Literal["grid", "image"] | None:


Mypy is complaining about the output type hint here:

pygmt/src/grdcut.py:108: error: Incompatible types in assignment (expression has type "Literal['grid', 'image'] | None", variable has type "Literal['grid', 'image']") [assignment]

Maybe we should just use Literal['grid', 'image'] and raise GMTInvalidInput if it is not a grid/image instead of returning None?

Maybe falling back to "grid" and then let GMT itself decide if it can process it or not.

In that case, with a binary image/grid output, maybe just make a function like def _is_image() -> bool?

grdcut can also cut a cube, so the _raster_kind function should return grid/image/cube if possible. Currently, the' cube' support is not yet finished (PR #3150). So, the name _raster_kind makes more sense.

Ah yes, and then there would be 3 try-except statements. Now wishing we could re-use the function from GMT C...

seisman · 2024-09-29T14:28:39Z

pygmt/clib/session.py

+    # 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.


It turns out the logic here is not 100% correct.

An image can be read into a grid container, but only the first band is read.

A grid can be read into an image container, as long as another libgdal package (maybe libgdal-netcdf or libgdal-hdf5) is installed. xref: Add libgdal-hdf5 dependency conda-forge/gmt-feedstock#293

A grid can be read into an image container, as long as another libgdal package (maybe libgdal-netcdf or libgdal-hdf5) is installed. xref:

The first try-except check is to try reading into a GMT_IMAGE struct, and return "image" if there are 3 bands, and "grid" if there is only 1 band. Even if grids can be read into GMT_IMAGE, we will still call it a grid if it has 1-band right?

Trying to refresh my memory about what's happening:

Reading a grid into GMT_GRID: Works!

Reading an image into GMT_GRID: Only the first band is read.

Reading a grid into GMT_IMAGE: It depends on libgdal-hdf5 being installed. Read as an single-band image if installed, otherwise raise errors.

Reading an image into GMT_IMAGE: works.

import pygmt pygmt.grdcut("@earth_relief_01d_p", region=[0, 10, 0, 10])

If libgdal-hdf5 is not installed, _raster_kind first read it as an image, which fails and reports following errors:

ERROR 4: `/home/runner/.gmt/server/earth/earth_relief/earth_relief_01d_p.grd' not recognized as being in a supported file format. It could have been recognized by driver HDF5, but plugin gdal_HDF5.so is not available in your installation. You may install it with 'conda install -c conda-forge libgdal-hdf5' Error: ession [ERROR]: Unable to open earth_relief_01d_p.grd. Error: ession [ERROR]: ERROR reading image with gdalread. pygmt-session (gmtapi_import_image): Could not read from file [earth_relief_01d_p.grd] [Session pygmt-session (43)]: Error returned from GMT API: GMT_IMAGE_READ_ERROR (22) [Session pygmt-session (43)]: Error returned from GMT API: GMT_IMAGE_READ_ERROR (22)

then it tries to read it as a grid, which succeeds. But we still see the error messages above. I believe we need to find a way to suppress these errors or find an alternative way to determine the data kind.

pygmt/tests/test_grdcut_image.py

seisman mentioned this pull request Mar 18, 2024

Wrap GMT's standard data type GMT_GRID for grids #2398

Merged

5 tasks

seisman and others added 2 commits March 20, 2024 11:00

Wrap GMT's standard data type GMT_IMAGE for images

bcf43f0

Initial implementation of to_dataarray method for _GMT_IMAGE class

a052a1a

Minimal xarray.DataArray output with data and coordinates, no metadata yet.

weiji14 reviewed Mar 22, 2024

View reviewed changes

pygmt/src/grdcut.py Outdated Show resolved Hide resolved

pygmt/clib/session.py Outdated Show resolved Hide resolved

pygmt/src/grdcut.py Outdated Show resolved Hide resolved

seisman changed the title ~~WIP: Initial try to support raster images using temporary files~~ WIP: pygmt.grdcut: Support both grids and images Mar 27, 2024

seisman force-pushed the datatypes/gmtgrid branch 3 times, most recently from d8f1160 to d0517e0 Compare April 1, 2024 07:11

seisman mentioned this pull request Apr 1, 2024

Wrap GMT_Inquire_VirtualFile to get the family of virtualfiles #3152

Merged

seisman force-pushed the datatypes/gmtgrid branch from d0517e0 to 4ba67df Compare April 1, 2024 11:32

seisman changed the base branch from datatypes/gmtgrid to main April 1, 2024 11:32

seisman mentioned this pull request Apr 16, 2024

Get rid of temporary files from pygmt functions and plotting methods #2730

Closed

31 tasks

seisman force-pushed the gmtimage branch from 27dfb33 to 84ce139 Compare April 16, 2024 07:23

pygmt.grdcut: Support both grid and image output

59d523c

seisman force-pushed the gmtimage branch from 84ce139 to 59d523c Compare April 16, 2024 14:44

Merge branch 'main' into datatypes/gmtimage

56a6d65

seisman added the enhancement Improving an existing feature label Apr 19, 2024

seisman added 2 commits April 19, 2024 10:22

Merge branch 'main' into gmtimage

3315324

Fix

cea3374

seisman mentioned this pull request Apr 19, 2024

Refactor the _load_remote_dataset function to load tiled and non-tiled grids in a consistent way #3120

Merged

seisman added 2 commits April 20, 2024 01:00

Refactor

80d9837

fix

22fba56

weiji14 added a commit that referenced this pull request Jun 18, 2024

Handle image kind in _load_remote_dataset function

6eabbf4

Special case `earth_day` to be loaded as GMT_IMAGE rather than GMT_GRID. Need to use `GMT_OUT|GMT_IS_REFERENCE` in virtualfile_out to avoid segfault, xref #3115.

weiji14 added 7 commits June 18, 2024 21:16

Merge branch 'main' into datatypes/gmtimage

f71e79c

Small typo fixes and add output type-hint for to_dataarray

4cce4a2

Fix mypy error using np.array([0, 1, 2]) instead of np.arange

e02b650

Parse name and data_attrs from grid/image header

f3d4b1f

Extra metadata from the _GMT_GRID_HEADER struct.

Transpose array to (band, y, x) order and add doctest for to_dataarray

4390136

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.

Set registration and gtype from header

5f25669

Get the registration and gtype info from the grid header and apply it to the GMT accessor attributes.

Print basic shape and padding info in _GMT_IMAGE doctest

a3c6c14

Trying to match some of the doctests in _GMT_GRID.

seisman added 8 commits September 20, 2024 22:25

Fix the incorrect value due to floating number conversion in sphinter…

15b8d53

…polate

Merge branch 'ctypesarray' into datatypes/gmtimage

8433e78

Update the to_dataarray method to match the codes in GMT_GRID

3e3a6f3

image data should has uint8 dtype

12ef40a

Further improve the tests

f64fbb8

Merge branch 'datatypes/gmtimage' into gmtimage

e9cb0a5

Merge branch 'main' into datatypes/gmtimage

4f2ae48

Add a note that currently only 3-band images are supported

d49afed

weiji14 mentioned this pull request Sep 28, 2024

Add function to load Blue Marble dataset #2235

Merged

6 tasks

weiji14 reviewed Sep 28, 2024

View reviewed changes

seisman and others added 5 commits September 28, 2024 15:28

Apply suggestions from code review

a97d0b3

Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com>

Merge branch 'main' into datatypes/gmtimage

f70bec0

Remove the old GMTGridID enums from pygmt/datatypes/header.py

2fd13fb

A minor fix

9972ba1

Merge branch 'datatypes/gmtimage' into gmtimage

ac6b7c3

weiji14 reviewed Sep 29, 2024

View reviewed changes

Base automatically changed from datatypes/gmtimage to main September 29, 2024 08:11

seisman added 5 commits September 29, 2024 19:02

Merge branch 'main' into gmtimage

7c32d41

Let _raster_kind return grid by default

9ec00be

Simplify the grdcut image tests

f3a2f8e

Add one more test for file in & file out

3c12e2b

Fix typos

f852b0d

seisman commented Sep 29, 2024

View reviewed changes

weiji14 reviewed Sep 30, 2024

View reviewed changes

pygmt/tests/test_grdcut_image.py Outdated Show resolved Hide resolved

seisman removed run/benchmark Trigger the benchmark workflow in PRs run/test-gmt-dev Trigger the GMT Dev Tests workflow in PR labels Sep 30, 2024

seisman added 3 commits September 30, 2024 13:05

Merge branch 'main' into gmtimage

5f7683c

Use the new load_blue_marble function

bb1a0b0

Drop the spatial_ref coord

584b5af

seisman self-assigned this Oct 30, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

pygmt.grdcut: Refactor to store output in virtualfiles for grids #3115

pygmt.grdcut: Refactor to store output in virtualfiles for grids #3115

seisman commented Mar 17, 2024 •

edited

Loading

weiji14 Sep 28, 2024

seisman Sep 28, 2024

weiji14 Sep 29, 2024

weiji14 Sep 29, 2024 •

edited

Loading

seisman Sep 29, 2024

weiji14 Sep 29, 2024

seisman Sep 29, 2024

weiji14 Sep 29, 2024

seisman Sep 29, 2024 •

edited

Loading

weiji14 Sep 30, 2024

seisman Sep 30, 2024

pygmt.grdcut: Refactor to store output in virtualfiles for grids #3115

Are you sure you want to change the base?

pygmt.grdcut: Refactor to store output in virtualfiles for grids #3115

Conversation

seisman commented Mar 17, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

weiji14 Sep 29, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

seisman Sep 29, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

seisman commented Mar 17, 2024 •

edited

Loading

weiji14 Sep 29, 2024 •

edited

Loading

seisman Sep 29, 2024 •

edited

Loading