`geeViz.fireLib` - wildland fire modeling on Earth Engine¶

Earth Engine is lazy, tile-parallel and stateless. Fire spread is sequential. Most Earth Engine fire projects founder on that seam, so this package draws the line explicitly:

Tier	Where	What
1	Earth Engine	Everything pixel-wise: fuels, terrain, Rothermel
2	Earth Engine	Cost-distance propagation, no timestep loop
3	Outside	FSim, FlamMap, FARSITE, ELMFIRE

Requires Earth Engine. Every cell below needs ee initialized.

import ee
import geeViz.geeView as gv
import geeViz.fireLib as fl

Map = gv.Map
aoi = ee.Geometry.Rectangle([-120.9, 44.1, -120.7, 44.3])  # Ochoco NF
print('Earth Engine ready')

1. Fuels and terrain¶

LANDFIRE surface and canopy fuels live in the community catalog (projects/sat-io/...), not the official LANDFIRE/ namespace - the official one carries vegetation and fire-regime products but not the fuel models. An id that looks right and does not exist is an easy mistake to make quietly.

fuels = fl.landfire_fuels(region=aoi)
terrain = fl.terrain_layers('USGS/SRTMGL1_003', region=aoi)
print('fuels  :', fuels.bandNames().getInfo())
print('terrain:', terrain.bandNames().getInfo())

northness / eastness exist because raw aspect must never be fed to a model or a statistic. 359 degrees and 1 degree are adjacent on the ground and maximally distant numerically; averaging them gives 180, which points the wrong way.

Check fuel coverage BEFORE trusting any result¶

The parameter table is a verified subset of Scott & Burgan’s 40, not the whole set, and an unlisted model produces a masked pixel. That is invisible in the output: a landscape can look calm simply because part of it was unmodellable.

Deliberately not closed by inventing values - a wrong fuel load scales predicted spread directly and yields a confidently wrong fire.

cov = fl.fuel_coverage(fuels, aoi)
print(f"covered: {cov['covered_fraction']*100:.1f}% of pixels")
print(f"{cov['models_in_table']} models in table, "
      f"{cov['models_present']} present in this AOI")
for m, px, frac in cov['missing'][:5]:
    print(f'  missing model {m}: {px:,} px ({frac*100:.1f}%)')

2. Rothermel surface rate of spread¶

A pure pixel-wise function of fuel bed, slope, wind and moisture - no neighbors, no iteration, no state. This is what Earth Engine is genuinely ideal for.

Note wind_adjustment: the factor converting 20-ft wind to midflame wind. It moves the answer more than almost any other input (0.1 in dense timber, 0.6 in open grass), which is why it is an explicit argument rather than a buried constant.

scenarios = [
    ('calm, moist',      2.0,  0.12),
    ('moderate',         8.0,  0.06),
    ('wind-driven, dry', 25.0, 0.03),
]
for label, wind, mois in scenarios:
    ros = fl.rate_of_spread(fuels, terrain,
                            wind_speed_20ft=wind, moisture_1h=mois)
    s = ros.reduceRegion(ee.Reducer.percentile([50, 90]), aoi, 90,
                         maxPixels=1e10, bestEffort=True).getInfo()
    print(f"{label:20s} p50={s['ros_ft_min_p50']:7.1f}  "
          f"p90={s['ros_ft_min_p90']:7.1f} ft/min")

Those should be strictly ordered. If they are not, something upstream is wrong - an early version of this module returned exactly zero everywhere because the effective-heating term was written exp(-138 * sigma) instead of exp(-138 / sigma). With sigma around 2000 that underflows to zero, the heat sink collapses, and the whole landscape reads as fireproof.

3. Spread without a timestep loop¶

The instinct is to dilate the burned area a few hundred times. On Earth Engine each neighborhood operation expands the footprint a tile must fetch by one pixel, so a hundred nested ones need a 100-pixel halo on every tile.

cumulativeCost does it in one call. Least-accumulated-cost from a source is minimum travel time - the same quantity FlamMap’s MTT computes, and the solution to the Eikonal equation.

ros = fl.rate_of_spread(fuels, terrain,
                        wind_speed_20ft=12.0, moisture_1h=0.05)
ros_ms = fl.ros_metric(ros)
ignition = ee.Geometry.Point([-120.80, 44.20])

arrival = fl.travel_time(ros_ms, ignition, max_distance_m=8000)

for hours in (1, 3, 6):
    burned = arrival.lte(hours * 3600).selfMask()
    a = (ee.Image.pixelArea().updateMask(burned)
         .reduceRegion(ee.Reducer.sum(), aoi, 90,
                       maxPixels=1e10, bestEffort=True).getInfo())
    print(f"t={hours}h  burned {(a.get('area') or 0)/1e4:>10,.0f} ha")

The 100 frames are 100 thresholds of ONE image¶

No halo growth, no graph depth, trivially parallel.

frames = fl.isochrones(arrival, n_frames=24, total_seconds=24*3600)
print('frames:', frames.size().getInfo())
print('first :', ee.Image(frames.first()).get('t_hours').getInfo(), 'h')

Calibration - why `geodeticDistance=True` is the default¶

Run against a uniform spread rate the answer is analytic: reaching a point d metres away at r m/s must take d/r seconds. Measured at 44.2 deg N:

