Use Cases Archives

El Niño 2026, atmospheric rivers, and California’s burn scars: one SQL engine, two data models

Posted on September 24, 2026October 4, 2026 by Jia Yu

Governor Gavin Newsom proclaimed a state of emergency on 21 September 2026 to prepare California for El Niño 2026. It names recent wildfire areas as particularly susceptible to mud and debris flows. Above Altadena, the Eaton fire left 42 basins that the US Geological Survey (USGS) rates high hazard for debris flows. USGS computes the rating for a 15-minute burst of 40 mm/h of rain (1.57 in/hour), and an atmospheric river can deliver that kind of deluge. The rating covers the basins on the hillside. The homes and roads downslope of the basins fall outside the rating, and the queries here count them.

WherobotsDB does the counting in SQL, reading every dataset in place. A join to the USGS basins returns 7,713 mapped building footprints inside the high-hazard basins or within 500 m (1,640 ft) of them. A raster vector join over elevation returns 388 km (241 mi) of mapped road and path.

Start with elevation, buildings and roads in the catalog

All four datasets sit in wherobots_open_data: the ASTER GDEM 30 m (98 ft) elevation model, Overture’s buildings and roads, and Regrid’s parcel sample. The companion notebook fetches the CAL FIRE perimeter and USGS basins. It also registers every shared view, so each query below runs as written, with no setup of its own. Each elevation row stores a footprint polygon, so WherobotsDB selects rows by location before it reads any pixel. Ten rasters cover the Eaton area.

Filter on the bounding box before the spatial test

Overture features carry a bbox struct beside the geometry. With a range filter on it before ST_Intersects, the first building count finished in 22.8 seconds. Without it, the run was cancelled at 7% after three minutes. The range check on four numbers is cheap, so the costly geometry test runs only on the rows that pass.

Show the code

-- scar holds the Eaton perimeter; the notebook fetches it from CAL FIRE and registers it as a view
WITH aoi AS (
  SELECT ST_Transform(
           ST_Buffer(ST_Transform(g, 'EPSG:4326', 'EPSG:32611'), 3000),
           'EPSG:32611', 'EPSG:4326') AS g
  FROM scar
),
bld AS (
  SELECT t.id, t.geometry AS geom
  FROM wherobots_open_data.overture_maps_foundation.buildings_building t
  WHERE t.bbox.xmin BETWEEN -118.20 AND -117.97   -- range scan first
    AND t.bbox.ymin BETWEEN  34.13  AND  34.27
)
SELECT count(*) FROM bld b JOIN aoi ON ST_Intersects(b.geom, aoi.g);

118,746 footprints pass the bounding box, 83,930 fall within 3 km (1.9 mi) of the perimeter, and 10,652 sit inside the scar.

Each of two questions below gets its own raster vector join. For buildings the answer is one elevation per footprint, so the first join reads the raster under every footprint. For roads, however, the answer is a length. So the second join turns the raster into a zone polygon first, then clips each road segment against it.

The two joins used below, drawn on a toy elevation grid with a 230 m (755 ft) line.

Buildings: read the elevation under each footprint

Which building footprints sit below the burn scar’s lowest point?

15,233 of the 83,930 footprints sit below that point, 223 m (732 ft). RS_ZonalStats runs the raster vector join per building, so every footprint gets one elevation value of its own. That value becomes one more column beside the building’s other attributes.

Show the code

-- scar, dem and inaoi (the 83,930 footprints within 3 km) are views from the notebook's setup cell
WITH scar_floor AS (
  SELECT min(RS_ZonalStats(d.rast, s.g, 1, 'min', true)) AS m FROM dem d, scar s
),
be AS (
  -- pool pixels across every raster a footprint touches: total value over pixel count
  SELECT i.id,
         sum(RS_ZonalStats(d.rast, i.geom, 1, 'sum', true))
           / NULLIF(sum(RS_ZonalStats(d.rast, i.geom, 1, 'count', true)), 0) AS elev
  FROM inaoi i JOIN dem d ON RS_Intersects(d.rast, i.geom)
  GROUP BY i.id
)
SELECT count(*) AS below_scar_floor
FROM be, scar_floor
WHERE be.elev < scar_floor.m;

Footprints that straddle two rasters have their pixels pooled before averaging, since the rasters meet under them. 380 of the 384 stay on the same side of the line. A 30 m (98 ft) pixel is wider than many houses, but every footprint still got an elevation value from the query.

Mean elevation under each of the 83,930 footprints within 3 km (1.9 mi) of the perimeter, in 20 m (66 ft) bands. Blue marks the buildings below the scar’s lowest point.

Roads: a raster vector join through a low-ground zone

Why do roads need a different join than buildings?

Roads need a different join because length, not a single point, decides the answer. A building gets one average elevation, and that single value places it above or below the line. A 1 km (0.6 mi) road dips below the line for a short stretch and stays above it for the rest. An average across the whole road would place it on one side, and the low stretch would go uncounted. The query takes three steps instead. RS_MapAlgebra marks every elevation pixel below 223 m (732 ft). RS_Polygonize turns the marked pixels into polygons, and ST_Union_Aggr merges them into one shape of low ground. ST_Intersection then cuts each road at the polygon’s edge, and the query adds up the length inside it. WherobotsDB builds the polygon once and tests every road against that same shape.

Show the code

-- aoi, dem and roads are views from the notebook's setup cell
WITH masked AS (
  SELECT RS_MapAlgebra(rast, 'D', 'out[0] = rast[0] < 223 ? 1 : 0;') AS m FROM dem
),
polys AS (
  SELECT t.p.geom AS g
  FROM masked LATERAL VIEW explode(RS_Polygonize(m, 1)) t AS p
  WHERE t.p.value = 1.0
),
low AS (
  SELECT ST_Intersection(u.g, a.g) AS g
  FROM (SELECT ST_Union_Aggr(ST_MakeValid(g)) AS g FROM polys) u, aoi a
)
SELECT r.class,
       sum(ST_Length(ST_Transform(ST_Intersection(r.geom, low.g),
                                  'EPSG:4326', 'EPSG:32611'))) / 1000 AS km
FROM roads r JOIN low ON ST_Intersects(r.geom, low.g)
GROUP BY r.class
ORDER BY km DESC;

Kilometres of Overture road and path inside the 16.8 km² (6.5 sq mi) of ground below the scar’s lowest point, within 3 km (1.9 mi) of the perimeter.

The low ground covers 16.8 km² (6.5 sq mi) and holds 388 km (241 mi) of road and path within 3 km (1.9 mi) of the scar. Of that, 305 km (190 mi) is drivable and 24 km (15 mi) is motorway. Lengths follow each mapped segment, so both carriageways of every divided road are counted in full.

Joining to what USGS already modelled

An elevation line is only a screen. USGS publishes the basins as a live feature service, so they join in as one more vector layer. Each building counts once per distance, however many basins it touches.

Show the code

-- high holds the 42 dissolved high-hazard basins; the notebook fetches them from USGS
WITH z AS (
  SELECT g AS inside,
         ST_Transform(ST_Buffer(ST_Transform(g, 'EPSG:4326', 'EPSG:32611'), 500),  'EPSG:32611', 'EPSG:4326') AS b500,
         ST_Transform(ST_Buffer(ST_Transform(g, 'EPSG:4326', 'EPSG:32611'), 1000), 'EPSG:32611', 'EPSG:4326') AS b1k
  FROM high
)
SELECT sum(CASE WHEN ST_Intersects(b.geom, z.inside) THEN 1 ELSE 0 END) AS inside,
       sum(CASE WHEN ST_Intersects(b.geom, z.b500)   THEN 1 ELSE 0 END) AS inside_or_500m,
       sum(CASE WHEN ST_Intersects(b.geom, z.b1k)    THEN 1 ELSE 0 END) AS inside_or_1km
FROM bld b, z
WHERE ST_Intersects(b.geom, z.b1k);

Overture building footprints inside, or within 500 m (1,640 ft) and 1 km (0.6 mi) of, the 42 Eaton basins that USGS rates high combined hazard.

696 footprints sit inside a high-hazard basin, 7,713 inside or within 500 m (1,640 ft), and 16,748 within 1 km (0.6 mi). The map below plots the 7,713 on Sentinel-2 imagery from 5 August 2026. They form a solid band along the scar’s lower edge, where the basins drain onto the valley floor. Nineteen months on, the scar is still bare.

Raster vector join result: 7,713 building footprints near 42 high-hazard debris-flow basins below the Eaton burn scar in Altadena — The 7,713 footprints (cyan) inside or within 500 m (1,640 ft) of the 42 high-hazard basins (red), on Sentinel-2 imagery from 5 August 2026.

Joining parcels to put a value on the screen

Building counts measure how many structures sit in each screen, not what they are worth. Regrid’s parcel sample in wherobots_open_data.partner_samples adds use codes, acreage and assessed values, and a join to the exposed footprints carries those onto both screens. The join runs in two steps so nothing is counted twice: first the distinct parcel and footprint pairs, keyed on Regrid’s stable ll_uuid, then one row per parcel for its use code, acreage and assessed value. A footprint that crosses a parcel line counts once in the total and once in each use it touches, so the rows sum slightly above the total.

Show the code

-- exposed holds the footprints from either screen; parcels is Regrid's sample with the same bbox filter as bld
WITH pairs AS (
  SELECT DISTINCT p.ll_uuid, p.usedesc, e.id AS bid
  FROM parcels p JOIN exposed e ON ST_Intersects(p.geometry, e.geom)
),
parcel AS (  -- one row per parcel, so acres and value are counted once
  SELECT ll_uuid, first(usedesc) AS usedesc, first(ll_gisacre) AS acres, first(parval) AS parval
  FROM parcels WHERE ll_uuid IN (SELECT ll_uuid FROM pairs) GROUP BY ll_uuid
),
by_use AS (
  SELECT usedesc, count(*) AS parcels, round(sum(acres)) AS acres, round(sum(parval) / 1e6, 1) AS assessed_musd
  FROM parcel GROUP BY usedesc
),
bld_use AS (SELECT usedesc, count(distinct bid) AS buildings FROM pairs GROUP BY usedesc)
SELECT u.usedesc, u.parcels, b.buildings, u.acres, u.assessed_musd
FROM by_use u JOIN bld_use b USING (usedesc) ORDER BY u.parcels DESC;

Inside or within 500 m (1,640 ft) of a high-hazard basin: 5,934 parcels, 7,709 footprints

