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)>
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 1 2 3
  * x            (x) float64 7.38e+05 7.38e+05 7.38e+05 ... 7.55e+05 7.55e+05
  * y            (y) float64 2.318e+06 2.318e+06 2.318e+06 ... 2.3e+06 2.3e+06
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:  Area
    scale_factor:   1.0
    add_offset:     0.0

xarray.DataArray

• band: 3
• y: 922
• x: 871

• 50 53 55 54 51 46 43 47 47 49 50 ... 63 57 49 60 50 54 59 55 52 58 58

  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: (4)
  • band
    (band)
    int64
    1 2 3

    array([1, 2, 3])

  • x
    (x)
    float64
    7.38e+05 7.38e+05 ... 7.55e+05

    axis :

    X

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    units :

    metre

    array([737998.046875, 738017.578125, 738037.109375, ..., 754951.171875,
           754970.703125, 754990.234375])

  • y
    (y)
    float64
    2.318e+06 2.318e+06 ... 2.3e+06

    axis :

    Y

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    units :

    metre

    array([2317998.046875, 2317978.515625, 2317958.984375, ..., 2300048.828125,
           2300029.296875, 2300009.765625])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 4N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-159],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32604"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 4N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -159.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 4N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-159],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32604"]]

    GeoTransform :

    737988.28125 19.53125 0.0 2318007.8125 0.0 -19.53125

    array(0)

• Indexes: (3)
  • band
    PandasIndex

    PandasIndex(Index([1, 2, 3], dtype='int64', name='band'))

  • x
    PandasIndex

    PandasIndex(Index([737998.046875, 738017.578125, 738037.109375, 738056.640625,
           738076.171875, 738095.703125, 738115.234375, 738134.765625,
           738154.296875, 738173.828125,
           ...
           754814.453125, 754833.984375, 754853.515625, 754873.046875,
           754892.578125, 754912.109375, 754931.640625, 754951.171875,
           754970.703125, 754990.234375],
          dtype='float64', name='x', length=871))

  • y
    PandasIndex

    PandasIndex(Index([2317998.046875, 2317978.515625, 2317958.984375, 2317939.453125,
           2317919.921875, 2317900.390625, 2317880.859375, 2317861.328125,
           2317841.796875, 2317822.265625,
           ...
           2300185.546875, 2300166.015625, 2300146.484375, 2300126.953125,
           2300107.421875, 2300087.890625, 2300068.359375, 2300048.828125,
           2300029.296875, 2300009.765625],
          dtype='float64', name='y', length=922))

• Attributes: (3)

  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()