Skip to content
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.

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

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

Merged
merged 10 commits into from
Aug 30, 2023
39 changes: 39 additions & 0 deletions examples/gallery/images/rgb_image.py
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

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

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

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
weiji14 marked this conversation as resolved.
Show resolved Hide resolved
image = rioxarray.open_rasterio(
filename="https://oin-hotosm.s3.us-east-1.amazonaws.com/64d6a49a19cb3a000147a65b/0/64d6a49a19cb3a000147a65c.tif",
Copy link
Member

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

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,
)
# Subset to area of Lāhainā in EPSG:32604 coordinates
image = image.rio.clip_box(minx=738000, maxx=755000, miny=2300000, maxy=2318000)
image

###############################################################################
# Plot the RGB imagery
weiji14 marked this conversation as resolved.
Show resolved Hide resolved
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
weiji14 marked this conversation as resolved.
Show resolved Hide resolved
frame=[r"WSne+tL@!a\225hain@!a\225, Hawai\140i on 9 Aug 2023", "af"],
Copy link
Member Author

Choose a reason for hiding this comment

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

Changed to using \140 instead of \047 as it looks a bit closer to the Okina ʻ letter (U+02BB) (see also #1474 (comment)). Now it looks like this:

Lāhainā, Hawaiʻi in Times-Roman font.

If anyone has ideas on how to use the actual Okina letter, let me know! I did try using ʻ directly, but it shows up like this:

Hawaiʻi with the Okina letter plotted incorrectly.

seisman marked this conversation as resolved.
Show resolved Hide resolved
)
fig.savefig("Lāhainā.png")
fig.show()