Use	Parcels	Buildings	Acres	Assessed value
Residential single	5,686	7,157	2,910	$4.69 B
Residential two units	72	186	34	$51 M
Government parcel	67	175	944	$4 M
Residential three units	29	86	14	$24 M
Utility pumping plants	15	33	25	$0.1 M
Apartments, five or more	14	82	10	$13 M
Miscellaneous	10	48	3,383	$0
Motion picture, radio, TV	8	52	868	$12 M

Below the scar’s lowest point, 223 m (732 ft): 10,513 parcels, 15,227 footprints

Use	Parcels	Buildings	Acres	Assessed value
Residential single	9,032	12,791	3,013	$8.49 B
Residential two units	373	988	77	$210 M
Apartments, five or more	172	498	72	$499 M
Commercial stores	129	187	50	$222 M
Commercial office buildings	105	184	62	$511 M
Residential three units	86	233	17	$56 M
Light manufacturing	84	123	27	$119 M
Residential four units	81	210	18	$66 M

The 7,713 footprints inside or within 500 m (1,640 ft) of a high-hazard basin sit on 5,934 parcels with a combined assessed value of $4.8 billion. Single-family homes account for $4.7 billion of it. The downhill screen is larger: the 15,233 footprints below the scar’s lowest point sit on 10,513 parcels assessed at $11.1 billion, and here the mix changes. Single-family homes are $8.5 billion of that total, and the remaining $2.6 billion is the commercial and apartment stock along the valley floor, led by $511 million in office buildings, $499 million in apartments and $222 million in stores. The two screens overlap, so the totals do not add. Every value comes from the assessor’s 2023 tax roll, so it describes the neighbourhood before the fire and before any rebuilding. Proposition 13 holds California assessed values well below market value, so these totals are a floor.

RS_Value reads the elevation under a point, so a cross-section is a line of generated points.

Show the code

-- a point every 0.00025 degrees, about 28 m, along 118.125° W, and the elevation under each
WITH pts AS (
  SELECT i, ST_SetSRID(ST_Point(-118.125, 34.265 - i * 0.00025), 4326) AS pt
  FROM (SELECT explode(sequence(0, 460)) AS i)
)
SELECT p.i, ST_Y(p.pt) AS lat, RS_Value(d.rast, p.pt, 1) AS elev
FROM pts p
JOIN wherobots_open_data.aster_gdem.v3_30m d ON ST_Intersects(d.footprint, p.pt)
ORDER BY p.i;

Terrain sampled every 28 m (92 ft) along 118.125° W, from the ridge north of the scar to the valley floor. Red bands mark the burn scar and the basins USGS rates high; cyan ticks mark buildings within 60 m (197 ft) of the line.

The line climbs to 1,577 m (5,174 ft) on the ridge. Then it crosses 2.4 km (1.5 mi) of high-hazard basin and 4.5 km (2.8 mi) of scar, down to 248 m (814 ft). 1,112 footprints lie within 60 m (197 ft) of it. The first of those buildings already sit inside the lower end of the high-hazard stretch.

Scope and limits of the exposure screen

These counts are an exposure screen. Neither they nor the USGS assessment predict where a debris flow travels or what it buries. USGS notes that most come within two years of a fire. So this winter falls near the end of that window, while the burned slopes are still recovering.

Three limits carry over. The elevation model is a surface model, so rooftops and treetops count as ground. Proximity is straight-line distance rather than flow path. Finally, the perimeter is simplified to 47 vertices, which changes its area by 0.4%. Flow routing would turn distance bands into runout paths. However, WherobotsDB has no built-in flow-direction function, so that takes a custom function over the rasters. It is the natural next step.

From one burn scar to the western United States

This screen covers one burn scar: the Eaton fire above Altadena. The same queries extend to every recent burn scar in California, and across the western United States, where USGS publishes debris-flow hazard assessments after major fires. Overture buildings and roads and the ASTER elevation model in wherobots_open_data already cover that whole region. Each new fire adds a perimeter and a set of basins to the same joins, once the study-area bounding box and UTM zone are updated to match. Governor Newsom declared the emergency ahead of El Niño rain this winter. The open question for every recently burned community is how many homes and roads sit below a high-hazard basin when that rain arrives.

One engine processed buildings, roads, elevation, and USGS hazard polygons, as vector and raster, with every file left where it sits. An AI agent grounded only in text, documents, databases, and the internet holds none of those layers, so it cannot say what sits below a basin. WherobotsDB supplies that physical-world context by joining them in place.

Imagery: Copernicus Sentinel data 2026. Elevation: ASTER GDEM is a product of METI and NASA. Buildings and roads: Overture Maps Foundation. Parcels: Regrid parcel sample. Fire perimeter: CAL FIRE FRAP, CC BY. Hazard basins: USGS Landslide Hazards Program, public domain.

Key takeaways

A raster vector join matches raster pixels to vector geometry by location. It is the case of a spatial join where one side is raster (elevation) and the other is vector (building and road geometry).
WherobotsDB runs both joins in SQL, reading Overture buildings and roads, the ASTER GDEM elevation model, Regrid parcels and the USGS hazard basins in place, as vector and raster, with no data movement.
7,713 building footprints sit inside the 42 Eaton high-hazard basins or within 500 m (1,640 ft). 696 sit inside a basin. 16,748 sit within 1 km (0.6 mi).
388 km (241 mi) of road and path sit in the 16.8 km² (6.5 sq mi) of ground below the scar’s lowest point. 305 km (190 mi) is drivable and 24 km (15 mi) is motorway.
Joined to Regrid parcels, the footprints near the basins sit on 5,934 parcels assessed at $4.8 billion on the 2023 tax roll.
A bounding-box filter before the spatial test cut the first building count to 22.8 seconds. Without it, the run was cancelled at 7% after three minutes.
The same queries extend to every recent burn scar in California and across the western United States, where USGS publishes debris-flow hazard assessments after major fires.

Get Started with Wherobots

Try Now

Measuring the Strait of Hormuz shutdown with Sentinel-2, WherobotsDB and Overture Maps

Posted on September 16, 2026October 4, 2026 by Jia Yu

Iran closed the Strait of Hormuz to commercial shipping on 28 February 2026. To measure what stopped moving, WherobotsDB read 719 Sentinel-2 scenes, 151 GB in total, straight from the public archive and fetched only about 1.5 GB covering three areas of interest. A near-infrared threshold found the hulls with no trained model, and a join against Overture Maps land polygons removed rocks and coastline. Traffic through the surveyed corridor fell 95% while loading and anchorage traffic held at 2025 levels, and the whole analysis cost about $63.

Sentinel-2 imagery is free and covers the planet every five days. Pointing the same job at another chokepoint takes two edited arguments.

What the satellites saw

We counted vessel hulls in Sentinel-2 imagery over three areas from January 2025 to August 2026 using a pixel-based reflectance method, while accounting for the sun’s angle over time:

area	2025	2026	change
Hormuz transit corridor	10.3 vessels	0.5	−95%
Fujairah anchorage, outside the strait	63.5	70.0	+10%
Kharg Island export terminal, inside the Gulf	7.3	9.9	+35%

Those figures are per satellite pass, averaged over the passes clear enough to cover the whole area.

We used imagery because the usual source cannot be trusted here. Ship traffic is normally tracked through AIS, which is the transponder every large vessel carries. In this crisis those transponders are being switched off or falsified. Near Fujairah and Khor Fakkan, roughly 470 vessels broadcast garbled or impossible positions inside a single 24-hour window. Windward counted 146 of 167 vessels in the strait area running dark on 5 May, and IMF PortWatch warns that its own numbers understate traffic. The optical imagery provides a signal (the hull) that is always on.

Vessel detections per satellite pass across the three areas, January 2025 to August 2026. The corridor line falls to near zero after the closure while the anchorage and the loading terminal return to 2025 levels.

How to count ships from space

Water absorbs almost all near-infrared light, so in Sentinel-2’s B08 band the sea is nearly black. Painted steel reflects it strongly. Over the Fujairah anchorage on August 7th, 2026 the water sat at 0.09 reflectance and the brightest pixel in the area read 0.87, so a vessel stands out by a wide margin. This means you don’t need a machine-learning model to identify a ship.

Wherobots queries the Sentinel-2 files directly in the public archive. The raster data source opens each scene in place and fetches only the necessary pixels.

The code below is an excerpt from the job script but the whole pipeline is available as a runnable notebook that opens in Wherobots Cloud, with every constant in one cell so you can point it at your own area of interest.

Show the code

AOI = ("POLYGON((56.42 25.10, 56.62 25.10, 56.62 25.30, "
       "56.42 25.30, 56.42 25.10))")
AOI_SQL = f"ST_SetSRID(ST_GeomFromWKT('{AOI}'), 4326)"

(sedona.read.format("raster")
    .option("tileWidth", "512").option("tileHeight", "512")
    .option("autoRescale", "true")      # off by default: without it you get raw DN
    .load(paths)
    .createOrReplaceTempView("tiles"))

# Background statistics on water only. The median centre ignores a bright
# minority of pixels. Pooling the per-tile standard deviations by their median
# discards a tile spoiled by cloud or land, but bright pixels inside a tile
# still inflate that tile's spread, so a crowded anchorage lifts the cutoff
# somewhat. That bias is conservative and falls on both years alike.
r = sedona.sql(f"""
    WITH s AS (
        SELECT zs.median AS med, zs.stddev AS sd
        FROM (SELECT RS_ZonalStatsAll(rast, {WATER_G}) AS zs FROM tiles)
        WHERE zs.count > 0
    )
    SELECT percentile_approx(med, 0.5) AS med, percentile_approx(sd, 0.5) AS sd
    FROM s
""").collect()[0]

threshold = r["med"] + max(SIGMA_K * r["sd"], 0.03)

2026-08-07  cov=100.0% cloud= 2.67%  med=0.0915 sd=0.0097 thr=0.1687  \
brightpx=   5636  vessels=  84  >=200m= 39

RS_ZonalStatsAll reports statistics for the pixels inside a shape. Using the median prevents a small number of bright pixels from raising the threshold; an average would be pulled upward by every ship in the frame. Adding a multiple of the spread gives a cutoff for that scene. The spread is pooled the same way, by taking the median of each tile’s standard deviation, which discards a tile spoiled by clouds. It is only a partial defense: ships inside a tile raise that tile’s spread, and with it the cutoff, so the count runs slightly conservative in a crowded anchorage. The sea is brighter in July than in January, so the cutoff has to move with it, which we account for in the analysis.

A short Python function applies that cutoff, and RS_Polygonize traces the outline of every bright patch.

Show the code

@sedona_vectorized_udf(return_type=RasterType())
def threshold_udf(rast: SedonaRaster) -> SedonaRaster:
    a = rast.as_numpy_masked()[0]
    out = np.zeros(a.shape, dtype=np.float64)
    with np.errstate(invalid="ignore"):
        out[np.isfinite(a) & (a >= threshold)] = 1.0
    return rast.with_bands(out[np.newaxis])


