Showcase — the 2023 Canadian wildfires (FIRMS + Sentinel-2 dNBR)¶
2023 was Canada's worst fire season on record — about 10–14 million hectares
burned (~4.7× the 1986–2022 average). This notebook combines NASA FIRMS
active-fire detections with a Sentinel-2 dNBR burn-severity map rendered via
pyramids' Dataset.plot().
What this notebook does¶
- FIRMS — VIIRS active-fire detections over Quebec at the June 2023 peak (point map + per-day counts + FRP distribution; FIRMS is a vector backend → matplotlib).
- dNBR — Sentinel-2 before/after burn severity, rendered with
Dataset.plot()and classified into USGS severity classes.
Needs
FIRMS_MAP_KEYfor FIRMS andGEE_SERVICE_ACCOUNT/GEE_SERVICE_KEY
pyramids-gis[viz]for the Sentinel-2 cells.
Setup¶
First the imports and a bit of noise-suppression. pyramids provides Dataset
(reading + plotting); earthlens provides the unified EarthLens entry point.
import os
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from dotenv import find_dotenv, load_dotenv
from loguru import logger
from pyramids.dataset import Dataset # pyramids' own reading + plotting
from earthlens import EarthLens
warnings.filterwarnings("ignore")
logger.remove()
Credentials¶
The notebook reads its credentials from a .env file at the repository root
(loaded with python-dotenv; your PATH is left untouched). Create that .env
with:
GEE_SERVICE_ACCOUNT="your-sa@your-project.iam.gserviceaccount.com"
GEE_SERVICE_KEY="/path/to/your/service-account-key.json"
FIRMS_MAP_KEY="your-free-NASA-FIRMS-map-key"
load_dotenv(find_dotenv(usecwd=True))
SERVICE_ACCOUNT = os.environ["GEE_SERVICE_ACCOUNT"]
SERVICE_KEY = os.environ["GEE_SERVICE_KEY"]
FIRMS_MAP_KEY = os.environ["FIRMS_MAP_KEY"]
if any(x is None for x in (SERVICE_ACCOUNT, SERVICE_KEY, FIRMS_MAP_KEY)):
raise ValueError("Missing credentials in .env")
1 · FIRMS active-fire detections¶
FIRMS serves active-fire pixels as a table, so earthlens treats it as a
vector backend: load() returns a GeoDataFrame, one row per detection. We
build the request first — source, sensor (variables), date window, and bounding
box over Quebec.
fires = EarthLens(
data_source="firms",
variables=["VIIRS_SNPP_SP"],
start="2023-06-02",
end="2023-06-08",
aoi=[-80.0, 48.0, -72.0, 53.0],
)
Authentication and loading are kept as separate steps so each is easy to read
and re-run: authenticate() resolves the FIRMS MAP_KEY, then load() fetches
the detections into memory.
fires.authenticate(api_key=FIRMS_MAP_KEY)
data = fires.load(progress_bar=False)
print(f"{len(data):,} fire detections")
Where the fires were¶
Each detection carries a fire-radiative-power (FRP) value; we colour the scatter
by log10(FRP) so both smouldering and intense pixels stay visible.
fig, ax = plt.subplots(figsize=(8, 6))
sc = ax.scatter(
data.longitude,
data.latitude,
c=np.log10(data.frp.clip(lower=1)),
cmap="inferno",
s=6,
alpha=0.6,
)
ax.set(
xlabel="longitude",
ylabel="latitude",
title=f"VIIRS active fire, 2–8 Jun 2023 ({len(data):,})",
)
fig.colorbar(sc, ax=ax, label="log10 FRP (MW)")
plt.tight_layout()
plt.show()
How many per day¶
Grouping the detections by acquisition date shows how the outbreak ramped over the week; the FRP summary quantifies the total radiative power.
per_day = (
pd.to_datetime(data.acq_datetime, utc=True).dt.date.value_counts().sort_index()
)
fig, ax = plt.subplots(figsize=(8, 5))
ax.bar([str(x) for x in per_day.index], per_day.values, color="firebrick")
ax.set(ylabel="detections", title="Detections per day")
plt.setp(ax.get_xticklabels(), rotation=30, ha="right")
plt.tight_layout()
plt.show()
print(f"median FRP {data.frp.median():.1f} MW | total {data.frp.sum() / 1e3:.1f} GW")
2 · Burn severity (dNBR)¶
NBR = (NIR − SWIR) / (NIR + SWIR); dNBR = NBR_before − NBR_after is the
standard burn-severity index. We pull a median Sentinel-2 composite for a
pre-fire and a post-fire window, difference them, and render the result with
Dataset.plot().
A helper for one NBR window¶
nbr() builds a GEE request for BURN_AOI, authenticates, loads the single
composite raster, and returns both the NBR array and the underlying Dataset (we
need its geotransform/EPSG later to wrap the dNBR back into a plottable raster).
Note the construct → authenticate → load steps are kept on separate lines.
BURN_AOI = [-77.0, 51.0, -76.0, 52.0]
def nbr(start, end):
"""Median Sentinel-2 NBR over BURN_AOI for [start, end]; returns (array, dataset)."""
scene = EarthLens(
data_source="gee",
cadence="raw",
dataset="COPERNICUS/S2_SR_HARMONIZED",
variables=["B8", "B12"],
aoi=BURN_AOI,
start=start,
end=end,
scale=80,
reducer="median",
)
scene.authenticate(service_account=SERVICE_ACCOUNT, service_key=SERVICE_KEY)
rasters = scene.load()
raster = rasters[0]
bands = np.asarray(raster.read_array(), dtype="float32")
nbr_arr = (bands[0] - bands[1]) / (bands[0] + bands[1] + 1e-6)
return nbr_arr, raster
Before vs after, and the difference¶
Run the helper for a May (pre-fire) and a September (post-fire) window, then difference them on their overlapping extent.
nbr_pre, raster_pre = nbr("2023-05-01", "2023-05-31")
nbr_post, _ = nbr("2023-09-01", "2023-09-30")
h = min(nbr_pre.shape[0], nbr_post.shape[0])
w = min(nbr_pre.shape[1], nbr_post.shape[1])
dnbr = nbr_pre[:h, :w] - nbr_post[:h, :w]
Render the dNBR raster¶
The dNBR is a single-band index array. Wrapping it back into a Dataset (reusing
the source geotransform/EPSG) lets pyramids colour-map and georeference it via
Dataset.plot().
dnbr_ds = Dataset.create_from_array(
arr=dnbr[None, :, :], geo=raster_pre.geotransform, epsg=raster_pre.epsg
)
fig, ax = plt.subplots(figsize=(7, 6))
dnbr_ds.plot(
band=0, ax=ax, cmap="RdYlGn_r", title="dNBR burn severity", vmin=-0.1, vmax=1.0
)
plt.tight_layout()
plt.show()
Area by burn-severity class¶
Binning dNBR into the USGS severity classes and multiplying pixel counts by the pixel area (80 m × 80 m) gives the burned area per class.
bounds = [-0.1, 0.1, 0.27, 0.44, 0.66, 1.3]
names = ["unburned", "low", "mod-low", "mod-high", "high"]
px_km2 = (80 * 80) / 1e6
areas = [
np.sum((dnbr >= bounds[i]) & (dnbr < bounds[i + 1])) * px_km2 for i in range(5)
]
fig, ax = plt.subplots(figsize=(7, 5))
ax.bar(names, areas, color=["#2c7bb6", "#abd9e9", "#ffffbf", "#fdae61", "#d7191c"])
ax.set_ylabel("area (km^2)")
ax.set_title("Area by burn-severity class")
plt.setp(ax.get_xticklabels(), rotation=30, ha="right")
plt.tight_layout()
plt.show()
print(
f"burned area in this tile: {sum(areas[1:]):,.0f} km^2 "
f"({areas[3] + areas[4]:,.0f} km^2 mod-high to high)"
)