Quantifying the severity of forest fires using LandTrendr

---

The LandTrendr algorithm (Kennedy et al., 2010) is used to detect patterns of disturbance and recovery in satellite image time series. It works by breaking the history of each pixel into straight-line segments.

Some segments are flat, showing little change over several years. These represent stable conditions. A disturbance, such as logging or fire, appears as a short segment with a steep slope, showing a sudden drop in the pixel value. Recovery, if it happens, is seen as a longer segment with a gentle slope, moving back toward the original value.

The algorithm can measure how large a disturbance was (magnitude), how long it lasted (duration), and how fast recovery occurred (slope). This makes it useful for tracking the timing, severity, and outcomes of changes in forests and other landscapes.

LandTrendr needs at least one image per year and works best with six or more years of data. By summarising changes over time, it helps scientists and managers understand both sudden events and slow trends in the environment.

In this example, I use the Normalized Burn Ratio (NBR) index, which is more effective than NDVI for mapping forest fire impacts. Here, I show disturbance in Kangaroo Island, Australia, caused by the 2019 forest fires. The algorithm detected a sharp drop in NBR (disturbance) followed by a slow upward trend (recovery). This helps measure the magnitude, duration, and speed of recovery for each pixel.

The LandTrendr algorithm detects trends in disturbance (a short period of change) and recovery (the longer process of returning to the original state) of pixels (Kennedy et al., 2010). It is based on the idea that a pixel’s history can be divided into several straight-line segments over time. Long periods with little change appear as flat lines. A disturbance shows as a short, often steeply sloped line. Recovery, such as regrowth after a forest fire, appears as a longer, gently sloped line moving back toward the original value.

For mapping forest fires, I use the Normalized Burn Ratio (NBR), which is more effective than NDVI. NBR is calculated as:

$$ \text{NBR} = \frac{\text{NIR} - \text{SWIR2}}{\text{NIR} + \text{SWIR2}} $$

where NIR is the near-infrared band (Sentinel-2 B8, Landsat B5) and SWIR2 is the shortwave infrared band 2 (Sentinel-2 B12, Landsat B7). Fires reduce NIR reflectance and increase SWIR2 reflectance, making changes easy to detect.

The algorithm needs one image for each year in the analysis and works best with at least six years of data. It is widely used to detect the magnitude, duration, and year of disturbances such as fires, logging, or storms.

Figure: Conceptual illustration of the LandTrendr algorithm. The black dotted line shows the original annual values of a spectral index (e.g., NBR, NDVI) for a single pixel. The red line is the fitted LandTrendr model, which breaks the pixel’s history into straight-line segments. Periods of stability have little slope, disturbances appear as short steep declines, and recovery is shown as a longer gentle increase. The algorithm can measure disturbance magnitude, duration, and rate of recovery. Image source: Mugiraneza T et al., (2020)

---

Import libraries and set working directories

# set working directory
import os
os.chdir("/media/roshan/data/03_projects/18_github_projects")
print("Current working directory:", os.getcwd())

Current working directory: /media/roshan/data/03_projects/18_github_projects

Connect to Earth Engine API

# import required libraries
import ee # Earth Engine Python API
import geemap # geemap for Earth Engine visualization and export

# --- Auth/Init ---
ee.Authenticate()
ee.Initialize()

Process Sentinel data from Earth Engine

# --- AOI & dates (Kangaroo Island example) ---
aoi = ee.Geometry.Rectangle([136.6, -36.2, 137.1, -35.7])
START = '2018-01-01'
END   = '2023-12-31'
MAX_CLOUD_PROB = 20  # cloud cover percentage

# --- Collections ---
s2_sr = (ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
         .filterBounds(aoi)
         .filterDate(START, END))

s2_cp = (ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY')
         .filterBounds(aoi)
         .filterDate(START, END))

# --- maskEdges: use 20m & 60m masks to clean 10/20m edges (doc-style) ---
def mask_edges(img):
    return (img
            .updateMask(img.select('B8A').mask())  # 20 m
            .updateMask(img.select('B9').mask()))  # 60 m

s2_sr = s2_sr.map(mask_edges)

# --- Join SR with cloud prob using saveFirst('cloud_mask') like the docs ---
join = ee.Join.saveFirst('cloud_mask')
cond = ee.Filter.equals(leftField='system:index', rightField='system:index')
s2_joined = ee.ImageCollection(join.apply(primary=s2_sr, secondary=s2_cp, condition=cond))