# with_bands() needs real pixels; the reader returns out-db references.
(sedona.sql("SELECT RS_AsInDB(rast) AS rast FROM tiles")
    .withColumn("mask", threshold_udf("rast"))
    .createOrReplaceTempView("masked"))

# RS_Polygonize returns every connected region in ONE array, so running it on a
# whole scene OOMs an executor -- per tile it is cheap and parallel.
sedona.sql(f"""
    SELECT t.p.geom AS g
    FROM masked LATERAL VIEW explode(RS_Polygonize(mask, 1)) t AS p
    WHERE t.p.value = 1.0
      AND ST_Area(t.p.geom) BETWEEN 2000 AND 200000
""").createOrReplaceTempView("raw_blobs")


# A vessel lying on a tile seam is polygonized twice, once either side of the
# cut. Union everything, then split it back apart: the halves become one
# polygon, while vessels that never touched stay separate. The area floor above
# runs before this merge, which has a cost the last section spells out.
sedona.sql("""
    WITH dissolved AS (
        SELECT explode(ST_Dump(ST_Union_Aggr(ST_MakeValid(g)))) AS g
        FROM raw_blobs
    ),
    sized AS (
        SELECT g, ST_Area(g) AS area_m2,
               2 * ST_MinimumBoundingRadius(g).radius AS length_m,
               ST_MinimumWidth(g) AS width_m
        FROM dissolved
        WHERE ST_Area(g) BETWEEN 2000 AND 200000
    )
    SELECT g AS geometry, area_m2, length_m, width_m
    FROM sized
    WHERE width_m >= 20
      AND length_m / width_m BETWEEN 2 AND 15
""").createOrReplaceTempView("candidates")

Sentinel-2 true color over the Fujairah anchorage on August 7th, 2026. A detected 280 m hull covers about 12,600 m², which is roughly two football pitches and matches the deck of a Suezmax tanker.

Ruling out rocks, land and glint

Rocky islets also reflect near-infrared light. The Quoin Islands sit directly in the shipping lane, and in an early run three of our first four detections turned out to be those rocks.

Wherobots holds imagery and map data in one engine, so checking detections against an independent source takes one query. Overture Maps publishes global land polygons. They already sit in the Wherobots open data catalog, and one join drops any detection that touches land.

Show the code

sedona.sql(f"""
    SELECT ST_Buffer(ST_Transform(geometry, 'EPSG:4326', 'EPSG:{epsg}'), 60) AS g
    FROM wherobots_open_data.overture_maps_foundation.base_land
    WHERE ST_Intersects(geometry, {AOI_SQL})
""").createOrReplaceTempView("land")

detections = sedona.sql("""
    SELECT d.geometry, d.area_m2, d.length_m
    FROM candidates d
    LEFT ANTI JOIN land l ON ST_Intersects(d.geometry, l.g)
""")

Two further checks run alongside it. Sun glint can brighten patches of open water, so we keep only detections shaped like ships, which filters out glint using the near-infrared band alone. A hull is typically four to eight times longer than it is wide and a glint patch is round, so the filter accepts anything between 2:1 and 15:1 and at least 20 m across. Those bounds are deliberately loose, wide enough that an unusual hull still passes while a round blob cannot. And a satellite pass sometimes covers only part of an area, so any date where less than 95% of the area was photographed is dropped. Counting a partly covered area would look like a decline in traffic that never happened. In practice nothing landed in between: every date behind the figures came in at 100% coverage, and the rest were dropped outright.

A study fusing Sentinel-1 radar with AIS measured a 97% decline across the two weeks either side of the closure. We measured 95% from optical imagery over a longer window, with no AIS involved.

Oversized shapes and the size bound

About one detection in forty measures longer than 400 m, a length no ship afloat reaches. The longest spans 1,310 m. These shapes cluster by date: at Fujairah, 15 of the 21 fall on a single date, and in the corridor, 13 of 24 fall on two, which fits cloud or sun glint on those passes. At Kharg, a shape about 1,120 m long recurs at the same spot on five dates, which points to a fixed terminal structure missing from the land polygons. An upper bound alongside the existing area floor removes all of them:

WHERE length_m BETWEEN 50 AND 400

The correction tightens the result. The corridor returned nine detections on 18 March 2026, one of the two clustered dates, and five of those exceed 400 m. Four hulls survive the bound, so the corridor reads emptier after the closure.

Every vessel we found

The job writes its detections as GeoParquet and then as PMTiles, a single file a browser can read directly from cloud storage without a tile server. Each month becomes its own layer, so a before-and-after comparison means publishing the two months that matter.

Show the code

from wherobots import vtiles

feats = sedona.sql("""
    SELECT geometry, aoi, obs_date, length_m, area_m2,
           date_format(to_date(obs_date), 'yyyy-MM') AS layer
    FROM det
""")

vtiles.generate_pmtiles(feats, f"{base}/hormuz_monthly.pmtiles")

The three maps below show one area each, with the 2025 month in blue and the 2026 month in red. Each caption gives the month’s total and a per-pass rate, which divides that total by the number of usable passes.

The corridor. May 2025 in blue against May 2026 in red. 46 detections across four clear passes in 2025, then two across one pass: about 11 vessels a pass, down to two. Open fullscreen.

The anchorage outside. April 2025 in blue against April 2026 in red. 200 detections over three passes, then 122 over two, which is roughly 67 a pass and then 61. April appears here because no May 2026 pass covered the anchorage. Open fullscreen.

The loading terminal inside the Gulf, May 2025 in blue against May 2026 in red. Loading held: 22 detections over two passes in 2025 and 29 over three in 2026, or 11 a pass then about 10. Open fullscreen.

All three areas at once, with every month switched on: 3,117 detections in 20 layers, one per month from January 2025 to August 2026. Kharg sits to the north-west, the strait and Fujairah to the south-east. Step through the months in the sidebar to see the corridor empty while the anchorage and the terminal keep their traffic. Open fullscreen.

Compute cost and data read

719 Sentinel-2 scenes cross these three areas over 20 months, which is 151 GB of imagery. The analysis read about 1.5 GB of that, because RS_Clip fetches only the parts of each file that overlap the three areas. The scenes stay in the public archive throughout.

The 20-month run took 26.8 minutes on a MEDIUM runtime, roughly 25 Spatial Units at $1.50 each in us-west-2, or about $37. Widening the same job to the whole Gulf of Oman and lower Persian Gulf, 368,875 km² at full coverage, added $23 and seven minutes of wall clock. Calibration and tiling took it to about $63. Wherobots Cloud reports the exact figure for every run in Workload History.

The coordinates and date range are parameters, so pointing the job at another chokepoint is an edit to two arguments. The notebook has them in its configuration cell, alongside the size and shape limits. Sentinel-2 passes overhead every five days, and a scheduled run keeps the count current independent of AIS broadcasts.

The surveyed corridor is 25 km wide within a strait roughly 40 km across, so the 95% decline describes that corridor. Whether some vessels shifted to routes outside it cannot be settled from these three areas. A second limit sits in the detector. The 2,000 m² area floor runs before the seam merge, so a small hull lying on a tile boundary can split into two pieces that each fall under the floor and go uncounted. A hull touches a seam on roughly 2% of passes, the loss bites only near the size floor, and it applies to every date equally, so the year on year comparisons survive it. The downloadable notebook applies the area floor after the merge instead.

Get Started with Wherobots

Try Now

How to score every building in a state for catastrophe risk: an exploration project

Posted on August 31, 2026October 4, 2026 by Ben Pruden

Part 1 of a series: What’s possible?

I didn’t start this project to build a real risk model, not one that can be used by an insurer tomorrow. But it is a workflow that an insurer can put into practice to create their own risk scores with Wherobots and their own expertise. The idea was simple: what can be done in a day with Wherobots MCP for Claude and the VSCode plugin? It turns out quite a lot, and the barrier to entry is a Claude (or other coding agent) Account and VS Code running on a laptop with access to a Wherobots pro tier connected.

Total costs across tokens and data processing was about $300. This example looks at the Insurance industry and scoring risk on a property and building level. I’ll be following up with additional examples on Drone Delivery hub / nest site selection as well as communication services network connectivity in the coming weeks.

Lets start with the concept.

Build a Risk Explorer Data Pipeline with AWS & Wherobots

Start Now

What does an insurer need to understand when evaluating the risk of a location?

Insurance works on a basic idea: what is the risk at a given location, and how should it be priced? For a property and casualty carrier that question repeats across every building and property asset in its portfolio. The data to answer it, in many cases, already exists spread across wildfire rasters, hail radar, flood layers, parcels, and road networks. Getting it into one scored, queryable, mappable form, updated on a recurring schedule is usually the hard part, because handling these data types is unsupported and often impossible in other compute environments.

This post walks that process end to end: every Overture building in Colorado, 2,771,126 of them, joined against the Parcels they sit upon from Regrid, and scored across five perils, from a statewide view down to a single rooftop. The data processing engine used is WherobotsDB on Wherobots Cloud in AWS. Our team at Wherobots maintains Apache Sedona, and WherobotsDB is 100% code compatible the OSS, but fully managed with more spatial functions and capabilities, grater performance, and enterprise support options. Three capabilities carry this data pipeline: spatial SQL over raster and vector datasets combined, a spatial operation that prunes on reads time, and PMTiles that WherobotsDB can deliver to enable a snappy browser experience on MapLibre.

The finished map is live on the Wherobots Website here. You can play with it inline below:

The Colorado catastrophe-risk explorer, built with WherobotsDB, Claude Opus 5 & VSCode via agentic app development.

The data fed into the application layer

Seven sources feed the app and the data pipeline:

Layer	Source
Building footprints	Overture Maps
Wildfire hazard	USFS Wildfire Risk to Communities (30 m raster)
Hail	NOAA radar, 403M observations since 2016
Flood	FEMA National Flood Hazard Layer
Tornado / wind	FEMA National Risk Index (tract)
Parcels, land use	Regrid
Fire-station access	Overture places (road-network distance computed in the pipeline)

Two structural decisions:

A count is a floor. No measurement means no record, so a building with no hail observation is not a building at zero hail risk. Every count reads as “at least“.
The ranking: by building count, not by dollars. Assessor value coverage runs near 63% and is county dependent. El Paso, the largest county at 290,630 buildings, with $0 reported by the county. A dollar ranking would identify who reports property values, not where actual risk value is identified. A P&C insurer would have a clear understanding of building / parcel / asset value, so this diverts from what would be the likely standard they would want. These values could easily be added to this data pipeline, however.

