Allow passing RGB xarray.DataArray images into grdimage #2590

weiji14 · 2023-07-03T06:02:07Z

Description of proposed changes

Enables grdimage to plot 3-band RGB images stored in an xarray.DataArray with a shape like (3, Y, X).

Example usage based on https://corteva.github.io/rioxarray/stable/examples/COG.html:

import pygmt
import rioxarray

image = rioxarray.open_rasterio(
    "https://oin-hotosm.s3.amazonaws.com/5d7dad0becaf880008a9bc88/0/5d7dad0becaf880008a9bc89.tif",
    masked=True,
    overview_level=4,
)
assert image.sizes == {"band": 3, "y": 312, "x": 688}

fig = pygmt.Figure()
fig.grdimage(grid=image, frame=True)
fig.show()

produces

Implementation currently saves the 3-band xarray.DataArray image to a temporary GeoTIFF before passing it to grdimage, similar to what tilemap did in #2394

pygmt/pygmt/src/tilemap.py

Lines 151 to 156 in e98efe1

    
           with GMTTempFile(suffix=".tif") as tmpfile: 
        
               raster.rio.to_raster(raster_path=tmpfile.name) 
        
               with Session() as lib: 
        
                   lib.call_module( 
        
                       module="grdimage", args=build_arg_string(kwargs, infile=tmpfile.name) 
        
                   )

Notes:

Conversion of xarray.DataArray to temporary GeoTIFF file requires rioxarray.
Only certain dtypes like uint8 and float32 are supported for plotting, so need to either warn users if other dtypes are passed in. Best scenario is if histogram equalization via grdhisteq Wrap grdhisteq #1433 or normalization to 0-255 via grdmix (Wrap grdmix #1437) is applied
This implementation will need to be refactored when GMT_IMAGE is wrapped in PyGMT so that temporary files can be avoided, see Integrate with rioxarray and implement GMT_IMAGE to read RGB GeoTIFFs #1555. The unit tests in this PR should still be valid though.

TODO:

Initial implementation using tempfile
Raise error on missing rioxarray dependency
Raise warning when unsupported dtypes are passed in
Maybe implement same logic for grdview too?

Closes #370, partially addresses #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 commands are:

/format: automatically format and lint the code
/test-gmt-dev: run full tests on the latest GMT development version

Saving 3-band xarray.DataArray images to a temporary GeoTIFF so that they can be plotted with grdimage.

github-actions · 2023-07-03T06:03:22Z

Summary of changed images

This is an auto-generated report of images that have changed on the DVC remote

Status	Path
added	pygmt/tests/baseline/test_grdimage_image.png

Image diff(s)

Added images

test_grdimage_image.png

Modified images

Path	Old	New

Report last updated at commit 15a7207

weiji14 · 2023-07-04T00:17:40Z

pygmt/tests/test_grdimage_image.py

+    geotiff = which(fname="@earth_day_01d_p", download="c")
+    with rioxarray.open_rasterio(filename=geotiff) as rda:
+        if len(rda.band) == 1:
+            with rasterio.open(fp=geotiff) as src:
+                df_colormap = pd.DataFrame.from_dict(
+                    data=src.colormap(1), orient="index"
+                )
+                array = src.read()
+
+                red = np.vectorize(df_colormap[0].get)(array)
+                green = np.vectorize(df_colormap[1].get)(array)
+                blue = np.vectorize(df_colormap[2].get)(array)
+                # alpha = np.vectorize(df_colormap[3].get)(array)
+
+            rda.data = red
+            da_red = rda.astype(dtype=np.uint8).copy()
+            rda.data = green
+            da_green = rda.astype(dtype=np.uint8).copy()
+            rda.data = blue
+            da_blue = rda.astype(dtype=np.uint8).copy()
+
+            xr_image = xr.concat(objs=[da_red, da_green, da_blue], dim="band")
+        assert xr_image.sizes == {"band": 3, "y": 180, "x": 360}
+        return xr_image


Should replace this with the load_blue_marble function from #2235 once done. Ideally, the rioxarray.open_rasterio function would be able to load this 3-band GeoTIFF from colorinterp directly without the hacky logic here.

weiji14 · 2023-07-04T00:18:56Z

pygmt/src/grdimage.py

+    with GMTTempFile(suffix=".tif") as tmpfile:
+        if hasattr(grid, "dims") and len(grid.dims) == 3:
+            grid.rio.to_raster(raster_path=tmpfile.name)
+            _grid = tmpfile.name
+        else:
+            _grid = grid


Should maybe put this logic in lib.virtualfile_from_data, since we'll need to reuse it for grdview, grdcut, and possibly other GMT modules that work with GMT_IMAGE.

Ok, created a new tempfile_from_image function at 89da916 (similar to tempfile_from_geojson). This can be replaced with lib.virtualfile_from_image once that is implemented.

