Skip to main content

Visualizing Spectra & Plotting

One of the outputs from this tutorial, a pixel's spectral curve!

One of the outputs from this tutorial, a pixels spectral curve!

Try it yourselfRuns free in your browser with no install or setup, and pulls the imagery straight from our Open Data Program.Or view the notebook source on GitHub
Open in Colab

This page details how to produce 2-D scatter plots and spectral curve plots with Wyvern data in Python. These plots can be used to understand the various materials, objects, and species of vegetation that are present in an image.

Getting started​

The notebook is self-contained and installs its own dependencies, so there is nothing to set up in advance.

  • Google Colab (easiest): open the notebook using the "Open in Colab" badge above, then run the cells top to bottom.
  • Local Jupyter: open the notebook in JupyterLab or VS Code (Python 3.10+) and run it.

Either way, the first cell installs the packages that are not already present:

%pip install -q numpy matplotlib requests rasterio

Image we're using​

We're using the image also used in Loading Data in Python. Check out the renderings of the image there, or on the Wyvern Open Data Catalog!

Spectral Plots​

These plots show an individual pixel's surface reflectance across the electromagnetic spectrum that Wyvern's imagery senses. The spectral curve can be used to visually identify the pixel, and can be used by machine learning (and other methods!) to classify the spectral curve as a particular material/species of vegetation. We will dive into classification methods in another tutorial.

Scatter Plots​

2-dimensional scatter plots showcase the distribution of pixels across different bands. This can be used to group spectrally similar pixels, or identify where a particular pixel lies on the plot.

Code​

First things first, let's import all of our packages.

import json
import numpy as np
import rasterio
import matplotlib.pyplot as plt

from matplotlib.patches import Rectangle

Download the Tutorial Files​

Let's download the GeoTIFF and its STAC JSON metadata!

Run the cell below to download these files. Once the download is complete, they'll appear in the same folder as your notebook!

import requests
from pathlib import Path

IMAGE_URL = "https://wyvern-data.com/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a.tiff"
JSON_URL  = "https://wyvern-odp.com/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a.json"

IMAGE_PATH = Path("wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a.tiff")
JSON_PATH  = Path("wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a.json")

downloads = [
    (IMAGE_URL, IMAGE_PATH),
    (JSON_URL, JSON_PATH),
]

for url, path in downloads:
    if path.exists():
        print(f"{path.name} File already here! Skipping download.")
    else:
        print(f"Downloading {path.name}...")
        response = requests.get(url, headers={"User-Agent": "Mozilla/5.0"})
        response.raise_for_status()
        path.write_bytes(response.content)

print("Downloads complete!")

Loading the Data​

We'll start by reading the scale factor from the metadata JSON file we downloaded earlier. Wyvern's L2A products represent surface reflectance, but the values are stored as scaled 16-bit integers to reduce file size.

# Open the JSON metadata file and extract the scale factor
with open(JSON_PATH) as json_file:
    metadata = json.load(json_file)

# The scale factor converts stored uint16 reflectance values into their floating-point representation
scale_factor = metadata["assets"]["Cloud optimized GeoTiff"]["raster:bands"][0]["scale"]

# if you don't have the JSON file, you can hardcode the scale factor:
# scale_factor = 0.0001

print(f"Scale factor: {scale_factor}")

Now let's open the GeoTIFF using rasterio. We'll replace any nodata pixels with np.nan so they are excluded from our calculations, then apply the scale factor to recover the floating-point reflectance values.

with rasterio.open(IMAGE_PATH) as image_file:
    # Convert to float so we can use np.nan and decimal values
    image_arr = image_file.read().astype(np.float32)
    nodata = image_file.nodata
    band_descriptions = image_file.descriptions

# Replace nodata pixels with NaN so they are excluded from calculations
if nodata is not None:
    image_arr[image_arr == nodata] = np.nan

# Apply the scale factor to recover floating-point reflectance values
image_arr *= scale_factor

print(f"Image shape: {image_arr.shape}")
print(f"Image bands: {', '.join(band_descriptions)}")

Building the Wavelength Map​

Just like in Loading Data in Python, we'll build a map from wavelength to band index using the eo:bands metadata, and a small helper to select bands by wavelength.

with open(JSON_PATH) as json_file:
    metadata = json.load(json_file)

wavelength_map = {}
for i, band in enumerate(metadata["assets"]["Cloud optimized GeoTiff"]["eo:bands"]):
    wavelength = round(band["center_wavelength"] * 1000)  # convert from μm to nm
    wavelength_map[wavelength] = i

print(wavelength_map)

Selecting Bands by Wavelength​

