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

30 m FBFM13 fuel models and elevation for the fire area come from LANDFIRE via Landfire.jl.

The clipped rasters are cached on this repository’s tutorial-data branch so the page renders reproducibly — the LANDFIRE Product Service is a live job queue and a failure there would otherwise abort the whole docs build. Set USE_LIVE_LANDFIRE=true to pull fresh data from the service instead.

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

const TUTORIAL_DATA = "https://raw.githubusercontent.com/RallypointOne/Elmfire.jl/tutorial-data"

"""
Fetch a cached tutorial input, or download it live from LANDFIRE.

`LF2025_FBFM13` covers only the SW and NW geoAreas, so the *latest* FBFM13 product
does not include Colorado — the version is pinned rather than taken as latest.
Elevation is published as `LF2020_Elev`.
"""
function landscape_layer(cached_name, layer)
    if get(ENV, "USE_LIVE_LANDFIRE", "false") == "true"
        products = Landfire.products(false; layer = layer)
        isempty(products) && error("No LANDFIRE product matching \"$layer\"")
        return get(Landfire.Dataset(products, aoi))
    end
    dest = joinpath(mktempdir(), cached_name)
    Downloads.download("$TUTORIAL_DATA/$cached_name", dest)
    return dest
end

fbfm13_tif = landscape_layer("marshall_fbfm13.tif", "LF2024_FBFM13")
elev_tif = landscape_layer("marshall_elev.tif", "LF2020_Elev")

println("Fuel models: $fbfm13_tif")
println("Elevation:   $elev_tif")
Fuel models: /tmp/jl_gxBQXf/marshall_fbfm13.tif
Elevation:   /tmp/jl_mqaLW7/marshall_elev.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

original_fuel_ids = copy(fuel_ids)

# Non-burnable display info (code, label, color) — used in both GIF and static plot
nb_info = [
    (91, "Urban → FBFM09",  RGBAf(0.45, 0.45, 0.45, 0.85)),
    (98, "Water",            RGBAf(0.26, 0.52, 0.80, 0.85)),
]

# Map LANDFIRE special codes:
#   91 (urban) → 9 (hardwood litter): slower-spreading surrogate for patchy WUI fuels
#   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] = 9  # urban → hardwood litter
    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:      382 x 253 cells (11.5 km x 7.6 km)
Cell size: 30.0 m (98.4 ft)
Elevation: 1618.0--1857.0 m
  FBFM01: 64 cells (0.1%)
  FBFM02: 43055 cells (44.5%)
  FBFM04: 119 cells (0.1%)
  FBFM05: 4860 cells (5.0%)
  FBFM06: 1266 cells (1.3%)
  FBFM08: 13974 cells (14.5%)
  FBFM09: 31181 cells (32.3%)
  NB (water/bare): 2127 cells (2.2%)

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 
# Perimeter GeoJSON, cached alongside the rasters. Set USE_LIVE_LANDFIRE=true to
# query the NIFC ArcGIS Feature Service directly.
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"

perimeter_url = get(ENV, "USE_LIVE_LANDFIRE", "false") == "true" ?
    nifc_url : "$TUTORIAL_DATA/marshall_perimeter.geojson"

buf = IOBuffer()
Downloads.download(perimeter_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 = "EPSG:3857"
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)

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")
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

Run Simulation

We simulate 360 minutes (6 hours) using the real landscape and observed weather conditions. We use the sustained wind speed (not gusts) as the input: the Rothermel model assumes steady-state surface fire spread, so sustained winds are more physically appropriate than the gust average used in some fire behavior tools.

Code 
fuel_table = create_standard_fuel_table(Float64)

# Weather from KBDU observations: sustained wind speed during peak period
weather = ConstantWeather(
    wind_speed_mph = avg_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,
        spread_rate_adj = 2.0,
        lb_cap = 4.0,
        dampening = SpreadRateDampeningConfig{Float64}(LINEAR_DAMPENING))
    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: 5233.0 acres
NIFC reported: 6080.0 acres

Animated Fire Progression

Five visual layers: terrain contours, burned area (time-of-arrival), simulated fire perimeter, the NIFC real fire perimeter (cyan dashed), and a static wind arrow field showing the constant WSW wind direction and relative speed.

Code 
dest_crs = "EPSG:3857"
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: non-burnable background (matches static plot)
for (code, label, color) in nb_info
    mask = Float64.(original_fuel_ids .== code)
    any(mask .> 0) || continue
    mask[mask .== 0] .= NaN
    heatmap!(ga, xs, ys, mask; colormap = [color, color], colorrange = (0.5, 1.5), nan_color = :transparent)
end

# Layer 3: 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 4: simulated fire perimeter (dynamic, red)
perim_obs = Observable(Point2f[])
scatter!(ga, perim_obs; color = :red, markersize = 3, strokewidth = 0)

