Skip to content

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

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

Merged
merged 10 commits into from
Aug 30, 2023
Merged

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

Changes from 8 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
9d38316
Add gallery example for plotting an RGB image from an xarray.DataArray
weiji14 Aug 25, 2023
60d1be3
Add abbreviation for Cloud-Optimized GeoTIFF
weiji14 Aug 26, 2023
b516c19
Set title font to Times-Roman
weiji14 Aug 26, 2023
e7b66e6
Change apostrophe to Okina symbol ʻ
weiji14 Aug 26, 2023
9085f65
Inline comment on map scale projection
weiji14 Aug 26, 2023
4fd90be
Use ¯ and ` instead of \225 and \140
weiji14 Aug 26, 2023
7586ef7
Merge branch 'main' into gallery/rgb_image
weiji14 Aug 30, 2023
3c0712d
Force loading dataarray into memory
weiji14 Aug 30, 2023
00b427b
Apply suggestions from code review
weiji14 Aug 30, 2023
2553cce
Lint
weiji14 Aug 30, 2023
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
39 changes: 39 additions & 0 deletions examples/gallery/images/rgb_image.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,39 @@
"""
RGB Image
---------
The :meth:`pygmt.Figure.grdimage` method can be used to plot Red, Green, Blue
(RGB) images, or any 3-band false color combination. Here, we'll use
:meth:`rioxarray.open_rasterio` to read a GeoTIFF file into an
Copy link
Member

@yvonnefroehlich yvonnefroehlich Aug 30, 2023

Choose a reason for hiding this comment

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

Is it possible to make this work as a link to https://corteva.github.io/rioxarray/stable/rioxarray.html#rioxarray.open_rasterio?

If this is not easy to do, then we leave it as is. At least in the code example below, the link works.

Copy link
Member Author

@weiji14 weiji14 Aug 30, 2023

Choose a reason for hiding this comment

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

Ah yes, good catch. I used the wrong intersphinx directive, should be :py:func:

Suggested change
:meth:`rioxarray.open_rasterio` to read a GeoTIFF file into an
:py:func:`rioxarray.open_rasterio` to read a GeoTIFF file into an

Copy link
Member

@yvonnefroehlich yvonnefroehlich Aug 31, 2023

Choose a reason for hiding this comment

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

Ah, great. Thanks!

:class:`xarray.DataArray` format, and plot it on a map.

The example below shows a Worldview 2 satellite image over
`Lāhainā, Hawaiʻi during the August 2023 wildfires
<https://en.wikipedia.org/wiki/2023_Hawaii_wildfires#L%C4%81hain%C4%81>`_.
Data is sourced from a Cloud-Optimized GeoTIFF (COG) file hosted on
`OpenAerialMap <https://map.openaerialmap.org>`_ under a
`CC BY-NC 4.0 <https://creativecommons.org/licenses/by-nc/4.0/>`_ license.
"""
import pygmt
import rioxarray

###############################################################################
# Read 3-band data from GeoTIFF into an xarray.DataArray object
with rioxarray.open_rasterio(
filename="https://oin-hotosm.s3.us-east-1.amazonaws.com/64d6a49a19cb3a000147a65b/0/64d6a49a19cb3a000147a65c.tif",
Copy link
Member

@seisman seisman Aug 26, 2023

Choose a reason for hiding this comment

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

This GeoTiff file is almost 1 GB, which is too big to download.

Copy link
Member Author

@weiji14 weiji14 Aug 26, 2023

Choose a reason for hiding this comment

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

The next line overview_level=5 gets a reduced resolution version. Original image has shape (y: 115712, x: 99328), but at overview_level=5, the shape is (y: 1808, x: 1552).

overview_level=5,
) as img:
# Subset to area of Lāhainā in EPSG:32604 coordinates
image = img.rio.clip_box(minx=738000, maxx=755000, miny=2300000, maxy=2318000)
image = image.load() # force loading dataarray into memory
Copy link
Member Author

@weiji14 weiji14 Aug 30, 2023

Choose a reason for hiding this comment

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

Note that the docs build for this example segfaults on both Readthedocs and Ubuntu, see e.g. https://github.com/GenericMappingTools/pygmt/actions/runs/5972759386/job/16203804119?pr=2641#step:7:82. Running this locally, I actually get a segfault too on the first run, but running cd doc && make all a second time works. Not quite sure what the problem is.

Running .load() should fix the segfault? See corteva/rioxarray#550 (comment), and we actually use dataarray.load() in the pygmt.load_dataarray() function in #1439.

Copy link
Member Author

@weiji14 weiji14 Aug 30, 2023

Choose a reason for hiding this comment

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

Yep, it works! See https://pygmt-dev--2641.org.readthedocs.build/en/2641/gallery/images/rgb_image.html

image

###############################################################################
# Plot the RGB imagery
fig = pygmt.Figure()
with pygmt.config(FONT_TITLE="Times-Roman"): # Set title font to Times-Roman
fig.grdimage(
grid=image,
projection="x1:100000", # map scale where 1 cm on map equals 1 km on the ground
frame=[r"WSne+tL@!a¯hain@!a¯, Hawai`i on 9 Aug 2023", "af"],
)
fig.show()