def get_band(wavelength_nm: int, return_index: bool = False):
    """
    Returns the image band closest to the requested wavelength.

    Args:
        wavelength_nm: Target wavelength in nanometres.
        return_index: If True, also return the band's index alongside the array.

    Returns:
        A 2D array containing reflectance values for the closest available wavelength band,
        or a (array, band_index) tuple if return_index is True.
    """
    available_wavelengths = list(wavelength_map.keys())
    closest_wavelength = min(available_wavelengths, key=lambda w: abs(w - wavelength_nm))
    band_index = wavelength_map[closest_wavelength]
    print(f"Requested {wavelength_nm}nm -- using {closest_wavelength}nm")

    if return_index:
        return image_arr[band_index], band_index

    return image_arr[band_index]

Spectral Curve Plot

Now, we're ready to plot our spectral curve. We'll select an arbitrary row/column coordinate in the image to select and plot. Feel free to change these values and explore different pixels in the image!

We'll use matplotlib's subplots function to easily plot two things in one plot. This will allow us to show the selected pixel's spectra, along with where that pixel lies in the image! This will help us to provide more context on the spectral curve we are looking at.

# To plot a spectral curve, we first need to pick a pixel in the image to plot!
PIXEL_ROW, PIXEL_COL = 1450, 1800  # Point on image to plot, feel free to change!

pixel_spectra = image_arr[:, PIXEL_ROW, PIXEL_COL]  # Grabbing pixel from image array

# Order our wavelengths by band index so they line up with pixel_spectra for the line plot
idx_to_wavelength = {i: wavelength for wavelength, i in wavelength_map.items()}
band_wavelengths = [idx_to_wavelength[i] for i in range(image_arr.shape[0])]

# Subplots allows us to plot two things at once! Neat!
fig, axs = plt.subplots(1, 2, figsize=(13, 5), width_ratios=[1.5, 1])

# First, plotting the spectra as a line chart
axs[0].plot(band_wavelengths, pixel_spectra, marker="o", color="#ff4848")
axs[0].set_title("Pixel Spectral Curve")
axs[0].set_xlabel("Band Wavelength (nm)")
axs[0].set_ylabel("Surface Reflectance")

# Secondly, plotting a buffered area around our pixel to visualize what we're looking at!
rgb_crop = np.dstack([
    get_band(650)[PIXEL_ROW - 200 : PIXEL_ROW + 200, PIXEL_COL - 200 : PIXEL_COL + 200],
    get_band(550)[PIXEL_ROW - 200 : PIXEL_ROW + 200, PIXEL_COL - 200 : PIXEL_COL + 200],
    get_band(503)[PIXEL_ROW - 200 : PIXEL_ROW + 200, PIXEL_COL - 200 : PIXEL_COL + 200],
]).astype(np.float32)

brightness_factor = 3  # Bring up the brightness a touch so we have a better visual
rgb_crop_bright = rgb_crop * brightness_factor
rgb_crop_bright = np.clip(rgb_crop_bright, 0, 1)  # keep values within 0-1

axs[1].imshow(rgb_crop_bright)

# We can also add a red box around the pixel we are looking at in the line chart!
pixel_rect = Rectangle(
    (195, 195), 10, 10, linewidth=2, edgecolor="#ff4848", facecolor="none"
)
axs[1].add_patch(pixel_rect)
axs[1].set_title("Pixel Location in Image (RGB)")

plt.tight_layout()
plt.show()
Spectral curve of selected pixel. Red box denotes where the pixel resides.

Spectral curve of selected pixel. Red box denotes where the pixel resides.

Spectral Curve Discussion​

This pixel shows there is vegetation present. The reflectance stays low and fairly flat through the visible range, since chlorophyll is absorbing most of that light, then jumps sharply around 680–730nm into the near-infrared. That jump is called the "red edge," and this makes healthy vegetation easy to pick out from soil or other surfaces!

Because this is L2A surface reflectance data, atmospheric effects have already been largely corrected for, so the curve above is an estimate of surface reflectance after atmospheric correction, much closer to the material's actual reflectance at this pixel than raw radiance would be. Keep in mind that comparing it against reference spectral libraries still requires accounting for differences in available bands and bandwidths between Wyvern's sensor and the library measurements, which we cover in more detail in the Detecting Rare Earth Elements tutorial.

Scatter Plot/Heatmap

Heatmap/scatter plots of pixel distributions are useful for discerning different groups of pixels within an image. For example, one group of pixels may be vegetation, while another group of pixels might be dirt or asphalt.

# Now, we're ready to plot some band combinations!
# We'll start out w/ Bands 649nm and 764nm
band_x, X_IDX = get_band(649, return_index=True)
band_y, Y_IDX = get_band(764, return_index=True)

# Ravel turns our 2D X-Y image array into a single dimensional array
band_x = band_x.ravel()
band_y = band_y.ravel()

# Here we're filtering out our NaNs from the array. We use a single joint mask so that if a
# pixel is NaN in either band, it's dropped from both, keeping band_x and band_y aligned pixel-for-pixel.
valid = ~np.isnan(band_x) & ~np.isnan(band_y)
band_x = band_x[valid]
band_y = band_y[valid]

