Lab 21 - Change Detection: Mapping Burn Severity

Before You Start

• Prerequisites: Spectral Indices module, Cloud Masking module, Change Detection module.
• Estimated time: 45 minutes
• You will need: GEE Code Editor access

Lab Instructions

The 2021 Caldor Fire burned over 221,000 acres in El Dorado County, California, threatening the Lake Tahoe basin. Let's map its burn severity using Landsat 8 imagery.

Part 1 - Load the pre-fire imagery

We'll start by creating a cloud-masked median composite of Landsat 8 imagery from before the fire. The Caldor Fire started in mid-August 2021, so we'll use May through July for the pre-fire period.

// Cloud masking function for Landsat 8 SR
function maskL8sr(image) {
  var qaMask = image.select('QA_PIXEL').bitwiseAnd(1 << 3).eq(0)
      .and(image.select('QA_PIXEL').bitwiseAnd(1 << 4).eq(0));
  return image.updateMask(qaMask)
      .multiply(0.0000275).add(-0.2);
}

// Define the study area (Caldor Fire region)
var fireArea = ee.Geometry.Point([-120.5, 38.7]);
var studyRegion = fireArea.buffer(25000); // 25 km buffer

// Load pre-fire imagery (May - July 2021)
var preFire = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
  .filterDate('2021-05-01', '2021-07-31')
  .filterBounds(studyRegion)
  .filter(ee.Filter.lt('CLOUD_COVER', 20))
  .map(maskL8sr)
  .median();

print('Pre-fire composite created');
Map.centerObject(studyRegion, 11);
Map.addLayer(preFire, {
  bands: ['SR_B4', 'SR_B3', 'SR_B2'],
  min: 0, max: 0.3
}, 'Pre-fire True Color');

Understanding the code:

• maskL8sr() removes cloud and shadow pixels using the QA_PIXEL band, then applies the scale/offset correction.
• .buffer(25000) creates a 25 km radius around the center point to capture the full fire perimeter.
• .median() combines multiple images into a single cloud-free composite. Think of it like stacking photos and picking the most typical value for each pixel.

Part 2 - Load the post-fire imagery

Now let's load imagery from after the fire was contained. The Caldor Fire was fully contained in October 2021, so we'll use October through December for the post-fire period.

// Load post-fire imagery (October - December 2021)
var postFire = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
  .filterDate('2021-10-01', '2021-12-31')
  .filterBounds(studyRegion)
  .filter(ee.Filter.lt('CLOUD_COVER', 20))
  .map(maskL8sr)
  .median();

Map.addLayer(postFire, {
  bands: ['SR_B4', 'SR_B3', 'SR_B2'],
  min: 0, max: 0.3
}, 'Post-fire True Color');

Understanding the code:

• The post-fire composite uses the same cloud masking and compositing approach as the pre-fire composite.
• Toggle between the pre-fire and post-fire layers in the Layers panel to visually compare the two scenes. The difference is dramatic.

Part 3 - Calculate NBR for both dates

The Normalized Burn Ratio uses NIR (Band 5) and SWIR (Band 7) to detect burn conditions. Here's the key insight: healthy vegetation has high NIR reflectance and low SWIR reflectance. Burned areas show the opposite pattern. NBR captures that contrast beautifully.

// Calculate NBR for pre-fire composite
// NBR = (NIR - SWIR2) / (NIR + SWIR2)
var preNBR = preFire.normalizedDifference(['SR_B5', 'SR_B7']).rename('preNBR');

// Calculate NBR for post-fire composite
var postNBR = postFire.normalizedDifference(['SR_B5', 'SR_B7']).rename('postNBR');

// Visualize NBR images
Map.addLayer(preNBR, {min: -1, max: 1,
  palette: ['red', 'yellow', 'green']}, 'Pre-fire NBR');
Map.addLayer(postNBR, {min: -1, max: 1,
  palette: ['red', 'yellow', 'green']}, 'Post-fire NBR');

Understanding the code:

• .normalizedDifference(['SR_B5', 'SR_B7']) computes (B5 - B7) / (B5 + B7), which is the standard NBR formula.
• Pre-fire NBR should be high (green) in forested areas. Post-fire NBR should drop sharply (red) where the fire burned. That drop is our signal.

Part 4 - Compute dNBR