The data pipeline: a medallion architecture for catastrophe risk scoring

The pipeline is a medallion architecture, and every layer is an Iceberg catalog table you can query on its own.

Bronze, raw as it lands, one table per source:

Table	What lands in it
`overture_buildings`	Overture footprints (2.77M)
`noaa_hail`	NOAA radar hail, 403M observations since 2016
`usfs_wildfire`	USFS Wildfire Risk to Communities, 30 m raster (COG)
`fema_nfhl`	FEMA National Flood Hazard Layer (flood zones)
`fema_nri`	FEMA National Risk Index, tract (tornado, wind)
`regrid_parcels`	Regrid parcels: reported parcel value, land use
`overture_places`	Overture fire-station locations

Silver, conformed: clipped to the Colorado boundary, coordinate systems reconciled, geometries validated: buildings, hail_grid, flood_zones, parcels, fire_stations, nri_tracts.

Gold, scored and tile-ready: building_perils (five peril scores per building), hex_rollup (H3 exposure and composite), parcels_scored.

Update frequency. Hail accumulation and Regrid parcels refresh monthly, so a monthly run keeps the map current. USFS wildfire and the FEMA National Risk Index refresh annually. The pipeline is defined with the Wherobots Python SDK and scheduled with Apache Airflow, one run per medallion stage.

Reading and clipping at scale with spatial predicate pushdown

The buildings layer is continental but the analysis is for Colorado. A spatial predicate filter pushdown on the read keeps makes it easy and cheap to run that state only:

-- Silver: every Colorado building footprint, clipped with a pushed-down spatial predicate
CREATE TABLE org_catalog.co_pc_risk_silver.buildings AS
SELECT b.id, b.geometry
FROM wherobots_open_data.overture_maps_foundation.buildings_building b
JOIN org_catalog.co_pc_risk_bronze.co_boundary aoi
  ON ST_Intersects(b.geometry, aoi.geometry);

WherobotsDB pushes ST_Intersects down to the GeoParquet files, reads each file’s bounding-box metadata, and skips the files that fall outside Colorado. Because WherobotsDB can easily index a dataset by proximity using say a Hilbert Curve index, this makes this type of operation even easier.

Fusing raster and vector in one query

This is a core WherobotsDB capability that this and many other use cases leverage. Wildfire hazard is a 30 m raster. Hail is millions of radar points. Flood is polygons. Buildings are polygons. One set-based query samples and joins all of them, per building:

-- Gold: score each building against the wildfire raster and the hazard vectors.
WITH wf AS (
  SELECT b.id, MAX(RS_ZonalStats(w.rast, b.geometry, 'max')) AS wildfire_score
  FROM org_catalog.co_pc_risk_silver.buildings     b
  JOIN org_catalog.co_pc_risk_silver.wildfire       w ON RS_Intersects(w.rast, b.geometry)
  GROUP BY b.id
),
hail AS (
  SELECT b.id, COUNT(*)  AS hail_events,
         MAX(h.size_in)  AS max_hail_in
  FROM org_catalog.co_pc_risk_silver.buildings     b
  LEFT JOIN org_catalog.co_pc_risk_silver.hail_obs h ON ST_Intersects(b.geometry, h.geometry)
  GROUP BY b.id
),
flood AS (
  SELECT b.id, ANY_VALUE(f.zone) AS flood_zone
  FROM org_catalog.co_pc_risk_silver.buildings       b
  LEFT JOIN org_catalog.co_pc_risk_silver.flood_zones f ON ST_Intersects(b.geometry, f.geometry)
  GROUP BY b.id
)
SELECT
  b.id,
  b.geometry,
  wf.wildfire_score,
  hail.hail_events,
  hail.max_hail_in,
  flood.flood_zone
FROM org_catalog.co_pc_risk_silver.buildings b
LEFT JOIN wf    ON b.id = wf.id
LEFT JOIN hail  ON b.id = hail.id
LEFT JOIN flood ON b.id = flood.id;

RS_ZonalStats reads the raster values under each footprint and returns key statistics. Each peril is computed in its own CTE so the joins stay one-to-many and the counts are exact. WherobotsDB executes RS_Intersects and ST_Intersects as spatial range joins, so the raster-to-building and point-to-building matches avoid a cross product. The run scores the full state in just a few minutes on one WherobotsDB medium runtime.

Five catastrophe risk scores per building

The fused query gives each building a raw hazard value per peril. Scoring turns those into five small integers on one row, so a building is one comparable object across every hazard. Each peril is scored on its own scale and the five sum to a 0-to-7 composite:

Peril	Range	What sets the top score
Wildfire	0-2	USFS hazard sampled under the footprint, banded high / moderate / low
Hail	0-2	Largest radar-detected hail size over the building: 2 at 2 inches or greater
Flood	0-1	Inside a FEMA Special Flood Hazard Area
Wind / tornado	0-1	FEMA National Risk Index rating for the tract
Access	0-1	Greater than five road miles of the nearest fire station

The perils are scored individually. A single blended number would reads like a hazard map while hail dominates the sum, because hail is Colorado’s most costly regular peril. Keeping five columns lets the map user and the copilot switch to one peril at a time and understand them specifically.

The assignment is set-based, one CASE per peril over the fused hazard values:

-- Gold: turn raw hazard values into five 0-based peril scores per building
SELECT
  id,
  geometry,
  -- wildfire: USFS hazard under the footprint, banded 0 / 1 / 2
  CASE WHEN wildfire_hazard >= hi_break THEN 2
       WHEN wildfire_hazard >= lo_break THEN 1 ELSE 0 END               AS wf_score,
  -- hail: largest radar-detected hail size over the building, in inches
  CASE WHEN max_hail_in >= 2.0 THEN 2
       WHEN max_hail_in >= 1.0 THEN 1 ELSE 0 END                        AS hail_score,
  -- flood: inside a FEMA Special Flood Hazard Area
  CASE WHEN in_sfha THEN 1 ELSE 0 END                                   AS flood_score,
  -- wind / tornado: FEMA National Risk Index rating for the tract
  CASE WHEN nri_wind_elevated THEN 1 ELSE 0 END                         AS wind_score,
  -- access: beyond five road miles (8047 m) of the nearest fire station
  CASE WHEN road_m IS NULL OR road_m > 8047 THEN 1 ELSE 0 END           AS access_score,
  -- composite: the five summed, domain 0-7
  ( CASE WHEN wildfire_hazard >= hi_break THEN 2
         WHEN wildfire_hazard >= lo_break THEN 1 ELSE 0 END
  + CASE WHEN max_hail_in >= 2.0 THEN 2
         WHEN max_hail_in >= 1.0 THEN 1 ELSE 0 END
  + CASE WHEN in_sfha THEN 1 ELSE 0 END
  + CASE WHEN nri_wind_elevated THEN 1 ELSE 0 END
  + CASE WHEN road_m IS NULL OR road_m > 8047 THEN 1 ELSE 0 END )       AS composite
FROM org_catalog.co_pc_risk_gold.building_hazards;

Two of these bands are exposed in a way worth highlighting.

Wildfire and hail are measured under each footprint; wind and flood are inherited. wf_score reads the USFS raster at the building, and hail_score reads the radar accumulation over it. wind_score comes from a FEMA National Risk Index rating at census-tract resolution, and flood_score from a FEMA zone polygon. A tract or a zone is far coarser than a rooftop, so those two scores are derived from an area where the detections exist, and that coarser grain is labeled in the schema, not intended to represent as a per-structure reading.

On access, a NULL distance means furthest, not missing. The access score flags every building beyond five road miles of a fire station. The nearest-station search runs to a 25 km radius; a building with no station inside that radius has a NULL distance and is the most remote in the state, so the rule scores NULL as 1. Writing road_m > 8047 alone drops exactly the worst-access cohort. That is why the score, not the raw distance, drives the map.

Run over all 2.77M buildings, the top wildfire band holds 307,679 buildings and the 2 in-and-up hail band holds 790,889. Both counts are floors: a building with no observation for a peril scores 0 for it, and reads as at least that exposed.

The accumulation view

The per-building table has fine grained details. A reinsurer often is concerned with concentration risk, and it’s easier to visualize when aggregated, so its helpful to roll the same scores up to an H3 grid:

-- Gold: roll per-building scores up to H3 for the accumulation view
SELECT
  ST_H3CellIDs(geometry, 7, false)[1] AS h3,
  COUNT(*)                            AS buildings,
  SUM(wildfire_score)                 AS wildfire_exposure,
  AVG(composite)                      AS mean_composite
FROM org_catalog.co_pc_risk_gold.building_perils
GROUP BY ST_H3CellIDs(geometry, 7, false)[1];

The same gold tables drive both the statewide hex view and the street-level building view. Nothing recomputes between zoom levels.

From a gold table to a map anyone can interact iwth

WherobotsDB builds the PMTiles as vector tiles and ships them to object storage, so a hand-written MapLibre GL JS page reads them by their S3 URL. No separate tiling tool, no tile server, no build step in the app.

# WherobotsDB builds and ships the PMTiles for the web app
from wherobots import WherobotsJob
WherobotsJob(script="s3://.../export_pmtiles.py", name="co-risk-export", runtime="small").submit()

The build-to-browser toolchain

The whole build was authored and shipped with one repeatable stack:

Build and version: Claude Code (Opus 5) inside VS Code, with the pipeline and the app committed to GitHub.
Compute and data: WherobotsDB on Wherobots Cloud runs the medallion pipeline; the Wherobots Python SDK defines each job and it is scheduled with Apache Airflow.
Tiles and hosting: WherobotsDB builds the PMTiles and ships them to Amazon S3; Vercel serves the static MapLibre GL JS page from its CDN.
Ask the map: the app includes a Claude copilot. Ask where is the worst wildfire exposure in El Paso County, and Claude switches the peril, recolors the map, and moves to show you the answer. It runs as a Vercel serverless function that calls the Anthropic API, with the analysis rules in its system prompt and the key held server-side.

What these catastrophe risk scores are, and what they are not

The five scores screen where to look. They are not a calibrated loss model, and nothing here is a determination of insurability. The bands, the floors, and the inherited-grain labels live in the pipeline and again in the copilot’s system prompt, so the map view and the copilot’s answers stay inside the same limits.

What the map surfaces

Hail is Colorado’s dominant peril, ahead of wildfire. In the top wildfire tier alone, 307,679 buildings screen as elevated: a starting list for where to focus. The access-to-services peril, powered by road-network distance to the nearest fire station, marks the insurability boundary, the cohort beyond a fire department’s reach. A P&C analyst can review all of it by toggling the map.

How to port this catastrophe risk scoring pipeline to any state or portfolio

