Plotting bands overlaid on composites#

Sometimes, we may want to overlay a part of a band (for example, the thermal band) on top of a composite. The overlay band can be partially transparent, cropped to a range of values, or both. As an example, let’s make a plot of the December 2015 eruption of Momotombo volcano, Nicaragua. We’ll overlay the thermal band (only pixels above 320 K) on top of an RGB composite to show the ongoing lava flow.

First, we’ll import the required packages and load the sample scene:

import xlandsat as xls
import matplotlib.pyplot as plt
import xarray as xr
import numpy as np

path = xls.datasets.fetch_momotombo()
scene = xls.load_scene(path)
# Fill the missing values due to the volcanic clouds to make it look nicer
scene = xls.interpolate_missing(scene)
scene

<xarray.Dataset>
Dimensions:   (easting: 467, northing: 333)
Coordinates:
  * easting   (easting) float64 5.44e+05 5.44e+05 ... 5.58e+05 5.58e+05
  * northing  (northing) float64 1.369e+06 1.369e+06 ... 1.379e+06 1.379e+06
Data variables:
    blue      (northing, easting) float16 -0.05554 -0.05481 ... 0.07788 0.07727
    green     (northing, easting) float16 0.04834 0.05066 ... 0.1056 0.1027
    red       (northing, easting) float16 0.006226 0.01105 ... 0.1163 0.1045
    nir       (northing, easting) float16 0.01483 0.02124 ... 0.2435 0.2367
    swir1     (northing, easting) float16 0.06189 0.06409 ... 0.2098 0.1976
    swir2     (northing, easting) float16 0.06519 0.06628 ... 0.1571 0.1527
    thermal   (northing, easting) float16 293.5 293.5 293.8 ... 311.0 312.0
Attributes: (12/19)
    Conventions:                CF-1.8
    title:                      Landsat 8 scene from 2015-12-05 (path/row=17/51)
    digital_object_identifier:  https://doi.org/10.5066/P9OGBGM6
    origin:                     Image courtesy of the U.S. Geological Survey
    landsat_product_id:         LC08_L2SP_017051_20151205_20200908_02_T1
    processing_level:           L2SP
    ...                         ...
    ellipsoid:                  WGS84
    date_acquired:              2015-12-05
    scene_center_time:          16:06:06.8773380Z
    wrs_path:                   17
    wrs_row:                    51
    mtl_file:                   GROUP = LANDSAT_METADATA_FILE\n  GROUP = PROD...

Now we can plot an RGB composite and thermal band separately to see that they have to show:

# Make the composite
rgb = xls.composite(scene, rescale_to=(0, 0.6))

# Histogram equalization for a better looking image
rgb = xls.equalize_histogram(rgb, clip_limit=0.02, kernel_size=200)

# Plot the RGB and thermal separately
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 12))

rgb.plot.imshow(ax=ax1)
scene.thermal.plot.imshow(ax=ax2, cmap="magma")

ax1.set_aspect("equal")
ax2.set_aspect("equal")

plt.show()

Notice that the lava flow is clearly visible as temperatures above 320 K in the thermal band but it’s difficult to see where the volcano and other landmarks are. Looking at the RGB composite, we can’t really make out the lava flow but we have a clear picture of where the volcano is and where the old lava flows are. A way to show the thermal data with the geographic context of the RGB is to overlay the two in a single plot.

To do so, we’ll first create a version of the thermal band that has all pixels with temperature below 320 K set to NaN (not-a-number). This is used to indicate to matplotlib that a pixel should be transparent. An easy way to do this is with the xarray.where function:

# If the condition is true, use the thermal values. If it's false, use nan
lava = xr.where(scene.thermal >= 320, scene.thermal, np.nan, keep_attrs=True)

fig, ax = plt.subplots(1, 1, figsize=(10, 6))

lava.plot.imshow(ax=ax, cmap="magma")

ax.set_aspect("equal")
plt.show()

Note

We used the keep_attrs=True parameter to tell xarray that it should keep the metadata from the original band in the lava-only version. This will preserve the information on units, procedence, etc. But be careful with this since it can lead to metadata being propagated when it’s no longer valid.

Now that we have an xarray.DataArray with the lava flow only, we can plot that on top of the RGB composite and add a bit of transparency using the alpha parameter of imshow.

fig, ax = plt.subplots(1, 1, figsize=(10, 6))

# RGB goes first so it's at the bottom
rgb.plot.imshow(ax=ax)
lava.plot.imshow(ax=ax, cmap="magma", alpha=0.6)

ax.set_aspect("equal")
plt.show()

With the plot above, all of the information we have available about the lava flow is displayed in a nice format.

xlandsat v0.4.0

Plotting bands overlaid on composites

Plotting bands overlaid on composites#