# --- Cloud mask using attached probability image (doc-style) ---
def mask_clouds(img):
    clouds = ee.Image(img.get('cloud_mask')).select('probability')
    is_clear = clouds.lt(MAX_CLOUD_PROB)
    return img.updateMask(is_clear)

# --- Build NBR-only collection (single band) ---
# Use 20 m bands: NBR = (B8A - B12) / (B8A + B12)
def add_nbr(img):
    nbr = img.normalizedDifference(['B8A', 'B12']).rename('NBR')
    return img.addBands(nbr)

s2_nbr = (s2_joined
          .map(mask_clouds)
          .map(add_nbr)
          .select(['NBR']))

# --- Annual NBR composites (median) with time_start set ---
def yearly_nbr(year):
    start = ee.Date.fromYMD(year, 1, 1)
    end   = ee.Date.fromYMD(year, 12, 31)
    nbr   = s2_nbr.filterDate(start, end).median().set('system:time_start', start.millis())
    return nbr.clip(aoi)

years = list(range(2018, 2024))
nbr_series = ee.ImageCollection([yearly_nbr(y) for y in years])

print('Annual images in NBR series:', nbr_series.size().getInfo())

Annual images in NBR series: 6

Running the LandTrendr Algorithm on NBR Time Series (2018–2024)

# Run LandTrendr temporal segmentation algorithm
lt = ee.Algorithms.TemporalSegmentation.LandTrendr(**{
    'timeSeries': nbr_series,           # Input time series (NBR images by year)
    'maxSegments': 6,                   # Max number of fitted segments
    'spikeThreshold': 0.9,              # Remove short-term noise spikes
    'vertexCountOvershoot': 3,          # Allow extra vertices before pruning
    'preventOneYearRecovery': True,     # Avoid classifying 1-year change as recovery
    'recoveryThreshold': 0.25,          # Minimum slope (positive) to define recovery
    'pvalThreshold': 0.05,              # Significance level for fitting segments
    'bestModelProportion': 0.75,        # Retain models within 75% of the best fit
    'minObservationsNeeded': 6          # Minimum years required in the time series
})

# Print the band names of the output LandTrendr image
print('LT bands:', lt.bandNames().getInfo())  # Expect: ['LandTrendr']

LT bands: ['LandTrendr', 'rmse']

Extracting LandTrendr Vertex Years and Fitted Values

def lt_vertices(lt_image):
    arr    = lt_image.select('LandTrendr')     # Select LandTrendr array output
    years  = arr.arraySlice(0, 0, 1)           # Slice band axis to get year values
    fitted = arr.arraySlice(0, 1, 2)           # Slice to get fitted index values
    isvert = arr.arraySlice(0, 2, 3)           # Slice to get vertex flags (0 or 1)
    
    # Mask to keep only vertices (flag=1) for both years and fitted values
    return (years.arrayMask(isvert), fitted.arrayMask(isvert))  

# Apply function to extract vertex years and fitted values
years_v, fitted_v = lt_vertices(lt)

Calculating LandTrendr Segment Metrics

def lt_segments(years_v, fitted_v):
    # Start and end years for each segment
    y0 = years_v.arraySlice(1, 0, -1)          # segment start years
    y1 = years_v.arraySlice(1, 1, None)        # segment end years
    
    # Start and end fitted values (index values) for each segment
    v0 = fitted_v.arraySlice(1, 0, -1)         # fitted value at start
    v1 = fitted_v.arraySlice(1, 1, None)       # fitted value at end

    # Segment metrics
    dur   = y1.subtract(y0)                    # duration (years)
    mag   = v1.subtract(v0)                    # magnitude of change (negative = loss)
    slope = mag.divide(dur)                    # rate of change per year
    rel   = mag.divide(v0.abs()).multiply(100) # relative change (%)

    # Stack attributes along axis 0; segments remain along axis 1
    seg = (y0.arrayCat(y1,0)
             .arrayCat(v0,0).arrayCat(v1,0)
             .arrayCat(dur,0).arrayCat(mag,0)
             .arrayCat(slope,0).arrayCat(rel,0))
    
    return seg  # band order: [y0, y1, v0, v1, dur, mag, slope, rel]

# Generate segment metrics from vertex years and fitted values
seg = lt_segments(years_v, fitted_v)

Function to calculate loss and recovery