# Layer 5: NIFC real fire perimeter (static, solid black)
lines!(ga, perim_x_crs, perim_y_crs; color = :black, linewidth = 1.5)

# Layer 6: ignition star with black outline
scatter!(ga, [ign_x], [ign_y]; marker = :star5, markersize = 15, color = :yellow,
    strokecolor = :black, strokewidth = 1.5)

# 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)

# Layer 7: wind field arrows (static — constant wind throughout simulation)
arrow_step = max(5, min(ncols, nrows) ÷ 8)
ax_idx = (arrow_step ÷ 2):arrow_step:ncols
ay_idx = (arrow_step ÷ 2):arrow_step:nrows
arrow_xs = vec([xs[i] for i in ax_idx, j in ay_idx])
arrow_ys = vec([ys[j] for i in ax_idx, j in ay_idx])
arrow_len = arrow_step * cellsize_m * 0.75
arrow_us = fill(wind_dx * arrow_len, length(arrow_xs))
arrow_vs = fill(wind_dy * arrow_len, length(arrow_ys))
arrows!(ga, arrow_xs, arrow_ys, arrow_us, arrow_vs;
    color = (:white, 0.75), arrowsize = 10, linewidth = 1.2)

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
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, avg_mph)) mph from $(round(Int, avg_dir)) deg (KBDU ASOS sustained)")
println("  Simulated area:   $burned_acres acres")
println("  NIFC reported:    $real_acres acres")
==================================================
Marshall Fire Simulation Summary
==================================================
  Grid:             382 x 253 cells (11.5 km x 7.6 km)
  Cell size:        30.0 m (98.4 ft)
  Duration:         360 min (6 hours)
  Wind:             43 mph from 256 deg (KBDU ASOS sustained)
  Simulated area:   5233.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 = "EPSG:3857"
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:    6036.0 acres (27141 cells)
  Simulated burned:   5233.0 acres (23531 cells)
  Area ratio (S/O):   0.87

  True positive:      15419 cells
  False positive:     8112 cells (commission)
  False negative:     11722 cells (omission)

  Sørensen coeff:     0.609
  Jaccard index:      0.437
  Commission rate:    34.5%
  Omission rate:      43.2%

Perimeter Comparison

Code 
# Fill interior holes in the burned mask before extracting the perimeter.
# Non-burnable cells (fuel code 256) inside the fire create unburned islands whose
# edges are picked up by the edge-detection, producing spurious interior "holes".
# Flood-fill from all grid edges to identify cells reachable from outside; any
# unburned cell not reachable from outside is interior and treated as burned.
function fill_holes(burned::AbstractMatrix{Bool})
    filled = copy(burned)
    outside = falses(size(burned))
    queue = Tuple{Int,Int}[]
    ncols, nrows = size(burned)
    for ix in 1:ncols, iy in 1:nrows
        if (ix == 1 || ix == ncols || iy == 1 || iy == nrows) && !burned[ix, iy]
            outside[ix, iy] = true
            push!(queue, (ix, iy))
        end
    end
    while !isempty(queue)
        ix, iy = pop!(queue)
        for (ddx, ddy) in ((1,0), (-1,0), (0,1), (0,-1))
            nx, ny = ix + ddx, iy + ddy
            if 1 <= nx <= ncols && 1 <= ny <= nrows && !burned[nx, ny] && !outside[nx, ny]
                outside[nx, ny] = true
                push!(queue, (nx, ny))
            end
        end
    end
    for ix in 1:ncols, iy in 1:nrows
        filled[ix, iy] = burned[ix, iy] || !outside[ix, iy]
    end
    return filled
end

sim_burned_filled = fill_holes(sim_burned)
sim_burned_float  = Float64.(sim_burned_filled)

dest_crs = "EPSG:3857"
fig = Figure(size = (1000, 560))
ga = GeoAxis(fig[1, 1]; source = proj4_str, dest = dest_crs,
    title = "Simulated vs NIFC Perimeter — Sørensen = $(round(sorensen, digits=3))")
hidedecorations!(ga)

# Base: terrain contours
contour!(ga, xs, ys, elevation; color = :gray80, linewidth = 0.4, levels = 10)

