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.

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.

pygmt/clib/session.py

weiji14 · 2024-06-18T09:01:46Z

pygmt/clib/session.py

+            with GMTTempFile(suffix=".tif") as tmpfile:
+                self.call_module("write", f"{vfname} {tmpfile.name} -Ti")
+                return xr.load_dataarray(tmpfile.name)


For 3-band images, we might need to use rioxarray.open_rasterio instead of xarray.load_dataarray? I've been thinking of modifying pygmt.io.load_dataarray to support both 1-band grids and 3-band images (see commit aa8a474 and fc060ea in #2235).

I'm a little confused. After #3128 is finished, we should be able to use self.read_virtualfile(vfname, kind=kind).contents.to_dataarray() to convert the GMT_IMAGE structure into a dataarray, instead of replying on temporary files.

In other words, we should merge changes of #3128 into this branch, right?

Ah I see, so this is just a placeholder for now. Yes, this PR will need to wait for #3128, though I'm thinking of cherry-picking some of the clib-related changes here in #3115 to #3128, or opening a separate PR for that.

The original idea of this PR is to make grdcut support both grids (via GMT_GRID) and images (via temporary files because GMT_IMAGE is not wrapped yet). I can continue this PR without waiting for #3128. Maybe it's easier.

That could work too. Finish up this PR #3115 with the placeholder code, and then we'll refactor things in #3128 to remove the tempfile path. Either way, we should get both of these more-or-less ready at the same time, so that we have both of them for PyGMT v0.13.0.

I've decided to finish this PR without waiting for PR #3128. Two tests are also added which can help test changes in PR #3128, too.

pygmt/datatypes/image.py

pygmt/src/grdcut.py

Co-authored-by: Wei Ji <[email protected]>

seisman · 2024-06-25T03:59:10Z

pygmt/src/grdcut.py

+        case "file":
+            realpath = which(grid, download="a")
+            if isinstance(realpath, list):
+                realpath = realpath[0]
+            outkind = "image" if realpath.endswith(".tif") else "grid"


When a file is given, there are two different ways to determine if the file is a grid or an image:

Based on file extensions

Read the file (header only if possible) via xarray or rasterio and check the data dims.

I guess option 2 is more robust and preferred.

Currently use option 1 because it's easy for coding.

Co-authored-by: Wei Ji <[email protected]>

pygmt/clib/session.py

seisman · 2024-07-08T03:06:05Z

pygmt/tests/test_grdcut_image.py

+            dtype=np.uint8,
+        ),
+        coords={
+            "band": [1, 2, 3],


An off-topic question:

In [3]: import rioxarray In [4]: rioxarray.open_rasterio("earth_day_01d_p.tif") Out[4]: <xarray.DataArray (band: 3, y: 180, x: 360)> Size: 194kB [194400 values with dtype=uint8] Coordinates: * band (band) int64 24B 1 2 3 * x (x) float64 3kB -179.5 -178.5 -177.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 spatial_ref int64 8B 0 Attributes: AREA_OR_POINT: Area STATISTICS_MAXIMUM: 255 STATISTICS_MEAN: 49.507654320988 STATISTICS_MINIMUM: 10 STATISTICS_STDDEV: 65.299010755209 STATISTICS_VALID_PERCENT: 100 _FillValue: 0 scale_factor: 1.0 add_offset: 0.0

The band here is [1, 2, 3], but tile_map, band is set to [0, 1, 2]:

pygmt/pygmt/datasets/tile_map.py

Line 160 in 41a0fe0

"band": np.array(object=[0, 1, 2], dtype=np.uint8), # Red, Green, Blue

Are there any specific reasons for the differences or any conventions about band?

Oh right, the 123 numbering comes from GDAL/rasterio. E.g. if you do gdalinfo ~/.gmt/server/earth/earth_day/earth_day_01d_p.tif, you'll get:

Driver: GTiff/GeoTIFF Files: /home/user/.gmt/server/earth/earth_day/earth_day_01d_p.tif Size is 360, 180 ... Band 1 Block=360x7 Type=Byte, ColorInterp=Undefined Min=10.000 Max=255.000 Minimum=10.000, Maximum=255.000, Mean=49.508, StdDev=65.299 NoData Value=0 Metadata: STATISTICS_MAXIMUM=255 STATISTICS_MEAN=49.507654320988 STATISTICS_MINIMUM=10 STATISTICS_STDDEV=65.299010755209 STATISTICS_VALID_PERCENT=100 Band 2 Block=360x7 Type=Byte, ColorInterp=Undefined Min=10.000 Max=255.000 Minimum=10.000, Maximum=255.000, Mean=49.911, StdDev=62.129 NoData Value=0 Metadata: STATISTICS_MAXIMUM=255 STATISTICS_MEAN=49.910941358025 STATISTICS_MINIMUM=10 STATISTICS_STDDEV=62.129078727248 STATISTICS_VALID_PERCENT=100 Band 3 Block=360x7 Type=Byte, ColorInterp=Undefined Min=10.000 Max=255.000 Minimum=10.000, Maximum=255.000, Mean=65.605, StdDev=44.820 NoData Value=0 Metadata: STATISTICS_MAXIMUM=255 STATISTICS_MEAN=65.605154320988 STATISTICS_MINIMUM=10 STATISTICS_STDDEV=44.819618148349 STATISTICS_VALID_PERCENT=100

The band here is [1, 2, 3], but tile_map, band is set to [0, 1, 2]:

Yeah, I simply forgot about GDAL's start from 1 convention when implementing #2125 😅. I'm guessing there's a historical reason for this, and I'm not sure if this is the exact reason why GDAL counts from 1, but optical satellite bands (e.g. from Landsat) are usually numbered Band 1, Band 2, etc, so that could be why? Unsure if we should re-number things, or just leave it.

I think we should follow the GDAL convention and let band indexing start from 1.

seisman · 2024-07-08T03:07:56Z

pygmt/tests/test_grdcut_image.py

+    )
+
+
+@pytest.mark.skipif(not _HAS_RIOXARRAY, reason="rioxarray is not installed")