The Colorado build is a template. Swap the boundary for your state or the whole country, swap the peril rasters and point clouds for the ones on your portfolio, and the pipeline holds: push the spatial predicate on the read, fuse raster and vector in one query, roll up to H3, let WherobotsDB build and ship the tiles down to the property parcel and building layer.

That is the point of an AI Context Engine for the Physical World: the physical-world data, prepared and scored at the scale a portfolio spans, in a form an analyst and an agent can both query.

Start on Wherobots Cloud, connect to the Wherobots Claude MCP, or open the Colorado risk explorer to see what I built, or check out the docs at docs.wherobots.com.

Get Started with Wherobots

Access Now

Key takeaways

The whole build took a day on a laptop with Claude and VS Code. The barrier to entry is a Claude account, VS Code, and a Wherobots Professional tier account. Part 1 of a series showing what one person can build in a day on this stack.
The project builds a scoring workflow an insurer would put into practice, not a calibrated risk model. It is a workflow an insurer runs to create their own risk scores with Wherobots and their own expertise. The five scores screen where to look. They do not determine insurability.
The layers already exist. Getting them into one scored, queryable form is the hard part. The data spreads across wildfire rasters, hail radar, flood layers, parcels, and road networks. WherobotsDB does the fusing in one query: spatial SQL over raster and vector, spatial predicate pushdown at read, and PMTiles the engine builds and ships to the browser.
The Colorado build is a template for any state or portfolio, and drone delivery and telecom examples will come next in the series. Swap the boundary for your state and the peril rasters for the ones on your portfolio.

Learning Wherobots by Building a National AI Data Center Suitability Report

Posted on August 27, 2026October 4, 2026 by George Chandeep Corea

This guest post is from a Wherobots user George Chandeep Corea, which covers an exploration of how he used Wherobots MCP with his preferred AI coding tools to build the interactive AI data center site suitability analysis tool below. The following quote is from his own biography.

Learning by Doing

Below is a working national suitability report for AI data centers. Move the weight sliders and the ranking reorders. Toggle the tailings-dam scenario and the buildable area changes. Everything in it was computed on Wherobots over 15.91 million authoritative geometries drawn from 16 national and state portals, then published as a zero-dependency interactive document that opens in any browser with no compute environment behind it.

Have a look first. The rest of this post explains how it was built and how I shaped it, as well as where I plan to go next.

Open the full report in a new tab →

I learn by doing. I built this on a private basis to contribute to the community discourse. For me it is all about the data; any opinions, perceived or intentional, are my own and have nothing to do with my official job roles. All data is from public sources.

It covers 17 candidate industrial sites across all 8 Australian states and territories, spanning the National Electricity Market and the SWIS grid.

Why I built this AI data center suitability model

For me, this is a community and conservation-minded project: not about saying no to development, but about helping identify where development can create the most value with the least harm. The point is to make siting decisions more transparent, more constructive, and easier to discuss across stakeholders.

AI data centers are resource-intensive. They depend on power, water, land, and supporting infrastructure, and they can also create real pressures on ecosystems and communities. I wanted to build something that makes those tradeoffs visible instead of hiding them behind a single score or reports that sit on shelves and don’t let users see it through there lenses.

Why I used Wherobots for spatial suitability at scale

I wanted a platform that could handle real spatial scale without forcing me into a slow, fragile workflow. I had been hearing about it on LinkedIn and Youtube by following Matt Forrest and wanted to learn by doing.

This project queried 15.91 million authoritative geometries across 16 national and state/territory spatial datasets. At the regional modeling scale, the underlying NSW stack ingested 1,751,315 geometries across multiple cloud spatial tables, evaluating 4.92 million spatial join combinations. Using cloud-native GeoParquet inside Wherobots Havasu (Spatially Aware Apache Iceberg) tables, the storage footprint compressed from ~2.9 GB raw equivalent to ~430.7 MB — an 85.2% reduction.

That scale matters because suitability analysis only becomes useful when it can be rerun quickly as assumptions change.

Verified runtime breakdown across 15.9M geometrics

Four different operations, four different numbers. They are not alternatives to each other: 2.4 s for the spatial SQL join execution (distributed joins + net developable area overlays across 1.75M+ geometries); 18.4 s → 3.2 s for the national scan across 15.91M geometries using Hilbert space-filling curve partitioning; 200.6 s for the cold end-to-end batch ETL (uncached ingest, GDA2020 reprojection, ST_MakeValid topology repair, Iceberg writes); and < 1 ms for the client-side What-If recalculation in the browser.

What I was trying to learn

I wanted to understand how far I could push a cloud spatial workflow while still producing something that a non-technical audience could review. The answer was to combine Wherobots for the heavy computation with a document-style HTML report for presentation.

The report is the public-facing layer of the analysis. It includes an interactive map, a ranked candidate leaderboard, benchmarking tables, provenance, methodology, and a what-if sandbox for changing ranking weights.

My recent work covers GDA2020 cadastral modernization at an enterprise level, environmental monitoring at a state level, and disaster modeling and high level coordination with FEMA. As the creator of *AuraSiting Crafter* (currently hunter_spatial crafter), I built an open-source multi-criteria engine benchmarking clean energy, water security, acoustic buffers, and developable land across 15.91 million Australian geometries using Wherobots MCP as the critical foundation.

George Chandeep Corea

Builder, Learner, Public Sector

Why the interactive suitability report matters

The HTML report is where the analysis becomes useful for public review. It is not just a visualization — it is a way to interrogate the assumptions behind the ranking.

The map is interactive, searchable, and configurable. Readers can change basemaps, add WMS/WFS services, and bring in external layers to compare candidate sites against their own data or public reference layers (some functionality coming soon!). That matters because spatial decisions are rarely made from one dataset alone. Adding data will make the report more useful for collaboration, because different stakeholders can overlay the information that matters to them: local context, community knowledge, infrastructure plans, environmental constraints, or jurisdictional boundaries.

The sandbox (already available) is especially valuable because readers can adjust the ranking weights and immediately see how the candidate ordering changes. That makes the tradeoffs tangible instead of theoretical.

Suitability report overview

The report overview combines summary KPI cards, an interactive What-If simulation sandbox, an interactive siting map with layer controls, and a ranked candidate leaderboard into a single reviewable spatial document.

How the suitability workflow works

The pipeline is straightforward:

Ingest authoritative layers — harvest from live government WFS/REST portals and national registers into cloud object storage
Apply exclusions and setbacks — deduct 30 m riparian corridors, 20 m high-pressure gas/water pipeline easements, high-value biodiversity zones, slope >5%, and tailings dam hazards for example.
Calculate network distances — topological winding distances (1.32× factor) to ≥132kV substations and wastewater outfalls for example
Evaluate sensitive receptors — continuous sigmoidal buffer penalty around schools, hospitals, and residential meshblocks for example
Benchmark against regional baselines — compare candidates against simulated baselines across interstate transition hubs (Latrobe Valley, Collie, Gladstone) for example. Only the hunter has detailed data. Sites outside the Development Application (DA) in LMCC has data simulated (not mocked or made up though) through analysing relevant real spatial data.
Publish a zero-dependency report — compile layers, audits, and metadata into a standalone interactive HTML document

The scoring model uses many weights:

Power grid proximity: 50%
Recycled water proximity: 30%
Parcel size: 20%
Others added/will be added as the project develops and

Those weights reflect the practical requirements of AI infrastructure while keeping environmental and resource considerations visible.

How the suitability score is calculated

Suitability = 0.40·S_power + 0.25·S_sensitive + 0.20·S_water + 0.15·S_size

Power grid proximity (40%). Proximity to ≥132kV transmission substations. Full points between 100 m and 500 m; linear decay to 0.0 at 5 km; a 0.7 penalty under 100 m for EMF safety and acoustic isolation.

Sensitive receptor buffer (25%). A continuous sigmoidal decay around a 500 m compliance threshold: S_sensitive(d) = 1 / (1 + e^(-0.01·(d - 500))). Hard exclusion under 300 m; acoustic-barrier buffer penalty 300–500 m (0.10–0.50); compliant 500 m–1.5 km (0.50–1.0); optimal workforce balance 1.5–5 km (1.0); linear commute decay beyond 5 km.

Recycled water proximity (20%). Distance to wastewater treatment plants for sustainable evaporative cooling, decaying linearly from 1.0 at ≤1 km to 0.0 at ≥10 km.

Developable parcel size (15%). Contiguous flat buildable land. Parcels ≥15 ha score 1.0; parcels under 3 ha get a 0.1 baseline, with linear interpolation between.

Slope calculations exclude land above 5% grade to avoid excessive earthworks.

What the analysis is based on

National input layers

The national model integrates 16 authoritative datasets across all 8 Australian jurisdictions: Geoscape National Cadastre & G-NAF (15,420,800), ABS 2021 Meshblocks & UCL (368,290), National Sensitive Receptors from ACARA / NHSD / OSM (47,510), Geoscience Australia electricity grid (4,820), BoM & GA surface water and outfalls (42,100), and state planning/cadastral portals — NSW SEED, QLD QSpatial, Vicmap, Landgate (27,725). Total: 15,911,245 geometries.

Those totals are the available universe published across the portals — 47,510 POIs and 15.4M cadastral parcels are what ACARA, NHSD, OpenStreetMap and Geoscape publish, not what this pipeline processed. The ingested and audited cohort is 1.75M regional features in the Hunter deep-dive, 368,290 ABS meshblocks partitioned nationally, 17 indexed candidate sites, and a ground-truth QA sample of 33 sensitive receptors (19 schools via ACARA, 14 hospitals via NHSD) sitting inside the candidate industrial zones across all 8 states, used to calibrate the sigmoidal buffer decay curve.

Of the 17 candidates, 4 are micro-sited in detail in the NSW Hunter; the remaining 13 are simulated baselines across interstate transition hubs (Latrobe Valley VIC, Collie WA, Gladstone QLD).

The earlier NSW regional stack: still the measured layer

Dataset	Feature count	Notes
NSW Transport Network (Rail)	275,421	Transport context
NSW Biodiversity Constraint Zones	262,258	Environmental exclusion
NSW Energy Grid Infrastructure	241,573	Power proximity
ABS Census Meshblocks	223,238	Demographic context
NSW Pipeline Corridors	197,247	Setback constraints
TfNSW Active Transport Pathways	188,576	Access and corridor context
NSW Hydrography & Waterways	181,501	Water and riparian context
ABS Regional Demographics	1,160	Regional benchmarking

Storage and compression on Havasu Iceberg

All spatial tables are cataloged under org_catalog.fgsdb.* on Wherobots Cloud and persisted in cloud object storage, in GDA2020 / MGA Zone 56 (EPSG:7856).