Putting the temporary GeoTIFF creation logic in a dedicated tempfile_from_image helper function, so that it can be reused by other GMT modules besides grdimage. Also ensure that an ImportError is raised when the `.rio` attribute cannot be found when rioxarray is not installed.

Refactor Figure.tilemap to use the same tempfile_from_image function that generates a temporary GeoTIFF file from the 3-band xarray.DataArray images.

Various updates from upstream GMT at GenericMappingTools/gmt#6258, GenericMappingTools/gmt@9069967, GenericMappingTools/gmt#7260.

Plotting a non-uint8 dtype xarray.DataArray works in grdimage, but the results will likely be incorrect. Warning the user about the incorrect dtype, and suggest recasting to uint8 with 0-255 range, e.g. using a histogram equalization function like skimage.exposure.equalize_hist.

weiji14 · 2023-07-08T07:09:11Z

pygmt/clib/session.py

+            if kind == "image" and data.dtype != "uint8":
+                msg = (
+                    f"Input image has dtype: {data.dtype} which is unsupported, "
+                    "and may result in an incorrect output. Please recast image "
+                    "to a uint8 dtype and/or scale to 0-255 range, e.g. "
+                    "using a histogram equalization function like "
+                    "skimage.exposure.equalize_hist."
+                )
+                warnings.warn(message=msg, category=RuntimeWarning, stacklevel=2)


Took a look at grdhisteq, but it doesn't appear to support 3-band image inputs, only 1-band grids, so suggesting to use skimage.exposure.equalize_hist instead. Could probably raise a feature request to upstream GMT to let grdhisteq support this too?

weiji14 · 2023-07-08T07:09:41Z

pygmt/helpers/utils.py

@@ -112,7 +112,7 @@ def data_kind(data=None, x=None, y=None, z=None, required_z=False):
    Returns
    -------
    kind : str
-        One of: ``'file'``, ``'grid'``, ``'matrix'``, ``'vectors'``.
+        One of: ``'file'``, ``'grid'``, ``image``, ``'matrix'``, ``'vectors'``.


TODO handle merge conflicts after #2493 is merged.

pygmt/helpers/utils.py

Co-authored-by: Dongdong Tian <[email protected]>

pygmt/helpers/tempfile.py

seisman

This PR looks good to me. Are there any more changes you want to make into this PR?

Co-authored-by: Dongdong Tian <[email protected]>

weiji14 · 2023-07-20T21:24:21Z

This PR looks good to me. Are there any more changes you want to make into this PR?

Was planning to wait for #2493 (comment) to be finalized so I can merge those changes on top of here, but I can merge this PR first too. Wanted to add some extra tests for grdview, which I'll try to do this weekend.

Edit: Hmm, and there seems to be a failing test on the Docs build at https://github.com/GenericMappingTools/pygmt/actions/runs/5615992349/job/15217416041?pr=2590#step:7:725 - AttributeError: 'DataArray' object has no attribute 'rio', need to see what's going on. Edit: fixed in d1d692b

weiji14 · 2023-08-07T20:44:16Z

On second thought, let's do the unit tests for grdview in a separate PR, because this one is already quite big and I'd also like to sync the grdview docstring with upstream GMT which can make the diff bigger. I've solved the merge conflicts with #2493 now @seisman, so maybe take another look before I merge this in.

Partially revert a773929 to fix `AttributeError: 'DataArray' object has no attribute 'rio'` in `pygmt/src/tilemap.py`, because the reproject from EPSG:3857 to OGC:CRS84 step requires it.

Allow passing RGB xarray.DataArray images into grdimage

7e4c411

Saving 3-band xarray.DataArray images to a temporary GeoTIFF so that they can be plotted with grdimage.

weiji14 added the enhancement Improving an existing feature label Jul 3, 2023

weiji14 added this to the 0.10.0 milestone Jul 3, 2023

weiji14 self-assigned this Jul 3, 2023

weiji14 mentioned this pull request Jul 3, 2023

Allow passing in Red, Green, Blue grids into grdimage #370

Closed

5 tasks

weiji14 commented Jul 4, 2023

View reviewed changes

weiji14 added 3 commits July 4, 2023 17:51

Let tilemap use tempfile_from_image func in virtualfile_from_data

a773929

Refactor Figure.tilemap to use the same tempfile_from_image function that generates a temporary GeoTIFF file from the 3-band xarray.DataArray images.

Lint to remove unneeded imports

8c89da2

weiji14 mentioned this pull request Jul 5, 2023

Add a tutorial explaining the generally accepted input types #1268

Open

weiji14 added 3 commits July 8, 2023 17:22

Merge branch 'main' into grdimage/rgb_dataarray

ba87dff