setting	east	north	diagonal
`geodeticDistance=False`	1.386	0.996	1.261
`geodeticDistance=True`	0.997	0.993	1.053

1/cos(44.2 deg) = 1.394. On an EPSG:4326 image a degree of longitude was being treated as a degree of latitude, so fire spread 39% too slowly east-west and correctly north-south - a directional error that looks entirely plausible on a map, and doubles at 60 deg N.

The same test settled the cost units: doubling the spread rate halves arrival times exactly, so cost accumulates per metre.

import math

lon, lat, D = -120.80, 44.20, 2000.0
pt = ee.Geometry.Point([lon, lat])
uniform = ee.Image.constant(1.0).rename('ros_m_s').clip(pt.buffer(6000))
arr_u = fl.travel_time(uniform, pt, max_distance_m=5000)

m_lat = 111320.0
m_lon = 111320.0 * math.cos(math.radians(lat))
for name, dx, dy in [('east', 1, 0), ('north', 0, 1),
                     ('NE', 0.7071, 0.7071)]:
    tgt = ee.Geometry.Point([lon + D*dx/m_lon, lat + D*dy/m_lat])
    o = arr_u.reduceRegion(ee.Reducer.first(), tgt.buffer(45), 30,
                           maxPixels=1e9).getInfo()
    v = next((x for x in o.values() if x is not None), None)
    print(f'  {name:6s} ratio to analytic: {v/D:5.3f}' if v
          else f'  {name}: None')

4. Wind shifts - iterate coarsely, not finely¶

cumulativeCost assigns cost per pixel, not per edge, so it cannot express ‘cheap downwind, expensive upwind’. A single call is directionally neutral - no elliptical head-fire elongation.

The workaround is to chain one call per wind period rather than one per timestep. Twenty calls for a five-day fire is comfortable, and it captures the wind shifts that drive the large runs. Within-block anisotropy is still lost: fine for planning and risk, not fine for operational head-fire prediction.

blocks = [
    {'wind_speed':  6.0, 'moisture_1h': 0.08, 'duration_s': 6*3600},
    {'wind_speed': 18.0, 'moisture_1h': 0.04, 'duration_s': 6*3600},
    {'wind_speed': 10.0, 'moisture_1h': 0.06, 'duration_s': 6*3600},
]

def burned_ha(img):
    a = (ee.Image.pixelArea().updateMask(img.selfMask())
         .reduceRegion(ee.Reducer.sum(), aoi, 90,
                       maxPixels=1e10, bestEffort=True).getInfo())
    return (a.get('area') or 0) / 1e4

for n in (1, 2, 3):
    p = fl.spread_with_wind_blocks(fuels, terrain, blocks[:n], ignition)
    print(f'after {n} block(s) ({n*6:>2}h): {burned_ha(p):>10,.0f} ha')

The third block adds nothing, and that is correct rather than broken: the 18 mph block already pushed the fire across most of this AOI, and fuels is clipped to it, so there is no burnable ground left to reach. Widen the AOI to see the third block do work.

One implementation note worth carrying into your own code: unmask(0) before the union is load-bearing. An image from paint() is masked everywhere except the painted geometry, and Earth Engine intersects masks on a binary operation - so perim.Or(grown) inherits the ignition point’s one-pixel mask and the union collapses back to the ignition. Measured while building this: grown covered 1,650 ha and the union of it returned 0.6 ha. A valid image, a plausible small number, and completely wrong.

5. Map it¶

Arrival time as an isochrone surface, with the fuels underneath.

Map.clearMap()
Map.addLayer(fuels.select('FBFM40'), {'autoViz': True},
             'LANDFIRE FBFM40', False)
Map.addLayer(terrain.select('slope'),
             {'min': 0, 'max': 45, 'palette': ['ffffff', '444444']},
             'Slope (deg)', False)
Map.addLayer(ros, {'min': 0, 'max': 200,
                   'palette': ['2b83ba', 'ffffbf', 'd7191c']},
             'Rothermel ROS (ft/min)', True)
Map.addLayer(arrival.divide(3600),
             {'min': 0, 'max': 12,
              'palette': ['d7191c', 'fdae61', 'ffffbf', '2b83ba']},
             'Arrival time (hours)', True)
Map.addLayer(ee.FeatureCollection([ee.Feature(ignition)]),
             {'strokeColor': '00FF00'}, 'Ignition', True)
Map.centerObject(aoi, 11)
Map.view()

What this package will not do¶

Stated so it can be designed around rather than discovered late:

No true elliptical head-fire spread. Cost is per pixel, not per edge. Chained wind blocks approximate direction; they do not reproduce Huygens wavelets.
No coupled fire-atmosphere behavior. Plume dynamics and downdraft-driven runs are WRF-Fire and QUIC-Fire territory.
No spotting. Ember transport is stochastic; cost-distance cannot express it.
No replacement for a project-level FSim run. For published burn probability, use USDA/WRC/v0 - it is FSim output already computed for CONUS, AK and HI at 30 m.

wrc = ee.ImageCollection(fl.RISK_ASSET)
print('WRC bands:', ee.Image(wrc.first()).bandNames().getInfo())
print('  BP=burn probability, CFL=conditional flame length,')
print('  FLEP4/8=flame length exceedance, WHP=hazard potential')

geeViz.fireLib - wildland fire modeling on Earth Engine¶