Table	Geometries	Uncompressed raw	GeoParquet footprint	Savings
macquarie_biodiversity_constraints	262,258	~580.0 MB	84.2 MB	85.5%
macquarie_energy_infrastructure	241,573	~420.0 MB	62.5 MB	85.1%
macquarie_transport_rail	275,421	~390.0 MB	58.1 MB	85.1%
macquarie_pipeline_corridors	197,247	~280.0 MB	41.8 MB	85.1%
macquarie_abs_meshblocks	223,238	~650.0 MB	98.4 MB	84.9%
macquarie_water_hydrography	181,501	~310.0 MB	44.6 MB	85.6%
macquarie_active_transport	188,576	~260.0 MB	39.2 MB	84.9%
abs_demographics	1,160	~12.0 MB	1.8 MB	85.0%
Total (8 tables)	1,751,315	~2.9 GB	~430.7 MB	85.2%

To avoid real-time WFS/FeatureServer REST API timeouts during Spark runs, every dataset is ingested and hosted as an optimized cloud spatial Iceberg table rather than fetched live.

The spatial SQL underneath

Building the net developable area mask — unioning riparian, pipeline, and rail buffers and subtracting them from the sub-precinct boundaries:

SELECT p.precinct_key,
       ST_Difference(p.geom, ST_Union_Aggr(c.geom)) AS net_developable_geom
FROM precinct_transform p
LEFT JOIN constraints c ON ST_Intersects(p.geom, c.geom)
GROUP BY p.precinct_key, p.geom

Computing nearest distances to transmission substations and wastewater treatment outfalls:

SELECT mb.mb_code21,
       MIN(ST_Distance(mb.mb_geom, ST_Transform(p.geometry, 'EPSG:4326', 'EPSG:7856'))) / 1000.0 AS dist_to_substation_km,
       MIN(ST_Distance(mb.mb_geom, ST_Transform(w.geometry, 'EPSG:4326', 'EPSG:7856'))) / 1000.0 AS dist_to_wwtw_km
FROM industrial_meshblocks mb
CROSS JOIN org_catalog.fgsdb.macquarie_energy_infrastructure p
CROSS JOIN org_catalog.fgsdb.macquarie_water_hydrography w
GROUP BY mb.mb_code21

How I think about the planning question

My background in conservation shaped this project. I’m not trying to use technology to say “no” to development. I’m trying to use it to ask better questions about where development belongs, where it creates value, and how to reduce unnecessary conflict between infrastructure, ecosystems, and communities.

That is why the report includes explicit constraint layers, benchmarking, and scenario controls. The aim is to support better judgment, not replace it.

Benchmarking and query speed mechanics

Benchmarking and speed mechanics breakdown detailing how spatial SQL queries execute across 1.75M+ geometries via metadata envelope pruning, Hilbert curve clustering, vectorized memory execution, and parallel distributed spatial joins.

Four things account for the query speed:

Metadata envelope pruning. Havasu Iceberg stores 2D bounding box envelopes directly inside Iceberg AVRO manifest files. Queries with spatial predicates such as ST_Intersects prune the large majority of irrelevant Parquet files at the metadata layer, before any raw disk bytes are scanned.

Hilbert curve spatial clustering. Geometries are sorted with 2D Hilbert space-filling curves during ingestion, so geographically adjacent features land in the same Parquet row groups and storage partitions. That removes random disk I/O seek overhead.

Vectorized memory execution. Apache Sedona operates directly on columnar GeoParquet WKB geometry buffers, avoiding serialization costs between Python, the Spark JVM, and native spatial drivers.

Parallel distributed spatial joins. Quad-tree and R-tree spatial indexes partition the query space across worker nodes, turning expensive O(N × M) cross-joins into O(N log M) parallel bucket joins.

Candidates are benchmarked against simulated local and regional baselines — Latrobe Valley in VIC, Collie in WA, Gladstone in QLD — to position NSW development opportunities within the wider national energy market transition.

How the what-if sandbox runs with no server

The simulation sandbox runs entirely in the browser: no server calls, no network latency, no cloud API charges during a slider session.

The heavy spatial work — topological winding distances, buffer overlaps, elevation head drops, thermodynamic decay rates — is computed at build time on Wherobots and embedded in the report’s JSON payload. Moving a weight slider then normalizes the raw values so the weights sum to 1.0 and recalculates suitability across every candidate record in JavaScript, in under a millisecond. Toggling the tailings dam safety switch swaps pad areas between the declared and de-declared cases (+15.2 ha unlocked) and updates map polygons and audit cards without refetching any GeoJSON. Candidates re-sort by active score, updating the leaderboard, marker radii, and popup badges.

The payoff is threefold: no per-query cloud compute cost during interactive sessions, full interactivity when the HTML report is emailed or opened offline, and slider manipulation that stays smooth without waiting on network round-trips.

What’s next for the suitability model

The following describes the author’s planned work, not shipped functionality.

Sensitive receptor scoring. The next version will score candidates against sensitive community receptors — schools and early childhood centers, hospitals and aged care, residential meshblocks, and workforce commute bands — to prevent acoustic, thermal, electromagnetic, and visual conflicts. The model uses a continuous sigmoidal penalty around a 500m critical setback threshold, with a workforce accessibility modifier that’s neutral in the 1.5km–5.0km commute band. In practice: under 300m is a critical exclusion, 300–500m carries a buffer penalty requiring an acoustic barrier, 500m–1.5km is compliant, 1.5–5.0km is the optimal community-and-workforce balance, and beyond 5km the score falls off for commute burden.

Grid and water policy analysis. Following the National Cabinet debate on AI data center power demand and regional grid security, the framework will evaluate proximity to 132kV, 330kV, and 500kV bulk transmission, flag candidates adjacent to retiring coal-fired stations that can reuse existing heavy transmission without expensive grid upgrades, map candidates against Renewable Energy Zones and firming assets for 24/7 clean energy matching, and restrict cooling supply strictly to recycled and wastewater sources — scoring zero for any site dependent on potable drinking water reserves or vulnerable aquifers. Candidates then sort into three tiers: grid-ready and fast-track eligible, conditional pending firming storage, or constrained by congestion and potable water reliance.

Open platform integration. To let the public explore these models interactively, hunter_spatial_crafter will integrate with opengeos/GeoLibre. The design point worth noting is zero duplication: Wherobots writes suitability layers to a central cloud bucket as GeoParquet and PMTiles, and GeoLibre’s in-browser DuckDB-WASM engine reads those exact same files over HTTP range requests — fetching only the row groups it needs, with no dataset conversion, no server-side copy, and no large client download.

Engineering efficiencies: reducing compute spend

Total cloud batch compute spend across dozens of iterative development, benchmark, and QA runs was ~$36 AUD (US$24.13) — roughly $1.03 per run across ~35 automated batch runs. Analyzing the execution profile shows how a production pipeline could cut that further:

Decouple heavy geometry joins from lightweight scoring. The pipeline splits into a compute-intensive geometric tier (ingest, GDA2020 reprojection, ST_MakeValid repair, 30 m/20 m buffers, ST_Difference overlays across millions of polygons) and a compute-light scoring tier (decay curves and weighted composites over precomputed distance attributes). Structured as a DAG with intermediate materialised GeoParquet stages, tuning a weight or the sigmoidal threshold never re-runs the heavy joins — only the downstream matrix recalculates.

Fingerprint sources and memoize snapshots. Baseline layers change infrequently. Content hashing (ETags, GeoParquet file hashes, Iceberg snapshot manifest IDs) lets untouched tables be skipped, reading straight from cached Havasu Iceberg partitions.

Process delta partitions. When state portals publish quarterly cadastral updates, Iceberg’s ACID snapshot metadata lets WherobotsDB (optimized and managed Apache Sedona) isolate only modified parcel geometries instead of full continental scans.

Offload interactive compute to the client. Compiling precomputed distance topologies into the standalone report puts millions of interactive public scenario evaluations at $0.00 cloud compute cost.

Applied together, these reduce continuous CI/CD pipeline cost from ~$36 AUD to under $5 AUD.

Why I think this is useful

This project taught me that Wherobots is more than a faster way to run spatial SQL. It is a way to make large-scale spatial analysis practical for real decision support. It gave me a path from raw geodata to a polished, reviewable spatial document — one that can be opened by a stakeholder, examined by a technical reviewer, and discussed in public.

For a personal project, that was the goal: learn the platform, test it at real scale, create something useful, and contribute to the broader discussion about responsible AI infrastructure.

If you’re interested in spatial analysis, cloud geospatial workflows, or responsible AI infrastructure, feel free to connect with me on LinkedIn.

Key takeaways

Wherobots computed a national suitability model over 15.91 million authoritative geometries drawn from 16 national and state portals, published as a zero-dependency interactive document that opens in any browser with no compute environment behind it.
Cloud-native GeoParquet inside Wherobots Havasu (Spatially Aware Apache Iceberg) tables compressed the storage footprint from ~2.9 GB raw equivalent to ~430.7 MB, an 85.2% reduction.
Runtime at national scale: 2.4 s for the spatial SQL join execution, 18.4 s down to 3.2 s for the national scan across 15.91M geometries using Hilbert space-filling curve partitioning, 200.6 s for the cold end-to-end batch ETL, and under 1 ms for the client-side What-If recalculation in the browser.
Total cloud batch compute spend across dozens of iterative development, benchmark, and QA runs was ~$36 AUD (US$24.13), roughly $1.03 per run across ~35 automated batch runs.
The model covers 17 candidate industrial sites across all 8 Australian states and territories, scored on the formula Suitability = 0.40·S_power + 0.25·S_sensitive + 0.20·S_water + 0.15·S_size.

Get Started with Wherobots

Try Now

From the Spokane firestorm to all of Washington: real-time wildfire monitoring for under $50 a pass

Posted on August 11, 2026October 4, 2026 by Jia Yu

On August 1, 2026, the most destructive fire event in Washington’s history swept into the Spokane area, driving 65,000 people from their homes. Wildfire response runs on one thing: knowing where the fire has burned and where it is growing. Every hour sooner that map exists, the better every decision downstream, from evacuation zones and containment lines to protecting homes, utilities, and insured assets.

The satellites already see every fire. What has been missing is real-time analysis: the moment a new pass lands in the archive, this pipeline turns it into a burn-severity map in 14 minutes, for about $20 of compute. One Python job on WherobotsDB, no downloads, no ingest, no GIS backlog. We built it on the Spokane firestorm, then pointed the same job at the entire state, and it found Washington’s other major fires (Sinlahekin, Modrite, Kaiser Canyon) by itself, while they were still burning. The pipeline keeps up with whatever imagery you feed it: the free Sentinel-2 constellation delivers a fresh statewide look every few days, and higher-cadence sources slot into the same job. Every number below comes from the production runs.