# This generates a two-dimensional histogram, or summary of our data statistics.
# It allows us to plot a "heatmap" of where points lie between two bands.
# We use percentiles to set the plot visualization range. This just
# keeps a few very bright or very dark pixels from stretching the axes
x_min, x_max = np.percentile(band_x, [0.01, 99.9])
y_min, y_max = np.percentile(band_y, [0.01, 99.9])

heatmap, xedges, yedges = np.histogram2d(
    band_x, band_y, bins=200, range=[[x_min, x_max], [y_min, y_max]]
)
extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]

# We then can plot our heatmap!
fig, ax = plt.subplots(1, 1)
ax.imshow(heatmap.T, extent=extent, origin="lower", cmap="viridis")
plt.xlabel(band_descriptions[X_IDX])
plt.ylabel(band_descriptions[Y_IDX])

# Since we highlighted where our selected pixel was on the image in the above spectral plot, let's
# also highlight where the pixel lies in our scatter plot!
box_size = 0.015
pixel_rect = Rectangle(
    (pixel_spectra[X_IDX] - box_size / 2, pixel_spectra[Y_IDX] - box_size / 2),
    box_size, box_size, linewidth=2, edgecolor="#ff4848", facecolor="none"
)
ax.add_patch(pixel_rect)

plt.title("2 Dimensional Band Scatter Plot")
plt.tight_layout()
plt.show()
649nm & 764nm Scatter Plot. Red box denotes pixel we looked at in the spectra plots.

649nm & 764nm Scatter Plot. Red box denotes pixel we looked at in the spectra plots.

# Let's now do our first and last bands!
band_x, X_IDX = get_band(503, return_index=True)
band_y, Y_IDX = get_band(799, return_index=True)

# Ravel turns our 2D X-Y image array into a single dimensional array
band_x = band_x.ravel()
band_y = band_y.ravel()

# Here we're filtering out our NaNs from the array. We use a single joint mask so that if a
# pixel is NaN in either band, it's dropped from both, keeping band_x and band_y aligned pixel-for-pixel.
valid = ~np.isnan(band_x) & ~np.isnan(band_y)
band_x = band_x[valid]
band_y = band_y[valid]

# This generates a two-dimensional histogram, or summary of our data statistics.
# It allows us to plot a "heatmap" of where points lie between two bands.
# We use percentiles to set the plot visualization range. This just
# keeps a few very bright or very dark pixels from stretching the axes
x_min, x_max = np.percentile(band_x, [0.01, 99.9])
y_min, y_max = np.percentile(band_y, [0.01, 99.9])

heatmap, xedges, yedges = np.histogram2d(
    band_x, band_y, bins=200, range=[[x_min, x_max], [y_min, y_max]]
)
extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]

# We then can plot our heatmap!
fig, ax = plt.subplots(1, 1)
ax.imshow(heatmap.T, extent=extent, origin="lower", cmap="viridis")
plt.xlabel(band_descriptions[X_IDX])
plt.ylabel(band_descriptions[Y_IDX])

# Since we highlighted where our selected pixel was on the image in the above spectral plot, let's
# also highlight where the pixel lies in our scatter plot!
box_size = 0.015
pixel_rect = Rectangle(
    (pixel_spectra[X_IDX] - box_size / 2, pixel_spectra[Y_IDX] - box_size / 2),
    box_size, box_size, linewidth=2, edgecolor="#ff4848", facecolor="none"
)
ax.add_patch(pixel_rect)

plt.title("2 Dimensional Band Scatter Plot")
plt.tight_layout()
plt.show()
503nm & 799nm Scatter Plot. Red box denotes pixel we looked at in the spectra plots.

503nm & 799nm Scatter Plot. Red box denotes pixel we looked at in the spectra plots.

As-generated these plots do not tell us much on their own. In both plots we can see a tight, bright cluster of pixels, this is the dirt/sand around the irrigated fields. The Y-axis in both plots represents a near-infrared band, while the X-axis represents a visible band. Vegetation generally has high near-infrared reflectance and lower visible reflectance, which produces the fainter diagonal spread of pixels extending away from that cluster toward higher Y (and lower X) values. We can see some distinct groups within that spread (feel free to zoom the plot in to this area!), these could correspond to different crop types, growth stages, or vegetation conditions. How cool!

Our vegetation pixel lands in the top left corner of both plots. Plants absorb most visible light for photosynthesis, keeping their visible band reflectance low, while reflecting near-infrared light, pushing their NIR band reflectance high. That combination of low visible and high NIR reflectance is what puts vegetation in that top left region.

We could extend these plots by highlighting areas of high-NDVI, or labeling areas with a shapefile and highlighting these areas to identify pixels that have a matching/similar spectra. We'll leave extending this notebook up to you, dear reader :)