Now let's subtract the post-fire NBR from the pre-fire NBR. This produces the differenced NBR. High positive values mean the fire hit hard. Values near zero mean nothing changed. Think of it like a "damage score" for every pixel.

// Compute dNBR (pre minus post)
var dNBR = preNBR.subtract(postNBR).rename('dNBR');

// Clip to study region
dNBR = dNBR.clip(studyRegion);

// Display the continuous dNBR
Map.addLayer(dNBR, {
  min: -0.5, max: 1.3,
  palette: ['green', 'lightyellow', 'yellow', 'orange', 'red', 'darkred']
}, 'dNBR');

// Print dNBR statistics
print('dNBR statistics:', dNBR.reduceRegion({
  reducer: ee.Reducer.minMax(),
  geometry: studyRegion,
  scale: 30,
  maxPixels: 1e9
}));

Understanding the code:

• preNBR.subtract(postNBR) computes the difference. Where fire burned, post-fire NBR is much lower than pre-fire NBR, producing a large positive dNBR value.
• .clip(studyRegion) restricts the output to our area of interest.

Part 5 - Classify burn severity

Here's the cool part. We'll apply the standard USGS burn severity thresholds to turn our continuous dNBR into a classified map. Each pixel gets assigned to one of seven severity classes, just like the professionals do it.

// Classify dNBR into USGS burn severity categories
var burnSeverity = ee.Image(0)
  .where(dNBR.lt(-0.25), 1)                           // Enhanced regrowth (high)
  .where(dNBR.gte(-0.25).and(dNBR.lt(-0.1)), 2)       // Enhanced regrowth (low)
  .where(dNBR.gte(-0.1).and(dNBR.lt(0.1)), 3)         // Unburned
  .where(dNBR.gte(0.1).and(dNBR.lt(0.27)), 4)         // Low severity
  .where(dNBR.gte(0.27).and(dNBR.lt(0.44)), 5)        // Moderate-low severity
  .where(dNBR.gte(0.44).and(dNBR.lt(0.66)), 6)        // Moderate-high severity
  .where(dNBR.gte(0.66), 7)                            // High severity
  .clip(studyRegion)
  .selfMask();

// Display classified map
Map.addLayer(burnSeverity, {
  min: 1, max: 7,
  palette: [
    '#38a800',   // 1 - Enhanced regrowth (high)
    '#6fc400',   // 2 - Enhanced regrowth (low)
    '#ffffbe',   // 3 - Unburned
    '#ffff00',   // 4 - Low severity
    '#ffaa00',   // 5 - Moderate-low severity
    '#ff5500',   // 6 - Moderate-high severity
    '#e60000'    // 7 - High severity
  ]
}, 'Burn Severity Classes');

Understanding the code:

• ee.Image(0) creates a blank image that serves as the starting canvas for classification.
• .where(condition, value) sets pixels to a class value wherever the condition is true. It's like painting by numbers.
• .selfMask() removes pixels with a value of 0 (any that were not classified).
• The palette assigns a color gradient from green (regrowth) through yellow (unburned) to red (high severity).

Part 6 - Calculate burned area

Finally, let's calculate the total area burned in each severity class. This is the type of summary statistic that appears in official fire damage reports.

// Calculate area for each burn class (in hectares)
// Focus on classes 4-7 (actual burn, not regrowth or unburned)
var burnedMask = burnSeverity.gte(4);
var burnedArea = burnedMask.multiply(ee.Image.pixelArea());

var totalBurned = burnedArea.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: studyRegion,
  scale: 30,
  maxPixels: 1e9
});

print('Total burned area (hectares):',
  ee.Number(totalBurned.get('constant')).divide(10000));

// Calculate high-severity area specifically (class 7)
var highSeverityMask = burnSeverity.eq(7);
var highArea = highSeverityMask.multiply(ee.Image.pixelArea());

var highTotal = highArea.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: studyRegion,
  scale: 30,
  maxPixels: 1e9
});

print('High severity area (hectares):',
  ee.Number(highTotal.get('constant')).divide(10000));

Understanding the code:

• ee.Image.pixelArea() creates an image where each pixel's value is its area in square meters.
• Multiplying by the binary mask keeps only the pixels we care about.
• ee.Reducer.sum() adds up all the pixel areas.
• Dividing by 10,000 converts from square meters to hectares.