Burn severity of the Spokane firestorm: Sentinel-2 dNBR, August 6 2026

The technique is dNBR (differenced Normalized Burn Ratio), the standard used by USGS and the EU’s Copernicus Emergency Management Service for burn severity mapping. Healthy vegetation is bright in near-infrared (Sentinel-2 band B8A) and dark in shortwave-infrared (B12); burned ground flips that. So:

NBR = (B8A − B12) / (B8A + B12), per pixel, per scene
dNBR = NBR_pre-fire − NBR_post-fire; big positive values mean severe burn

The input is 56 cloud-optimized GeoTIFFs (~2.6 GB) in the public Sentinel-2 archive on S3, and none of it is downloaded or ingested. WherobotsDB’s out-of-database rasters resolve lazily, so only the ~9% of pixels inside our area of interest ever cross the wire.

Which scenes detected the fire?

WherobotsDB’s STAC reader loads a STAC collection as a DataFrame, with spatial and temporal filter pushdown. We point it at Element84’s Earth Search catalog and filter to a bounding box covering all three fires, a pre-fire window (July 10–31) and a post-fire window (August 2 onward). Each STAC asset comes with a ready-to-use out-db raster column.

Show the code

AOI = "POLYGON((-117.75 47.55, -117.20 47.55, -117.20 47.95, -117.75 47.95, -117.75 47.55))"
AOI_SQL = f"ST_SetSRID(ST_GeomFromWKT('{AOI}'), 4326)"

raw = (
    sedona.read.format("stac")
    .option("itemsLimitMax", "1000")
    .load("https://earth-search.aws.element84.com/v1/collections/sentinel-2-c1-l2a")
)
raw.createOrReplaceTempView("stac_items")

scenes = sedona.sql(f"""
    SELECT id,
           element_at(split(id, '_'), 2) AS tile,
           datetime,
           CASE WHEN datetime <= TIMESTAMP '2026-07-31 23:59:59'
                THEN 'pre' ELSE 'post' END AS phase,
           `eo:cloud_cover` AS cloud_cover,
           assets['nir08'].rast AS nir,    -- B8A, 20 m
           assets['swir22'].rast AS swir   -- B12, 20 m
    FROM stac_items
    WHERE ST_Intersects(geometry, {AOI_SQL})
      AND datetime >= TIMESTAMP '2026-07-10 00:00:00'
      AND datetime <= TIMESTAMP '2026-08-07 23:59:59'
      AND (datetime <= TIMESTAMP '2026-07-31 23:59:59'
           OR datetime >= TIMESTAMP '2026-08-02 00:00:00')
""")
scenes.createOrReplaceTempView("scenes")
scenes.select("id", "tile", "datetime", "phase", "cloud_cover").orderBy("tile", "datetime").show(5)

+------------------------------+------+-----------------------+-----+-----------+
|id                            |tile  |datetime               |phase|cloud_cover|
+------------------------------+------+-----------------------+-----+-----------+
|S2C_T11TMN_20260712T185514_L2A|T11TMN|2026-07-12 19:01:16.775|pre  |99.997228  |
|S2B_T11TMN_20260714T184425_L2A|T11TMN|2026-07-14 18:51:18.916|pre  |88.490695  |
|S2A_T11TMN_20260714T190034_L2A|T11TMN|2026-07-14 19:01:32.559|pre  |89.065939  |
|S2B_T11TMN_20260717T185644_L2A|T11TMN|2026-07-17 19:01:16.249|pre  |12.527606  |
|S2C_T11TMN_20260719T184029_L2A|T11TMN|2026-07-19 18:51:20.68 |pre  |3.3E-4     |
+------------------------------+------+-----------------------+-----+-----------+
only showing top 5 rows

28 scenes match, across two MGRS tiles (T11TMN to the south, T11UMP to the north). With Sentinel-2A, -2B, and -2C all flying, Spokane is revisited every two to three days. That cadence is what lets the pipeline track a fire’s status over time.

How burned is a pixel?

Each scene’s B8A and B12 arrive as separate single-band rasters. We clip both to the AOI with RS_Clip, which keeps every downstream operation working on a ~4 M-pixel window instead of a full 30 M-pixel tile. The band math itself is a vectorized Python UDF: the function takes one SedonaRaster argument per raster column, reads pixels as NumPy arrays with nodata masked to NaN, and returns a new raster with with_bands, which keeps the CRS and transform intact (a follow-up RS_SetBandNoDataValue marks the sentinel for the statistics downstream). WherobotsDB feeds the UDF through Apache Arrow in batches, so plain NumPy runs at cluster scale, and everything in the SciPy ecosystem is available inside the function when an analysis needs it.

The raster reader already applies the COG’s raster:bands scale and offset, so pixel values arrive as surface reflectance (≈0–1). Use them directly; no DN conversion is needed. If an index map ever looks suspiciously flat, a one-second RS_Value probe at a known point tells you why.

Show the code

import numpy as np
from pyspark.sql.functions import col, expr
from sedona.spark.raster import SedonaRaster
from sedona.spark.sql.functions import sedona_vectorized_udf
from sedona.spark.sql.types import RasterType

NODATA = -9999.0

@sedona_vectorized_udf(return_type=RasterType())
def nbr_udf(nir: SedonaRaster, swir: SedonaRaster) -> SedonaRaster:
    """NBR = (B8A - B12) / (B8A + B12).

    as_numpy_masked() turns nodata into NaN, so invalid pixels ride through
    the arithmetic; values <= 0 (clip fill, dark pixels) join them.
    """
    n = nir.as_numpy_masked()[0]
    s = swir.as_numpy_masked()[0]
    with np.errstate(divide="ignore", invalid="ignore"):
        nbr = np.where((n <= 0) | (s <= 0), np.nan, (n - s) / (n + s))
    return nir.with_bands(np.where(np.isnan(nbr), NODATA, nbr)[np.newaxis])

clipped = sedona.sql(f"""
    SELECT id, tile, datetime, phase, cloud_cover,
           RS_Clip(nir, 1, {AOI_SQL}, false, 0.0) AS nir_aoi,
           RS_Clip(swir, 1, {AOI_SQL}, false, 0.0) AS swir_aoi
    FROM scenes
""").filter("nir_aoi IS NOT NULL AND swir_aoi IS NOT NULL")

nbr = clipped.select(
    "id", "tile", "datetime", "phase", "cloud_cover",
    nbr_udf(col("nir_aoi"), col("swir_aoi")).alias("nbr_raw"),
).withColumn(
    "nbr", expr(f"RS_SetBandNoDataValue(nbr_raw, {NODATA}d)")
).drop("nbr_raw")
nbr.persist()
nbr.createOrReplaceTempView("nbr")

The value <= 0 guard handles nodata, clip fill, and the offset-negative dark pixels (water, deep shadow) in one condition. From here, one RS_ZonalStatsAll per scene gives a mean-NBR time series over the AOI: a quick health check on the imagery now, and a recovery tracker later.

Show the code

timeseries = sedona.sql(f"""
    SELECT id, tile, date(datetime) AS date, phase, cloud_cover,
           zs.mean AS mean_nbr, zs.count AS valid_px
    FROM (SELECT *, RS_ZonalStatsAll(nbr, {AOI_SQL}) AS zs FROM nbr)
    ORDER BY tile, date
""").localCheckpoint()   # materialize the 28-row table and cut its lineage
timeseries.createOrReplaceTempView("nbr_ts")

Mean NBR per clear scene, July 17 to August 6, with ignition marked. — **Mean NBR per clear scene, July 17 to August 6**, with ignition marked. Eight clear scenes in three weeks: that revisit rate is what turns a one-off map into a monitor.

Which two scenes? (the boring part that matters)

dNBR needs one pre-fire and one post-fire scene per tile, and the naive choice (“latest scene on each side”) bit us: Sentinel-2 sometimes acquires two overlapping products on the same day, and the latest pre-fire scene for T11TMN was a partial swath covering only a third of the tile. The result was a nodata band right through the Old Trails burn scar.

The selection rule: (1) drop scenes with 20% cloud or more, (2) gate on near-full swath coverage using the zonal pixel counts we already computed, and (3) only then prefer recency. Two window functions express the whole rule in SQL: max(valid_px) OVER each tile and phase sets the coverage bar, and row_number() takes the best-ranked scene on each side. The localCheckpoint() above matters here; it hands the planner a small, static table for this self-join while the raster work stays lazy.

Show the code

pairs = sedona.sql("""
    WITH scored AS (
        SELECT tile, phase, id, date, valid_px,
               max(valid_px) OVER (PARTITION BY tile, phase) AS max_px
        FROM nbr_ts
        WHERE cloud_cover IS NULL OR cloud_cover < 20
    ),
    ranked AS (
        SELECT *, row_number() OVER (
            PARTITION BY tile, phase
            ORDER BY (valid_px >= 0.95 * max_px) DESC,
                     date DESC, valid_px DESC, id DESC) AS rn
        FROM scored
    )
    SELECT pre.tile,
           pre.id AS pre_id, pre.date AS pre_date,
           post.id AS post_id, post.date AS post_date
    FROM ranked pre
    JOIN ranked post
      ON pre.tile = post.tile
     AND pre.phase = 'pre' AND post.phase = 'post'
     AND pre.rn = 1 AND post.rn = 1
""")
pairs.createOrReplaceTempView("pairs")

if pairs.count() == 0:
    raise ValueError("no usable scene pairs: every candidate was cloud- or coverage-filtered")

@sedona_vectorized_udf(return_type=RasterType())
def dnbr_udf(pre: SedonaRaster, post: SedonaRaster) -> SedonaRaster:
    """dNBR = pre-fire NBR - post-fire NBR; NaN stays NaN through the subtraction."""
    diff = pre.as_numpy_masked()[0] - post.as_numpy_masked()[0]
    return pre.with_bands(np.where(np.isnan(diff), NODATA, diff)[np.newaxis])

dnbr = sedona.sql("""
    SELECT p.tile, p.pre_date, p.post_date,
           pre.nbr AS pre_nbr, post.nbr AS post_nbr
    FROM pairs p
    JOIN nbr pre  ON pre.id  = p.pre_id
    JOIN nbr post ON post.id = p.post_id
""").select("tile", "pre_date", "post_date",
            dnbr_udf(col("pre_nbr"), col("post_nbr")).alias("dnbr_raw"))