# Anderson 13 FBFM colors (standard LANDFIRE symbology approximation)
fbfm_colors = Dict(
     1 => RGBAf(0.99, 0.99, 0.60, 0.9),   # short grass — pale yellow
     2 => RGBAf(0.80, 0.90, 0.55, 0.9),   # timber grass — light yellow-green
     3 => RGBAf(0.55, 0.85, 0.20, 0.9),   # tall grass — bright green
     4 => RGBAf(0.90, 0.55, 0.15, 0.9),   # chaparral — orange
     5 => RGBAf(0.45, 0.70, 0.35, 0.9),   # brush — olive green
     6 => RGBAf(0.65, 0.50, 0.25, 0.9),   # dormant brush — tan-brown
     7 => RGBAf(0.50, 0.40, 0.20, 0.9),   # southern rough — dark brown
     8 => RGBAf(0.80, 0.72, 0.58, 0.9),   # compact timber litter — beige
     9 => RGBAf(0.65, 0.50, 0.35, 0.9),   # hardwood litter — warm brown
    10 => RGBAf(0.45, 0.33, 0.20, 0.9),   # timber litter — dark brown
    11 => RGBAf(0.78, 0.75, 0.68, 0.9),   # light slash — light gray
    12 => RGBAf(0.58, 0.55, 0.50, 0.9),   # medium slash — medium gray
    13 => RGBAf(0.38, 0.35, 0.30, 0.9),   # heavy slash — dark gray
)
fbfm_names = Dict(
     1 => "FBFM01: Short Grass",    2 => "FBFM02: Timber Grass",
     3 => "FBFM03: Tall Grass",     4 => "FBFM04: Chaparral",
     5 => "FBFM05: Brush",          6 => "FBFM06: Dormant Brush",
     7 => "FBFM07: Southern Rough", 8 => "FBFM08: Compact Timber Litter",
     9 => "FBFM09: Hardwood Litter",10 => "FBFM10: Timber Litter",
    11 => "FBFM11: Light Slash",    12 => "FBFM12: Medium Slash",
    13 => "FBFM13: Heavy Slash",
)

fbfm_legend_elems = LegendElement[]
fbfm_legend_labels = String[]
for id in 1:13
    mask = Float64.(fuel_ids .== id)
    any(mask .> 0) || continue
    mask[mask .== 0] .= NaN
    c = fbfm_colors[id]
    heatmap!(ga, xs, ys, mask; colormap = [c, c], colorrange = (0.5, 1.5), nan_color = :transparent)
    push!(fbfm_legend_elems, PolyElement(color = c, strokecolor = :transparent))
    push!(fbfm_legend_labels, fbfm_names[id])
end

# Nonburnable cells (fuel_ids == 256) — barriers that split the fire front
nb256_color = RGBAf(0.15, 0.15, 0.15, 0.85)
nb256_mask = Float64.(fuel_ids .== 256)
nb256_mask[nb256_mask .== 0] .= NaN
heatmap!(ga, xs, ys, nb256_mask; colormap = [nb256_color, nb256_color], colorrange = (0.5, 1.5), nan_color = :transparent)

# LANDFIRE non-burnable codes on top of FBFMs
nb_legend_elems = LegendElement[]
nb_legend_labels = String[]
for (code, label, color) in nb_info
    mask = Float64.(original_fuel_ids .== code)
    any(mask .> 0) || continue
    mask[mask .== 0] .= NaN
    heatmap!(ga, xs, ys, mask; colormap = [color, color], colorrange = (0.5, 1.5), nan_color = :transparent)
    push!(nb_legend_elems, PolyElement(color = color, strokecolor = :transparent))
    push!(nb_legend_labels, label)
end

# Simulated: red fill + contour outline (marching squares — no diagonal gaps)
burned_fill = copy(sim_burned_float)
burned_fill[burned_fill .== 0] .= NaN
heatmap!(ga, xs, ys, burned_fill;
    colormap = [RGBAf(0.7,0,0,0.3), RGBAf(0.7,0,0,0.3)], colorrange = (0.5, 1.5), nan_color = :transparent)
contour!(ga, xs, ys, sim_burned_float; levels = [0.5], color = RGBf(0.7, 0, 0), linewidth = 3)

# NIFC perimeter: solid black
lines!(ga, perim_x_crs, perim_y_crs; color = :black, linewidth = 2)

# Ignition star with black outline
scatter!(ga, [ign_x], [ign_y]; marker = :star5, markersize = 15, color = :yellow,
    strokecolor = :black, strokewidth = 1.5)

Legend(fig[1, 2],
    [fbfm_legend_elems..., nb_legend_elems...,
     PolyElement(color = RGBAf(0.7,0,0,0.3), strokecolor = RGBf(0.7,0,0), strokewidth = 1),
     LineElement(color = :black, linewidth = 2),
     MarkerElement(marker = :star5, color = :yellow, markersize = 12, strokecolor = :black, strokewidth = 1.5)],
    [fbfm_legend_labels..., nb_legend_labels..., "Simulated", "NIFC", "Ignition"];
    framevisible = false)

save("marshall_fire_perimeter.png", fig; px_per_unit = 3)
fig

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 = "EPSG:3857"
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