Case Study: Marshall Fire

Reconstruction using real landscape, fuel, and weather data

On December 30, 2021, the Marshall Fire ignited in Boulder County, Colorado, driven by sustained winds of 50–70+ mph from the west-southwest. Fueled primarily by dry grasslands, the fire burned roughly 6,000 acres and destroyed over 1,000 structures in the communities of Louisville and Superior, making it the most destructive wildfire in Colorado history.

This tutorial reconstructs the fire using real data from federal sources:

Layer	Source
Fuel models (FBFM13)	LANDFIRE via Landfire.jl
Elevation (30 m DEM)	LANDFIRE via Landfire.jl
Fire perimeter	NIFC WFIGS
Wind observations	KBDU ASOS (Boulder Municipal Airport)

Code

# Load Julia startup.jl if LANDFIRE_EMAIL not set (Quarto may skip startup)
if !haskey(ENV, "LANDFIRE_EMAIL")
    startup = joinpath(first(DEPOT_PATH), "config", "startup.jl")
    isfile(startup) && include(startup)
end

using Elmfire
using Landfire
using ArchGDAL
using GeoFormatTypes
using Downloads
using JSON3
using CairoMakie, GeoMakie
using Random

CairoMakie.activate!(visible = false)  # non-interactive backend — prevents GUI windows during render
Random.seed!(2021)

TaskLocalRNG()

Download Real Data

Landscape: LANDFIRE Fuel Models and Elevation

We use Landfire.jl to download 30 m FBFM13 fuel models and elevation for the fire area.

Download LANDFIRE FBFM13 and elevation

# Bounding box covering the Marshall Fire (with buffer west of ignition)
aoi = "-105.260 39.924 -105.126 39.992"

# Download FBFM13 (Anderson 13 fuel models)
fbfm13_prods = Landfire.products(layer="FBFM13", conus=true)
fbfm13_ds = Landfire.Dataset(fbfm13_prods, aoi)
fbfm13_tif = get(fbfm13_ds)

# Download elevation
elev_prods = Landfire.products(layer="ELEV", conus=true)
elev_ds = Landfire.Dataset(elev_prods, aoi)
elev_tif = get(elev_ds)

println("Fuel models: $fbfm13_tif")
println("Elevation:   $elev_tif")