Update docstring of grdimage with upstream GMT 6.4.0

819cc8b

Various updates from upstream GMT at GenericMappingTools/gmt#6258, GenericMappingTools/gmt@9069967, GenericMappingTools/gmt#7260.

weiji14 commented Jul 8, 2023

View reviewed changes

weiji14 changed the title ~~WIP: Allow passing RGB xarray.DataArray images into grdimage~~ Allow passing RGB xarray.DataArray images into grdimage Jul 8, 2023

weiji14 added the needs review This PR has higher priority and needs review. label Jul 8, 2023

seisman reviewed Jul 9, 2023

View reviewed changes

pygmt/helpers/utils.py Outdated Show resolved Hide resolved

Update pygmt/helpers/utils.py

fabcb09

Co-authored-by: Dongdong Tian <[email protected]>

seisman reviewed Jul 20, 2023

View reviewed changes

pygmt/helpers/tempfile.py Outdated Show resolved Hide resolved

seisman approved these changes Jul 20, 2023

View reviewed changes

Set intersphinx link for xr.DataArray in docstring

ba9a2a3

Co-authored-by: Dongdong Tian <[email protected]>

weiji14 mentioned this pull request Aug 4, 2023

Better handling of optional virtual files (e.g., shading in Figure.grdimage) #2493

Merged

7 tasks

Merge branch 'main' into grdimage/rgb_dataarray

40eb80a

weiji14 marked this pull request as ready for review August 7, 2023 20:42

weiji14 added final review call This PR requires final review and approval from a second reviewer and removed needs review This PR has higher priority and needs review. labels Aug 7, 2023

Add back rioxarray import and ImportError check to tilemap.py

d1d692b

Partially revert a773929 to fix `AttributeError: 'DataArray' object has no attribute 'rio'` in `pygmt/src/tilemap.py`, because the reproject from EPSG:3857 to OGC:CRS84 step requires it.

Lint to add missing line

15a7207

seisman approved these changes Aug 8, 2023

View reviewed changes

weiji14 merged commit cf022a2 into main Aug 8, 2023
8 of 14 checks passed

weiji14 deleted the grdimage/rgb_dataarray branch August 8, 2023 02:59

weiji14 removed the final review call This PR requires final review and approval from a second reviewer label Aug 8, 2023

weiji14 mentioned this pull request Aug 25, 2023

Add gallery example for plotting an RGB image from an xarray.DataArray #2641

Merged

7 tasks

weiji14 mentioned this pull request Jan 7, 2024

Fix fixture_xr_image to open earth_day_01d_p directly with rioxarray #2963

Merged

7 tasks

weiji14 mentioned this pull request Sep 29, 2024

WIP: Implement virtualfile_from_image #3468

Draft

10 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Allow passing RGB xarray.DataArray images into grdimage #2590

Allow passing RGB xarray.DataArray images into grdimage #2590

weiji14 commented Jul 3, 2023 •

edited

Loading

github-actions bot commented Jul 3, 2023 •

edited

Loading

Added images

Modified images

weiji14 Jul 4, 2023

weiji14 Jul 4, 2023

weiji14 Jul 4, 2023

weiji14 Jul 8, 2023

weiji14 Jul 8, 2023

seisman left a comment

weiji14 commented Jul 20, 2023 •

edited

Loading

weiji14 commented Aug 7, 2023 •

edited

Loading

	with GMTTempFile(suffix=".tif") as tmpfile:
	raster.rio.to_raster(raster_path=tmpfile.name)
	with Session() as lib:
	lib.call_module(
	module="grdimage", args=build_arg_string(kwargs, infile=tmpfile.name)
	)

Allow passing RGB xarray.DataArray images into grdimage #2590

Allow passing RGB xarray.DataArray images into grdimage #2590

Conversation

weiji14 commented Jul 3, 2023 • edited Loading

github-actions bot commented Jul 3, 2023 • edited Loading

Summary of changed images

Image diff(s)

Added images

Modified images

weiji14 Jul 4, 2023

Choose a reason for hiding this comment

weiji14 Jul 4, 2023

Choose a reason for hiding this comment

weiji14 Jul 4, 2023

Choose a reason for hiding this comment

weiji14 Jul 8, 2023

Choose a reason for hiding this comment

weiji14 Jul 8, 2023

Choose a reason for hiding this comment

seisman left a comment

Choose a reason for hiding this comment

weiji14 commented Jul 20, 2023 • edited Loading

weiji14 commented Aug 7, 2023 • edited Loading

weiji14 commented Jul 3, 2023 •

edited

Loading

github-actions bot commented Jul 3, 2023 •

edited

Loading

weiji14 commented Jul 20, 2023 •

edited

Loading

weiji14 commented Aug 7, 2023 •

edited

Loading