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

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.

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.

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.

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

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.

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:

@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):

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.

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:

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.

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.