dnbr = dnbr.withColumn("dnbr", expr(f"RS_SetBandNoDataValue(dnbr_raw, {NODATA}d)")).drop("dnbr_raw")
dnbr.persist()
dnbr.createOrReplaceTempView("dnbr")

How bad, and over how many acres?

USGS thresholds turn dNBR into severity classes (low ≥ 0.10, moderate-low ≥ 0.27, moderate-high ≥ 0.44, high ≥ 0.66). One more UDF pass classifies the raster, and RS_CountValue turns class pixel counts into acres. At 20 m resolution each pixel is 400 m², about a tenth of an acre:

Show the code

@sedona_vectorized_udf(return_type=RasterType())
def classify_udf(dnbr: SedonaRaster) -> SedonaRaster:
    """USGS dNBR classes: 1 low, 2 moderate-low, 3 moderate-high, 4 high.

    NaN compares False against every threshold, so nodata lands in class 0.
    """
    v = dnbr.as_numpy_masked()[0]
    with np.errstate(invalid="ignore"):
        out = np.select(
            [v >= 0.66, v >= 0.44, v >= 0.27, v >= 0.10],
            [4, 3, 2, 1],
            default=0)
    return dnbr.with_bands(out.astype(np.float64)[np.newaxis])

classified = dnbr.withColumn("classified", classify_udf(col("dnbr")))
classified.createOrReplaceTempView("classified")

ACRES = 400.0 / 4046.8564224

severity = sedona.sql(f"""
    SELECT tile, pre_date, post_date,
           ROUND(RS_CountValue(arr, 1) * {ACRES}, 1) AS low_acres,
           ROUND(RS_CountValue(arr, 2) * {ACRES}, 1) AS mod_low_acres,
           ROUND(RS_CountValue(arr, 3) * {ACRES}, 1) AS mod_high_acres,
           ROUND(RS_CountValue(arr, 4) * {ACRES}, 1) AS high_acres
    FROM (
        SELECT tile, pre_date, post_date,
               RS_BandAsArray(classified, 1) AS arr
        FROM classified
    )
""")
severity.show(truncate=False)

+------+----------+----------+---------+-------------+--------------+----------+
|tile  |pre_date  |post_date |low_acres|mod_low_acres|mod_high_acres|high_acres|
+------+----------+----------+---------+-------------+--------------+----------+
|T11TMN|2026-07-29|2026-08-06|6017.8   |3335.0       |1788.9        |502.9     |
|T11UMP|2026-07-29|2026-08-06|4372.5   |1663.2       |677.6         |193.8     |
+------+----------+----------+---------+-------------+--------------+----------+

Five days after ignition, the burn reads like this (the two tiles overlap in the 47.7–47.9°N band where the fires sit, so the rows must not be summed):

Tile	Pair	Low	Mod-low	Mod-high	High	Moderate+
T11TMN (south: NW Spokane, Nine Mile, Airway Heights)	Jul 29 → Aug 6	6,018	3,335	1,789	503	5,627
T11UMP (north: Suncrest, Mead)	Jul 29 → Aug 6	4,373	1,663	678	194	2,535

The map at the top of this post is the dNBR raster classified with these thresholds. The job also writes it out as GeoTIFFs with RS_AsGeoTiff and the raster writer, so the same output drops straight into QGIS. The main scar straddles the Spokane River around Nine Mile Falls and the Indian Trail bluffs, with the high-severity core on the slopes above the river. The low class (0.10–0.27) is the noisiest; at these thresholds it also picks up harvest and irrigation changes in cropland, so the moderate-and-worse figure is the one to quote. Official incident reports put the combined fires at ~10,000 acres of perimeter, which is expected to exceed spectrally-burned pixel area.

Burned acres by severity class, five days after ignition, per Sentinel-2 tile. — **Burned acres by severity class, five days after ignition**, per Sentinel-2 tile.

The same statements, one hundred times the area

Nothing in the pipeline knows it is about Spokane; the AOI is a WKT string. We swapped the 30 km box for all of Washington plus the Idaho panhandle (42 Sentinel-2 tiles, 690 candidate scenes, three UTM zones) and let the engine do what engines do:

**SCL-masked dNBR wall-clock: region scale is a parameter.** `RS_TileExplode` splits each scene into 512-pixel sub-tiles, so thousands of small tasks fill the cluster. That pattern is what keeps the statewide run this fast.

**457 detections in one run: the big ones are all named fires.** Dot area scales with detected acres (Jul 27 to Aug 6 change). Gray dashed dots are cropland-flagged by ESA WorldCover, 403 of 457; the 54 survivors are the fire candidates, and every news-confirmed fire is among them. The dashed box is the query AOI.

Explore every detection, live

The map below is the real output of the statewide run. Our production monitoring job, the same pipeline described above extended with burn-cluster detection and the cropland flag, ends by writing its detections as a single PMTiles archive with one call, vtiles.generate_pmtiles(detections, s3_path). A browser reads that file directly from S3 with HTTP range requests, no tile server, no export pipeline, rendered here in the Wherobots visualization tool. Click any polygon for its properties. If the viewer below is not working, try fullscreen mode here.

Did it work? Check the news

A detector is only as good as its false-positive rate, so we cross-checked the largest region-wide detections against news coverage of Washington’s August 2026 fire outbreak, the one that triggered the state’s first-ever “particularly dangerous situation” fire-weather warning.

Detection (Aug 6 imagery)	Identity	Cropland flag
69,022 ac near Tonasket	Sinlahekin Fire: 141,411 ac by Aug 10; ours is the Jul 27 → Aug 6 growth slice	clear: fire candidate
40,485 ac near Inchelium	Modrite Fire: “over 40,000 acres” in the same window	clear: fire candidate
Two clusters near Nespelem	Kaiser Canyon Fire flanks: ignited Jul 16, before our baseline, so only fresh growth appears	clear: fire candidate
4,512 ac northwest of Spokane	The Spokane complex this post began with, rediscovered without being told it exists	clear: fire candidate
46,947 ac in the Palouse	No fire in the news; wheat harvest between the two dates	flagged: cropland
11,073 ac Skagit delta · 6,219 ac Columbia Basin	No fires in the news; farmland and irrigation change	flagged: cropland

Every news-confirmed fire lands in the clear column; every no-news detection lands in the flagged column. At region scale the flag does serious work: 403 of 457 raw detections are cropland-majority, leaving 54 fire candidates from one extra join against ESA WorldCover, in the same engine that did everything else.

What didn’t happen here

Most remote-sensing pipelines spend their time and money on steps this analysis skipped:

No ingest. The STAC catalog was queried live and the COGs stayed where ESA/Element84 publish them. There is no “download the scenes” step, no staging bucket of copies to manage, and adding the next satellite pass to the analysis costs nothing but a date-range change.
No raster database to load. The out-db raster model means RS_Clip, the NBR UDF, and RS_ZonalStatsAll executed against windowed range-reads of the source files: ~2.6 GB of input resolved to only the ~9% of pixels inside the AOI.
Open formats on both ends. Results landed back in the lakehouse as GeoParquet and GeoTIFF. The severity table is queryable by any engine that reads Parquet, and the dNBR map opens directly in QGIS. Nothing is locked in a proprietary store.
One engine for rasters and vectors. The same SQL session can join these burn pixels against Iceberg/Havasu tables, like Overture Maps buildings in the Wherobots Open Data catalog, which is where this analysis goes next.

That combination (query-in-place over open cloud archives, spatial SQL over both raster and vector, and open outputs) is the lakehouse difference: the time from “a fire happened” to “a severity map exists” is measured in minutes of compute, not days of data wrangling.

From one-off analysis to standing watch

The analysis lives in one Python job script submitted with the Wherobots Airflow provider. New Sentinel-2 scenes land every two to three days, so a scheduled DAG keeps the severity map and the recovery time series current with zero manual steps:

run = WherobotsRunOperator(
    task_id="run_dnbr_job",
    name="spokane_fire_dnbr_{{ ds_nodash }}",
    region=Region.AWS_US_WEST_2,   # same region as the Sentinel-2 COGs
    runtime=Runtime.MEDIUM,
    run_python={"uri": f"{SCRIPT_S3}", "args": ["--post-end", "{{ ds }}"]},
    poll_logs=True,
)

End to end (STAC discovery, 28 NBR rasters, pair selection, dNBR, classification, acreage, Parquet + GeoTIFF outputs) takes about 14 minutes on a MEDIUM runtime and costs about $20: roughly 13 Spatial Units at $1.50 each in us-west-2. Wherobots Cloud reports the exact cost and Spatial Unit consumption of every run in Workload History, so the price of a standing watch is a line item you can read, per run, per day. The statewide watch is the same arithmetic: the Washington plus North Idaho detection run behind the live map above came to about $45 (30 Spatial Units on a LARGE runtime, 14 minutes for 42 tiles across three UTM zones).

Where to take it

Mask smoke and cloud per-pixel with the SCL band instead of scene-level metadata.
Vectorize the burn perimeter with RS_Polygonize on the classified raster.
Count affected structures by joining the perimeter against Overture Maps buildings, already hosted in the Wherobots Open Data catalog.
Track recovery: the same NBR time series that found the drop will show vegetation green-up over the coming seasons.

Everything here runs on any Wherobots organization. The imagery is public, and the code above is the entire method. If you want to try it on a fire (or flood, or storm) near you, create a Wherobots Cloud organization and point the STAC reader at your own AOI.

Key takeaways

On August 1, 2026, the Spokane firestorm drove 65,000 people from their homes. The pipeline turns a new Sentinel-2 pass into a burn-severity map in about 14 minutes for about $20 of compute: one Python job on WherobotsDB, with no downloads, no ingest, and no GIS backlog.
Input is 56 cloud-optimized GeoTIFFs (~2.6 GB) in the public Sentinel-2 archive on S3. WherobotsDB out-of-database rasters resolve lazily, so only about 9% of pixels inside the area of interest ever cross the wire. End to end on a MEDIUM runtime costs about $20 (roughly 13 Spatial Units at $1.50 each in us-west-2).
The same job pointed at Washington plus the Idaho panhandle (42 Sentinel-2 tiles, 690 candidate scenes, three UTM zones) finished in 14 minutes on a LARGE runtime for about $45 (30 Spatial Units). It produced 457 detections; 403 were cropland-majority via ESA WorldCover, leaving 54 fire candidates that included every news-confirmed fire in the window.
dNBR uses USGS thresholds (low ≥ 0.10, moderate-low ≥ 0.27, moderate-high ≥ 0.44, high ≥ 0.66). Five days after ignition, tile T11TMN showed 5,627 acres at moderate-and-worse; the low class also picks up harvest and irrigation change, so the post quotes moderate-and-worse rather than summing all classes.