RGB image

The 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 rioxarray.open_rasterio to read a GeoTIFF file into an 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. Data is sourced from a Cloud-Optimized GeoTIFF (COG) file hosted on OpenAerialMap under a CC 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",
    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 the DataArray into memory
image  # noqa: B018

<xarray.DataArray (band: 3, y: 922, x: 871)> Size: 2MB
array([[[ 50,  53,  55, ..., 244, 243, 242],
        [ 49,  52,  53, ..., 244, 243, 242],
        [ 49,  53,  56, ..., 243, 242, 241],
        ...,
        [  0,   0,   0, ...,  27,  32,  30],
        [  0,   0,   0, ...,  23,  28,  28],
        [  0,   0,   0, ...,  20,  24,  23]],

       [[ 57,  60,  62, ..., 240, 239, 238],
        [ 56,  59,  60, ..., 240, 239, 238],
        [ 53,  57,  60, ..., 239, 238, 237],
        ...,
        [  0,   0,   0, ...,  42,  48,  46],
        [  0,   0,   0, ...,  36,  43,  43],
        [  0,   0,   0, ...,  33,  39,  40]],

       [[ 67,  70,  72, ..., 237, 236, 235],
        [ 66,  69,  70, ..., 237, 236, 235],
        [ 65,  69,  72, ..., 236, 235, 234],
        ...,
        [  0,   0,   0, ...,  61,  64,  62],
        [  0,   0,   0, ...,  55,  62,  62],
        [  0,   0,   0, ...,  52,  58,  58]]], dtype=uint8)
Coordinates:
  * band         (band) int64 24B 1 2 3
  * x            (x) float64 7kB 7.38e+05 7.38e+05 ... 7.55e+05 7.55e+05
  * y            (y) float64 7kB 2.318e+06 2.318e+06 ... 2.3e+06 2.3e+06
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:  Area
    scale_factor:   1.0
    add_offset:     0.0

Plot the RGB imagery:

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

Total running time of the script: (0 minutes 1.682 seconds)

Gallery generated by Sphinx-Gallery