def lt_greatest_loss_recovery(seg, metric='mag'):
    """
    Finds the greatest loss (negative change) and greatest recovery (positive change)
    from LandTrendr segments in a single pass.

    Parameters:
        seg    : ee.Image (stacked segment metrics from lt_segments)
        metric : str, one of {'mag','slope','rel'}

    Returns:
        ee.Image with 10 bands:
        loss_year, loss_mag, loss_dur, loss_slope, loss_rel,
        rec_year,  rec_mag,  rec_dur,  rec_slope,  rec_rel
    """
    IDX = {'y0':0,'y1':1,'v0':2,'v1':3,'dur':4,'mag':5,'slope':6,'rel':7}

    def pick_extreme(is_loss=True):
        target = seg.arraySlice(0, IDX[metric], IDX[metric]+1)  # [1, nSeg]

        if is_loss:
            extreme_val = target.arrayReduce(ee.Reducer.min(), [1])
        else:
            extreme_val = target.arrayReduce(ee.Reducer.max(), [1])

        extreme_1d = extreme_val.arrayProject([0])
        seglen     = target.arrayLength(1)
        extreme_rep = extreme_1d.arrayRepeat(1, seglen)
        is_extreme  = target.eq(extreme_rep)  # boolean mask

        def pick(i):
            attr   = seg.arraySlice(0, i, i+1)
            masked = attr.arrayMask(is_extreme)
            return masked.arrayReduce(ee.Reducer.firstNonNull(), [1])

        def squeeze(arr_img, name):
            return arr_img.arrayFlatten([['row0'], ['col0']]).rename(name)

        year  = squeeze(pick(IDX['y1']),  'loss_year' if is_loss else 'rec_year').toInt16()
        mag   = squeeze(pick(IDX['mag']), 'loss_mag'  if is_loss else 'rec_mag')
        dur   = squeeze(pick(IDX['dur']), 'loss_dur'  if is_loss else 'rec_dur')
        slope = squeeze(pick(IDX['slope']),'loss_slope' if is_loss else 'rec_slope')
        rel   = squeeze(pick(IDX['rel']), 'loss_rel'  if is_loss else 'rec_rel')

        return ee.Image.cat([year, mag, dur, slope, rel])

    # Get loss (most negative) and recovery (most positive)
    loss_img = pick_extreme(is_loss=True)
    rec_img  = pick_extreme(is_loss=False)

    return loss_img.addBands(rec_img)

Apply the function to the segment

loss_recovery = lt_greatest_loss_recovery(seg, metric='mag')
print('Loss+Recovery bands:', loss_recovery.bandNames().getInfo()) 

Loss+Recovery bands: ['loss_year', 'loss_mag', 'loss_dur', 'loss_slope', 'loss_rel', 'rec_year', 'rec_mag', 'rec_dur', 'rec_slope', 'rec_rel']

Export the image to google drive

# Export all bands from 'loss' image as one multi-band GeoTIFF
loss_recovery = loss_recovery.toFloat().clip(aoi).unmask()  # ensure no NaNs, clip to AOI
task = ee.batch.Export.image.toDrive(
    image=loss_recovery,
    description='LT_greatest_loss_recovery_export',
    folder='GEE_exports',
    fileNamePrefix='lt_loss_recovery_metrics',
    region=aoi,
    scale=100,             # downscale for faster export
    maxPixels=1e13
)
task.start()
print("Export started. Check your Earth Engine Tasks tab.")


# Need to manually download the exported file from Google Drive
# or use geemap to download it directly if you have geemap installed
# (but geempap was not working in this case)

Export started. Check your Earth Engine Tasks tab.

Plot the results

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

tif_path = "lt_loss_recovery_metrics.tif"

# bands to plot (in this order)
band_names = ["loss_mag", "loss_slope", "rec_mag", "rec_slope"]

# Map names → band index in the TIFF (1-based)
# Adjust if your export order differs.
band_idx = {
    "loss_year": 1, "loss_mag": 2, "loss_dur": 3, "loss_slope": 4, "loss_rel": 5,
    "rec_year":  6, "rec_mag":  7, "rec_dur":  8, "rec_slope":  9, "rec_rel":  10
}

with rasterio.open(tif_path) as src:
    nd = src.nodata
    # read only the requested bands
    bands = [src.read(band_idx[name]) for name in band_names]

# mask helper
def valid_mask(arr, nd):
    return (~np.isclose(arr, nd)) if nd is not None else np.isfinite(arr)

# robust limits for continuous bands
def robust_limits(b, m):
    vals = b[m]
    if vals.size == 0:
        return None, None
    lo = np.nanpercentile(vals, 2)
    hi = np.nanpercentile(vals, 98)
    if np.isfinite(lo) and np.isfinite(hi) and lo != hi:
        return float(lo), float(hi)
    return None, None