[ Info: Downloading dataset to /home/runner/.julia/scratchspaces/16f38fd6-0379-4569-89bb-2d3d56e50de8/landfire_data/job_6231172435630567030.zip
[ Info: Job submitted.  View job messages at: https://lfps.usgs.gov/job/66172eb1-fcee-4a17-a6bb-235281342e2d
[ Info: Submitted job with ID: 66172eb1-fcee-4a17-a6bb-235281342e2d.  Checking job every 5 seconds.
.....[ Info: Job Status: Pending
.....[ Info: Job Status: Pending
.....[ Info: Job Status: Executing
.....[ Info: Job Status: Executing
.....[ Info: Job Status: Executing
.....[ Info: Job Status: Succeeded
[ Info: Extracting dataset to /home/runner/.julia/scratchspaces/16f38fd6-0379-4569-89bb-2d3d56e50de8/landfire_data/job_6231172435630567030

7-Zip (a) 25.01 (x64) : Copyright (c) 1999-2025 Igor Pavlov : 2025-08-03
 64-bit locale=C.UTF-8 Threads:4 OPEN_MAX:65536, ASM

Scanning the drive for archives:
1 file, 33307 bytes (33 KiB)

Extracting archive: /home/runner/.julia/scratchspaces/16f38fd6-0379-4569-89bb-2d3d56e50de8/landfire_data/job_6231172435630567030.zip
--
Path = /home/runner/.julia/scratchspaces/16f38fd6-0379-4569-89bb-2d3d56e50de8/landfire_data/job_6231172435630567030.zip
Type = zip
Physical Size = 33307

Everything is Ok

Files: 6
Size:       32329
Compressed: 33307
[ Info: Downloading dataset to /home/runner/.julia/scratchspaces/16f38fd6-0379-4569-89bb-2d3d56e50de8/landfire_data/job_15346698323572098831.zip
[ Info: Job submitted.  View job messages at: https://lfps.usgs.gov/job/014d04e8-9b6d-4662-a5f0-2c91a4a99b9f
[ Info: Submitted job with ID: 014d04e8-9b6d-4662-a5f0-2c91a4a99b9f.  Checking job every 5 seconds.
.....[ Info: Job Status: Pending
.....[ Info: Job Status: Executing
.....[ Info: Job Status: Succeeded
[ Info: Extracting dataset to /home/runner/.julia/scratchspaces/16f38fd6-0379-4569-89bb-2d3d56e50de8/landfire_data/job_15346698323572098831

7-Zip (a) 25.01 (x64) : Copyright (c) 1999-2025 Igor Pavlov : 2025-08-03
 64-bit locale=C.UTF-8 Threads:4 OPEN_MAX:65536, ASM

Scanning the drive for archives:
1 file, 88411 bytes (87 KiB)

Extracting archive: /home/runner/.julia/scratchspaces/16f38fd6-0379-4569-89bb-2d3d56e50de8/landfire_data/job_15346698323572098831.zip
--
Path = /home/runner/.julia/scratchspaces/16f38fd6-0379-4569-89bb-2d3d56e50de8/landfire_data/job_15346698323572098831.zip
Type = zip
Physical Size = 88411

Everything is Ok

Files: 6
Size:       87433
Compressed: 88411
Fuel models: /home/runner/.julia/scratchspaces/16f38fd6-0379-4569-89bb-2d3d56e50de8/landfire_data/job_6231172435630567030/j264635b71cc04f64adc8da3bb4191137.tif
Elevation:   /home/runner/.julia/scratchspaces/16f38fd6-0379-4569-89bb-2d3d56e50de8/landfire_data/job_15346698323572098831/ja67b6f7bbd2543e089373a2bf9035795.tif

Process Landscape into Simulation Inputs

Both rasters are in the same projected CRS (NAD83 CONUS Albers, 30 m cells), so they align perfectly.

Read rasters and compute slope/aspect

# Read fuel raster and map nonburnable codes to 256
# ArchGDAL.read returns (width, height) = (ncols, nrows), matching Elmfire convention
fuel_ids = ArchGDAL.read(fbfm13_tif) do ds
    Int.(ArchGDAL.read(ArchGDAL.getband(ds, 1)))
end

# Map LANDFIRE special codes:
#   91 (urban) → 2 (short grass): the Marshall Fire burned through residential areas
#     with dry grass, so treating urban as short grass is a reasonable approximation
#   93 (agriculture) → 1 (short grass): crop stubble burns like short grass
#   98 (water), 99 (bare ground) → 256 (nonburnable)
for i in eachindex(fuel_ids)
    v = fuel_ids[i]
    if v == 91
        fuel_ids[i] = 2  # urban → short grass
    elseif v == 93
        fuel_ids[i] = 1  # agriculture → short grass
    elseif v <= 0 || v >= 98
        fuel_ids[i] = 256  # nodata (-9999), water, bare → nonburnable
    end
end

# Read elevation (meters)
elevation, gt = ArchGDAL.read(elev_tif) do ds
    data = Float64.(ArchGDAL.read(ArchGDAL.getband(ds, 1)))
    gt = ArchGDAL.getgeotransform(ds)
    data, gt
end

ncols, nrows = size(fuel_ids)
cellsize_m = gt[2]                       # 30 m
cellsize = cellsize_m / 0.3048           # convert to feet

# Compute slope and aspect from elevation (both in meters for consistent units)
slope, aspect = compute_slope_aspect(elevation, cellsize_m)

# Print summary
println("Grid:      $ncols x $nrows cells ($(round(ncols * cellsize_m / 1000, digits=1)) km x $(round(nrows * cellsize_m / 1000, digits=1)) km)")
println("Cell size: $(round(cellsize_m, digits=1)) m ($(round(cellsize, digits=1)) ft)")
println("Elevation: $(round(extrema(elevation)[1], digits=0))--$(round(extrema(elevation)[2], digits=0)) m")

# CRS projection for GeoMakie (native Albers → lon/lat display)
wkt = ArchGDAL.read(fbfm13_tif) do ds
    ArchGDAL.getproj(ds)
end
proj4_str = ArchGDAL.toPROJ4(ArchGDAL.importWKT(wkt))

# Native CRS coordinate arrays (cell centers in Albers meters)
xs = [gt[1] + (i - 0.5) * gt[2] for i in 1:ncols]
ys = [gt[4] + (j - 0.5) * gt[6] for j in 1:nrows]

# Fuel type distribution
for v in sort(unique(fuel_ids))
    n = count(==(v), fuel_ids)
    name = v == 256 ? "NB (water/bare)" : "FBFM$(lpad(v, 2, '0'))"
    println("  $name: $n cells ($(round(100n / length(fuel_ids), digits=1))%)")
end

Grid:      383 x 253 cells (11.5 km x 7.6 km)
Cell size: 30.0 m (98.4 ft)
Elevation: -9999.0--1857.0 m
  FBFM01: 64 cells (0.1%)
  FBFM02: 73165 cells (75.5%)
  FBFM04: 121 cells (0.1%)
  FBFM05: 4865 cells (5.0%)
  FBFM06: 1265 cells (1.3%)
  FBFM08: 13956 cells (14.4%)
  FBFM09: 1080 cells (1.1%)
  NB (water/bare): 2383 cells (2.5%)

Fire Perimeter from NIFC

We fetch the official Marshall Fire perimeter from the WFIGS Interagency Fire Perimeters REST API.

Fetch NIFC fire perimeter and reproject to grid coordinates

# Fetch perimeter GeoJSON from NIFC ArcGIS Feature Service
nifc_url = "https://services3.arcgis.com/T4QMspbfLg3qTGWY/arcgis/rest/services/" *
    "WFIGS_Interagency_Perimeters/FeatureServer/0/query?" *
    "where=poly_IncidentName%3D%27Marshall%27+AND+attr_POOState%3D%27US-CO%27" *
    "&outFields=poly_GISAcres,attr_InitialLatitude,attr_InitialLongitude" *
    "&f=geojson"

buf = IOBuffer()
Downloads.download(nifc_url, buf)
geojson = JSON3.read(String(take!(buf)))
feat = geojson.features[1]

println("NIFC Marshall Fire: $(round(feat.properties.poly_GISAcres, digits=0)) acres")

# Extract perimeter ring (lon, lat)
ring = feat.geometry.coordinates[1][1]
perim_lons = [c[1] for c in ring]
perim_lats = [c[2] for c in ring]

# Reproject perimeter from WGS84 to LANDFIRE CRS, then convert to grid coords
# Note: EPSG:4326 uses (lat, lon) axis order in GDAL 3+
target_crs = GeoFormatTypes.WellKnownText(GeoFormatTypes.CRS(), wkt)
source_crs = GeoFormatTypes.EPSG(4326)

perim_col = Float64[]
perim_row = Float64[]
perim_x_crs = Float64[]
perim_y_crs = Float64[]
for (lon, lat) in zip(perim_lons, perim_lats)
    x, y = ArchGDAL.reproject((lat, lon), source_crs, target_crs)
    push!(perim_x_crs, x)
    push!(perim_y_crs, y)
    push!(perim_col, (x - gt[1]) / gt[2] + 1)  # 1-based Julia indexing
    push!(perim_row, (y - gt[4]) / gt[6] + 1)
end

# Ignition point: documented origin near 288 Marshall Road, Boulder County
# (NIFC attr_InitialLongitude/attr_InitialLatitude are unreliable for this event)
ign_lon = -105.231
ign_lat = 39.955

ign_x, ign_y = ArchGDAL.reproject((ign_lat, ign_lon), source_crs, target_crs)
ign_col = clamp(round(Int, (ign_x - gt[1]) / gt[2] + 1), 1, ncols)
ign_row = clamp(round(Int, (ign_y - gt[4]) / gt[6] + 1), 1, nrows)

println("Perimeter: $(length(perim_lons)) vertices")
println("Ignition:  col=$ign_col, row=$ign_row ($ign_lon, $ign_lat)")

NIFC Marshall Fire: 6080.0 acres
Perimeter: 599 vertices
Ignition:  col=84, row=139 (-105.231, 39.955)

Static Overview

Four-panel summary of the real input data.

Code

dest_crs = "+proj=longlat +datum=WGS84"
fig = Figure(size = (900, 700))

# Panel 1: Elevation
ga1 = GeoAxis(fig[1, 1]; source = proj4_str, dest = dest_crs, title = "Elevation (m)")
contourf!(ga1, xs, ys, elevation; colormap = :terrain, levels = 12)
hidedecorations!(ga1)

# Panel 2: Fuel models
ga2 = GeoAxis(fig[1, 2]; source = proj4_str, dest = dest_crs, title = "LANDFIRE Fuel Models (FBFM13)")
heatmap!(ga2, xs, ys, Float64.(fuel_ids); colormap = [:lightgreen, :darkgreen, :sienna, :gray])
hidedecorations!(ga2)

# Panel 3: Slope
ga3 = GeoAxis(fig[2, 1]; source = proj4_str, dest = dest_crs, title = "Slope (degrees)")
contourf!(ga3, xs, ys, slope; colormap = :YlOrRd, levels = 8)
hidedecorations!(ga3)

# Panel 4: NIFC perimeter + ignition
ga4 = GeoAxis(fig[2, 2]; source = proj4_str, dest = dest_crs, title = "NIFC Fire Perimeter")
lines!(ga4, perim_x_crs, perim_y_crs; color = :red, linewidth = 2)
scatter!(ga4, [ign_x], [ign_y]; marker = :star5, markersize = 15, color = :orange)
hidedecorations!(ga4)

fig

Weather: KBDU ASOS Observations

Wind observations from Boulder Municipal Airport (KBDU) during the fire event on December 30, 2021. Data from the Iowa Environmental Mesonet.

The fire burned under extreme conditions: sustained winds of 35–50 mph from the WSW (250–270 degrees) with gusts to 75 mph. We compute an effective wind speed (average of sustained and gust) from the peak burning period (18:00–21:00 UTC), which better captures the influence of gusty conditions on fire spread.

KBDU ASOS wind observations (Dec 30, 2021)

# Real KBDU ASOS observations during peak fire period (18:00-22:00 UTC)
# Source: Iowa Environmental Mesonet ASOS download
# Columns: UTC hour, wind direction (deg), sustained speed (kt), gust (kt)
kbdu_obs = [
    17.25  290  21  44
    17.58  300  28  48
    17.92  300  16  42
    18.25  290  31  50
    18.58  270  28  50
    18.92  260  39  56
    19.25  250  44  61
    19.58  250  35  62
    19.92  240  35  63
    20.25  250  39  57
    20.58  250  41  55
    20.92  250  40  58
    21.25  250  43  65
    22.58  250  37  60
    22.92  250  35  55
    23.25  250  36  48
    23.58  250  28  47
    23.92  240  26  56
]

# Average over peak burning period (18:00-22:00 UTC)
peak = kbdu_obs[4:13, :]  # rows where hour is 18.25-20.92
avg_dir = round(sum(peak[:, 2]) / size(peak, 1), digits=0)
avg_kt = round(sum(peak[:, 3]) / size(peak, 1), digits=1)
avg_mph = round(avg_kt * 1.151, digits=0)
avg_gust_kt = round(sum(peak[:, 4]) / size(peak, 1), digits=1)
avg_gust_mph = round(avg_gust_kt * 1.151, digits=0)

# Effective wind: average of sustained and gust (common in fire behavior modeling)
eff_mph = round((avg_mph + avg_gust_mph) / 2, digits=0)

println("Peak burning period (18:00-21:00 UTC = 11AM-2PM MST):")
println("  Wind direction: $(avg_dir) degrees (WSW)")
println("  Sustained:      $(avg_kt) kt = $(avg_mph) mph")
println("  Gusts:          $(avg_gust_kt) kt = $(avg_gust_mph) mph")
println("  Effective:      $(eff_mph) mph (mean of sustained and gust)")

Peak burning period (18:00-21:00 UTC = 11AM-2PM MST):
  Wind direction: 256.0 degrees (WSW)
  Sustained:      37.5 kt = 43.0 mph
  Gusts:          57.7 kt = 66.0 mph
  Effective:      54.0 mph (mean of sustained and gust)

Run Simulation

We simulate 360 minutes (6 hours) using the real landscape and observed weather conditions. The effective wind speed (average of sustained and gust) captures the influence of gusts on fire spread better than sustained speed alone.

Code

fuel_table = create_standard_fuel_table(Float64)

# Weather from KBDU observations: effective wind speed during peak period
weather = ConstantWeather(
    wind_speed_mph = eff_mph,
    wind_direction = avg_dir,
    M1 = 0.03,    # 3% — extremely dry (December, chinook winds)
    M10 = 0.05,
    M100 = 0.08,
    MLH = 0.30,   # Low live moisture (dormant season)
    MLW = 0.60
)

state = FireState(ncols, nrows, cellsize)

# Ignite near the NIFC-reported origin (small circle to account for
# positional uncertainty in the ignition coordinates)
ignite_circle!(state, ign_col, ign_row, 3.0, 0.0)

# Capture snapshots for animation
n_frames = 40
dt_frame = 360.0 / n_frames  # 9 min per frame
frames_burned = [copy(state.burned)]
frames_toa = [copy(state.time_of_arrival)]
frame_times = [0.0]

for i in 1:n_frames
    t_start = (i - 1) * dt_frame
    t_stop = i * dt_frame
    simulate!(state, fuel_ids, fuel_table, weather,
        slope, aspect, t_start, t_stop;
        dt_initial = 0.5)
    push!(frames_burned, copy(state.burned))
    push!(frames_toa, copy(state.time_of_arrival))
    push!(frame_times, t_stop)
end

println("Simulation complete: $(length(frames_burned)) frames captured")
println("Burned area: $(round(get_burned_area_acres(state), digits=0)) acres")
println("NIFC reported: $(round(feat.properties.poly_GISAcres, digits=0)) acres")

Simulation complete: 41 frames captured
Burned area: 9391.0 acres
NIFC reported: 6080.0 acres

Animated Fire Progression

Five visual layers per frame: terrain contours, burned area (time-of-arrival), simulated fire perimeter, the NIFC real fire perimeter (cyan dashed), and wind particles.

Code

dest_crs = "+proj=longlat +datum=WGS84"
fig = Figure(size = (800, 500))

ga = GeoAxis(fig[1, 1];
    source = proj4_str, dest = dest_crs,
    title = "Marshall Fire — t = 0 min"
)
hidedecorations!(ga)

# Layer 1: terrain contours (static base)
contour!(ga, xs, ys, elevation; color = :gray70, linewidth = 0.4, levels = 10)

# Layer 2: time-of-arrival heatmap (dynamic)
toa_obs = Observable(fill(NaN, ncols, nrows))
hm = heatmap!(ga, xs, ys, toa_obs;
    colormap = :YlOrRd, colorrange = (0, 360), nan_color = :transparent)
Colorbar(fig[1, 2], hm; label = "Time (min)")

# Layer 3: simulated fire perimeter (dynamic)
perim_obs = Observable(Point2f[])
scatter!(ga, perim_obs; color = :orangered, markersize = 3, strokewidth = 0)

# Layer 4: NIFC real fire perimeter (static)
lines!(ga, perim_x_crs, perim_y_crs; color = :cyan, linestyle = :dash, linewidth = 1.5)

# Layer 5: wind particles (dynamic)
n_particles = 80
x_min, x_max = xs[1], xs[end]
y_lo, y_hi = min(ys[1], ys[end]), max(ys[1], ys[end])
wx_crs = [x_min + rand() * (x_max - x_min) for _ in 1:n_particles]
wy_crs = [y_lo + rand() * (y_hi - y_lo) for _ in 1:n_particles]
wind_obs = Observable([Point2f(wx_crs[k], wy_crs[k]) for k in eachindex(wx_crs)])
scatter!(ga, wind_obs; color = (:white, 0.6), markersize = 3, strokewidth = 0)

# Wind vector components (from WSW = 250 deg → blows ENE)
wind_rad = avg_dir * pi / 180
wind_dx = -sin(wind_rad)
wind_dy = -cos(wind_rad)

record(fig, "marshall_fire.gif", eachindex(frame_times); framerate = 4) do i
    ga.title = "Marshall Fire — t = $(round(Int, frame_times[i])) min"

    burned = frames_burned[i]
    toa = frames_toa[i]

    # Update TOA display (NaN = transparent for unburned cells)
    toa_mat = fill(NaN, ncols, nrows)
    for ix in 1:ncols, iy in 1:nrows
        if burned[ix, iy] && toa[ix, iy] >= 0
            toa_mat[ix, iy] = toa[ix, iy]
        end
    end
    toa_obs[] = toa_mat

    # Update fire perimeter (edge cells in CRS coordinates)
    pts = Point2f[]
    for ix in 1:ncols, iy in 1:nrows
        if burned[ix, iy]
            for (ddx, ddy) in ((1,0), (-1,0), (0,1), (0,-1))
                nx, ny = ix + ddx, iy + ddy
                if nx < 1 || nx > ncols || ny < 1 || ny > nrows || !burned[nx, ny]
                    push!(pts, Point2f(xs[ix], ys[iy]))
                    break
                end
            end
        end
    end
    perim_obs[] = pts

    # Update wind particles
    for k in eachindex(wx_crs)
        wx_crs[k] += rand(3:6) * cellsize_m * wind_dx
        wy_crs[k] += rand(3:6) * cellsize_m * wind_dy
        if wx_crs[k] > x_max || wx_crs[k] < x_min
            wx_crs[k] = x_min + rand() * cellsize_m * 5
            wy_crs[k] = y_lo + rand() * (y_hi - y_lo)
        end
        wy_crs[k] = clamp(wy_crs[k], y_lo, y_hi)
    end
    wind_obs[] = [Point2f(wx_crs[k], wy_crs[k]) for k in eachindex(wx_crs)]
end

Summary

Code

burned_acres = round(get_burned_area_acres(state), digits=0)
real_acres = round(feat.properties.poly_GISAcres, digits=0)

println("=" ^ 50)
println("Marshall Fire Simulation Summary")
println("=" ^ 50)
println("  Grid:             $ncols x $nrows cells ($(round(ncols * cellsize_m / 1000, digits=1)) km x $(round(nrows * cellsize_m / 1000, digits=1)) km)")
println("  Cell size:        $(round(cellsize_m, digits=0)) m ($(round(cellsize, digits=1)) ft)")
println("  Duration:         360 min (6 hours)")
println("  Wind:             $(round(Int, eff_mph)) mph from $(round(Int, avg_dir)) deg (KBDU ASOS effective)")
println("  Simulated area:   $burned_acres acres")
println("  NIFC reported:    $real_acres acres")

==================================================
Marshall Fire Simulation Summary
==================================================
  Grid:             383 x 253 cells (11.5 km x 7.6 km)
  Cell size:        30.0 m (98.4 ft)
  Duration:         360 min (6 hours)
  Wind:             54 mph from 256 deg (KBDU ASOS effective)
  Simulated area:   9391.0 acres
  NIFC reported:    6080.0 acres

Final State

Two-panel comparison of the simulated burned area vs the NIFC-reported perimeter and the time-of-arrival map.

Code

dest_crs = "+proj=longlat +datum=WGS84"
fig = Figure(size = (900, 450))

toa_final = copy(state.time_of_arrival)
toa_final[toa_final .< 0] .= NaN

# Panel 1: burned area with NIFC perimeter overlay
ga1 = GeoAxis(fig[1, 1]; source = proj4_str, dest = dest_crs, title = "Simulated vs Real Perimeter")
burned_float = Float64.(state.burned)
burned_float[burned_float .== 0] .= NaN
heatmap!(ga1, xs, ys, burned_float; colormap = :YlOrRd, nan_color = :transparent, colorrange = (0, 1))
lines!(ga1, perim_x_crs, perim_y_crs; color = :cyan, linestyle = :dash, linewidth = 2)
scatter!(ga1, [ign_x], [ign_y]; marker = :star5, markersize = 15, color = :white)
hidedecorations!(ga1)

# Panel 2: time-of-arrival
ga2 = GeoAxis(fig[1, 2]; source = proj4_str, dest = dest_crs, title = "Time of Arrival (min)")
hm = heatmap!(ga2, xs, ys, toa_final; colormap = :viridis, nan_color = :transparent)
Colorbar(fig[1, 3], hm; label = "min")
lines!(ga2, perim_x_crs, perim_y_crs; color = :cyan, linestyle = :dash, linewidth = 2)
hidedecorations!(ga2)

fig

Quantitative Validation

Comparing the simulated burned area against the NIFC perimeter cell-by-cell. We rasterize the observed perimeter polygon using a ray-casting point-in-polygon test, then compute standard binary classification metrics.

The primary metric is the Sørensen coefficient (also called the Dice score), defined as:

\[ \text{SC} = \frac{2 |A \cap B|}{|A| + |B|} \]

where \(A\) is the set of simulated burned cells and \(B\) is the set of observed burned cells. SC ranges from 0 (no overlap) to 1 (perfect agreement). It is the standard metric for fire spread model validation (Duff et al. 2018; Filippi et al. 2014) because it penalizes both commission errors (cells the model burned but the fire did not) and omission errors (cells the fire burned but the model missed), while being insensitive to the large number of true-negative (unburned) cells that dominate the domain.

We also report the Jaccard index (\(|A \cap B| / |A \cup B|\)), which is a monotonic transformation of SC and is sometimes preferred in spatial analysis.

Rasterize NIFC perimeter and compute overlap metrics

# Ray-casting point-in-polygon test
function point_in_polygon(px, py, poly_x, poly_y)
    n = length(poly_x)
    inside = false
    j = n
    for i in 1:n
        yi, yj = poly_y[i], poly_y[j]
        xi, xj = poly_x[i], poly_x[j]
        if ((yi > py) != (yj > py)) &&
           (px < (xj - xi) * (py - yi) / (yj - yi) + xi)
            inside = !inside
        end
        j = i
    end
    return inside
end

# Rasterize: for each grid cell, check if its center falls inside the NIFC polygon
observed_burned = falses(ncols, nrows)
for ix in 1:ncols, iy in 1:nrows
    if point_in_polygon(xs[ix], ys[iy], perim_x_crs, perim_y_crs)
        observed_burned[ix, iy] = true
    end
end

sim_burned = state.burned

# Binary classification
tp = count(sim_burned .& observed_burned)       # true positive
fp = count(sim_burned .& .!observed_burned)      # false positive (commission)
fn = count(.!sim_burned .& observed_burned)      # false negative (omission)
tn = count(.!sim_burned .& .!observed_burned)    # true negative

sorensen = 2tp / (2tp + fp + fn)
jaccard = tp / (tp + fp + fn)
commission_rate = fp / (tp + fp)   # fraction of simulated area that is false
omission_rate = fn / (tp + fn)     # fraction of observed area that was missed

obs_acres = count(observed_burned) * (cellsize_m^2 / 4046.86)
sim_acres = get_burned_area_acres(state)

println("=" ^ 50)
println("Validation Metrics")
println("=" ^ 50)
println("  Observed burned:    $(round(obs_acres, digits=0)) acres ($(count(observed_burned)) cells)")
println("  Simulated burned:   $(round(sim_acres, digits=0)) acres ($(count(sim_burned)) cells)")
println("  Area ratio (S/O):   $(round(sim_acres / obs_acres, digits=2))")
println()
println("  True positive:      $tp cells")
println("  False positive:     $fp cells (commission)")
println("  False negative:     $fn cells (omission)")
println()
println("  Sørensen coeff:     $(round(sorensen, digits=3))")
println("  Jaccard index:      $(round(jaccard, digits=3))")
println("  Commission rate:    $(round(100 * commission_rate, digits=1))%")
println("  Omission rate:      $(round(100 * omission_rate, digits=1))%")

==================================================
Validation Metrics
==================================================
  Observed burned:    6043.0 acres (27172 cells)
  Simulated burned:   9391.0 acres (42226 cells)
  Area ratio (S/O):   1.55

  True positive:      17872 cells
  False positive:     24354 cells (commission)
  False negative:     9300 cells (omission)

  Sørensen coeff:     0.515
  Jaccard index:      0.347
  Commission rate:    57.7%
  Omission rate:      34.2%

Spatial Agreement Map

Green = true positive (both agree the cell burned), red = commission error (simulated but not observed), blue = omission error (observed but not simulated).

Code

# Encode agreement: 0=TN, 1=TP (green), 2=FP/commission (red), 3=FN/omission (blue)
agreement = fill(NaN, ncols, nrows)
for ix in 1:ncols, iy in 1:nrows
    s, o = sim_burned[ix, iy], observed_burned[ix, iy]
    if s && o
        agreement[ix, iy] = 1.0  # TP
    elseif s && !o
        agreement[ix, iy] = 2.0  # FP (commission)
    elseif !s && o
        agreement[ix, iy] = 3.0  # FN (omission)
    end
end

dest_crs = "+proj=longlat +datum=WGS84"
fig = Figure(size = (700, 500))
ga = GeoAxis(fig[1, 1]; source = proj4_str, dest = dest_crs,
    title = "Spatial Agreement (Sørensen = $(round(sorensen, digits=3)))")
hidedecorations!(ga)

heatmap!(ga, xs, ys, agreement;
    colormap = [:green, :red, :royalblue],
    colorrange = (1, 3),
    nan_color = :transparent)

lines!(ga, perim_x_crs, perim_y_crs; color = :black, linestyle = :dash, linewidth = 1)

# Legend
Legend(fig[1, 2], [
    MarkerElement(marker = :rect, color = :green, markersize = 15),
    MarkerElement(marker = :rect, color = :red, markersize = 15),
    MarkerElement(marker = :rect, color = :royalblue, markersize = 15)
], ["True Positive", "Commission", "Omission"])

fig

Data Sources

LANDFIRE FBFM13 & Elevation: LANDFIRE Product Service v2.5.0, downloaded via Landfire.jl
Fire Perimeter: WFIGS Interagency Fire Perimeters, ArcGIS Feature Service
Wind Observations: KBDU ASOS, Boulder Municipal Airport, Dec 30, 2021
Fire Summary: NWS Boulder - Marshall Fire