The pytest.mark.skipif marker can be removed after #3128 is implemented.

seisman · 2024-07-08T03:08:40Z

pygmt/tests/test_grdcut_image.py

+    path = which("@earth_day_01d", download="a")
+    raster = rioxarray.open_rasterio(path)


These lines can be replaced by load_blue_marble after #2235.

pygmt/tests/test_grdcut_image.py

Co-authored-by: Wei Ji <[email protected]>

seisman mentioned this pull request Mar 18, 2024

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

Merged

5 tasks

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

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 reviewed Jun 20, 2024

View reviewed changes

seisman and others added 5 commits June 20, 2024 21:52

Apply suggestions from code review

3a24ebd

Co-authored-by: Wei Ji <[email protected]>

Merge branch 'main' into gmtimage

4eee7e6

Address reviewer's comments

5e390d4

Fix GMT_OUT

003383d

Merge branch 'main' into gmtimage

606ac7e

seisman commented Jun 25, 2024

View reviewed changes

seisman added 3 commits July 7, 2024 15:57

Merge branch 'main' into gmtimage

c6cdcc8

Revert changes for _GMT_IMAGE

377941a

Use rioxarray.open_rasterio for loading images

20617f5

Update pygmt/clib/session.py

1c8312c

Co-authored-by: Wei Ji <[email protected]>

seisman commented Jul 8, 2024

View reviewed changes

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

seisman added 6 commits July 8, 2024 09:20

Use rioxarray.open_rasterio in a context manager

3913430

Merge branch 'main' into gmtimage

812a225

Merge remote-tracking branch 'origin/gmtimage' into gmtimage

90bd29e

Fix mypy errors

ab77187

Move grdcut image tests to a separate test file

6f3e474

Fix copy & paste errors

5b07dd9

seisman commented Jul 8, 2024

View reviewed changes

weiji14 reviewed Jul 8, 2024

View reviewed changes

pygmt/tests/test_grdcut_image.py Show resolved Hide resolved

Run codspeed benchmark for test_grdcut_image_dataarray

31272ab

Co-authored-by: Wei Ji <[email protected]>

seisman marked this pull request as draft July 8, 2024 05:36

seisman removed the needs review This PR has higher priority and needs review. label Jul 8, 2024

This was referenced Jul 8, 2024

POC: Wrap GMT_Read_Data and read datasets/grids/images into GMT data container #3318

Closed

load_tile_map: Fix the raster band indexing, should start from 1 #3322

Merged

This was referenced Jul 19, 2024

GMT_IMAGE: Implement the GMT_IMAGE.to_dataarray method #3128

Draft

Add libgdal-hdf5 dependency conda-forge/gmt-feedstock#293

Draft

seisman removed this from the 0.13.0 milestone Aug 4, 2024

seisman added 3 commits August 5, 2024 16:47

Merge branch 'main' into gmtimage

5a09329

Add the raster_kind function to determine the raster kind

279595b

Simplify the grdcut function

7def4b5

seisman added this to the 0.14.0 milestone Aug 5, 2024

seisman removed this from the 0.14.0 milestone Sep 5, 2024

seisman added 3 commits September 19, 2024 22:24

Merge branch 'main' into gmtimage

be175d8

Merge branch 'datatypes/gmtimage' into gmtimage

2825eae

Remove the workaround for images

bf9275c

seisman changed the base branch from main to datatypes/gmtimage September 20, 2024 08:43

seisman changed the title ~~pygmt.grdcut: Refactor to store output in virtualfiles for grids [Still use temporary files for images]~~ pygmt.grdcut: Refactor to store output in virtualfiles for grids Sep 20, 2024

Merge branch 'datatypes/gmtimage' into gmtimage

e9cb0a5

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 Jun 18, 2024

seisman Jun 20, 2024

weiji14 Jun 20, 2024

seisman Jun 21, 2024

weiji14 Jun 21, 2024

seisman Jul 7, 2024

seisman Jun 25, 2024

seisman Jul 7, 2024

seisman Jul 8, 2024

weiji14 Jul 8, 2024

seisman Jul 8, 2024

seisman Jul 8, 2024

seisman Jul 8, 2024

		)


		@pytest.mark.skipif(not _HAS_RIOXARRAY, reason="rioxarray is not installed")

		path = which("@earth_day_01d", download="a")
		raster = rioxarray.open_rasterio(path)

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

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

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 commented Mar 17, 2024 •

edited

Loading