# choose vmin/vmax per band
vmins, vmaxes = [], []
for i, name in enumerate(band_names):
    b = bands[i]
    m = valid_mask(b, nd)
    # for slopes, use robust symmetric limits around 0
    if "slope" in name:
        lo, hi = robust_limits(b, m)
        if lo is None or hi is None:
            vmins.append(None); vmaxes.append(None)
        else:
            mx = max(abs(lo), abs(hi))
            vmins.append(-mx); vmaxes.append(mx)
    else:
        lo, hi = robust_limits(b, m)
        vmins.append(lo); vmaxes.append(hi)

# colormaps
cmaps = {
    "loss_mag":  "RdBu", 
    "loss_slope":"PuOr",
    "rec_mag":   "YlGn",
    "rec_slope": "PuOr",
}

# plot 2x2
fig, axes = plt.subplots(2, 2, figsize=(12, 9), constrained_layout=True)
axes = axes.ravel()

for i, name in enumerate(band_names):
    b = bands[i]
    img = np.where(valid_mask(b, nd), b, np.nan)  # set nodata to NaN
    im = axes[i].imshow(img, cmap=cmaps.get(name, "viridis"),
                        vmin=vmins[i], vmax=vmaxes[i])
    axes[i].set_title(name)
    axes[i].axis("off")
    cbar = fig.colorbar(im, ax=axes[i], shrink=0.8)
    cbar.ax.set_ylabel(name, rotation=270, labelpad=12)

plt.show()

---

Does areas with higher loss magnitude also experience higher recovery (in this example only short-term recovery)

import rasterio
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import linregress

tif_path = "lt_loss_recovery_metrics.tif"

# Band positions (adjust if needed)
band_idx = {
    "loss_mag": 2,  # loss magnitude band number
    "rec_mag":  7   # recovery magnitude band number
}

with rasterio.open(tif_path) as src:
    nd = src.nodata
    loss_mag = src.read(band_idx["loss_mag"])
    rec_mag  = src.read(band_idx["rec_mag"])

# Mask nodata
mask = np.isfinite(loss_mag) & np.isfinite(rec_mag)
if nd is not None:
    mask &= (loss_mag != nd) & (rec_mag != nd)

# Flatten & invert loss magnitude (negative → positive)
x = -loss_mag[mask].flatten()  # bigger = more loss
y = rec_mag[mask].flatten()    # bigger = more recovery

# Regression stats
slope, intercept, r_value, p_value, std_err = linregress(x, y)
print(f"Regression: y = {slope:.3f}x + {intercept:.3f}, r = {r_value:.2f}, p = {p_value:.2e}")

# Scatter with regression line
plt.figure(figsize=(8, 6))
sns.regplot(x=x, y=y, scatter_kws={'s': 2, 'alpha': 0.3}, line_kws={'color': 'red'}, ci=None)
plt.xlabel("Loss Magnitude (absolute)")
plt.ylabel("Recovery Magnitude")
plt.title("Relationship between Loss and Recovery Magnitudes")
plt.grid(True, linestyle='--', alpha=0.5)
plt.show()

Regression: y = 0.321x + 0.104, r = 0.67, p = 0.00e+00

Relationship Between Loss Magnitude and Recovery Magnitude

• Regression equation: ( y = 0.321x + 0.104 )
• Correlation coefficient (r): 0.67
• p-value: ( < 0.001 )
• Interpretation:
  • Areas with greater loss magnitude tend to show greater recovery magnitude.
  • Slope < 1 indicates recovery is weaker than loss within the observation period.
  • Positive intercept suggests some baseline recovery even in areas with minimal loss.
• Ecological reasoning:
  • Large disturbances can create conditions favorable for regrowth (e.g., increased light, nutrient pulses).
  • Recovery is constrained by ecosystem resilience, seed availability, and time since loss.
• Caveats:
  • Potential LandTrendr fitting artefacts.
  • Recovery may be underestimated if observation window is short after loss.
  • Relationship may vary by vegetation type and climate zone.

---

References

Kennedy, R. E., Yang, Z., & Cohen, W. B. (2010). Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrendr—Temporal segmentation algorithms. Remote Sensing of Environment, 114(12), 2897-2910.

Mugiraneza, T., Nascetti, A., & Ban, Y. (2020). Continuous monitoring of urban land cover change trajectories with landsat time series and landtrendr-google earth engine cloud computing. Remote Sensing, 12(18), 2883.