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.

Change Detection Using AlphaEarth Foundations (Part 2)

Posted on May 6, 2026October 4, 2026 by Ben Pruden

Author: Len Strnad

AlphaEarth Foundations embeddings gave us an initial view of temporal change in our earlier post on agricultural dynamics, but that first pass also exposed a limitation: raw embedding distances do not live on a consistent scale across pixels or land-cover types. As a result, it is hard to compare how unusual a given year looks from place to place, even when the relative signal is informative.

In this follow-up, we ask a more targeted question: can we score each annual embedding by how much it stands out from the rest of its local time series? To do that, we compute leave-one-out (LOO) medoid scores in embedding space and apply robust z-score scaling through time for each pixel. This gives us a more stable and interpretable way to identify unusual years across a short annual record.

As in the earlier AlphaEarth Foundations agriculture post, the goal is not to present a definitive change-detection benchmark. Instead, we show a practical workflow for surfacing unusual temporal states in AEF embeddings. We run that workflow against the global AEF Zarr mosaic hosted on Source Cooperative, build a multiscale output for exploration, and inspect several examples around Boulder, Colorado.

Global Zarr AEF Mosaic

For this follow-up we use the global AlphaEarth Foundations Zarr mosaic hosted on Source Cooperative. Thanks to Taylor Geospatial and Source Cooperative for making the dataset available in a form that is easy to access lazily in a single python function call!

The mosaic stores annual AEF embeddings at 10-meter resolution in EPSG:4326, with nine temporal snapshots from 2017 through 2025. Working from one global store keeps the workflow simple: we can subset the region of interest, compute scores, and publish a derived multiscale product without assembling separate local inputs.

import xarray as xr

ds = xr.open_zarr(
    "s3://us-west-2.opendata.source.coop/tge-labs/aef-mosaic/",
    consolidated=False,
    storage_options={"anon": True},
)
ds

/Users/len/code/wherobots-notebook-blogs/.pixi/envs/default/lib/python3.14/site-packages/zarr/core/group.py:3559: ZarrUserWarning: Object at README.md is not recognized as a component of a Zarr hierarchy.
  warnings.warn(
/Users/len/code/wherobots-notebook-blogs/.pixi/envs/default/lib/python3.14/site-packages/zarr/core/group.py:3559: ZarrUserWarning: Object at .checkpoint is not recognized as a component of a Zarr hierarchy.
  warnings.warn(

<xarray.Dataset> Size: 4PB
Dimensions:     (time: 9, band: 64, y: 1859584, x: 4009984)
Coordinates:
  * time        (time) int32 36B 2017 2018 2019 2020 2021 2022 2023 2024 2025
  * band        (band) object 512B 'A00' 'A01' 'A02' 'A03' ... 'A61' 'A62' 'A63'
  * y           (y) float64 15MB 83.69 83.69 83.69 ... -83.36 -83.36 -83.36
  * x           (x) float64 32MB -180.0 -180.0 -180.0 ... 180.2 180.2 180.2
Data variables:
    embeddings  (time, band, y, x) int8 4PB dask.array<chunksize=(1, 64, 256, 256), meta=np.ndarray>
Attributes: (12/15)
    proj:code:               EPSG:4326
    spatial:dimensions:      ['y', 'x']
    spatial:transform:       [8.983111749910169e-05, 0.0, -180.0, 0.0, -8.983...
    spatial:transform_type:  affine
    spatial:bbox:            [-180.0, -83.36280346631479, 180.2213438735178, ...
    spatial:shape:           [1859584, 4009984]
    ...                      ...
    geoemb:model:            https://developers.google.com/earth-engine/datas...
    geoemb:source_data:      https://source.coop/tge-labs/aef/v1/annual/
    geoemb:data_type:        int8
    geoemb:gsd:              8.983111749910169e-05
    geoemb:quantization:     {'method': 'signed_square', 'original_dtype': 'f...
    zarr_conventions:        [{'uuid': 'f17cb550-5864-4468-aeb7-f3180cfb622f'...

At first glance, 3.81 PiB and 1,024,049,664 chunks in S3 sounds impractically large. For this demo, though, a few details keep the workflow manageable:

only chunks that intersect land are present
the dataset is sharded into nine spatial chunks, which reduces the effective object count
we subset Colorado rather than operating on the full global mosaic

RasterFlow Tools Used

RasterFlow is Wherobots’ serverless inference engine for Earth Observation data. It builds inference-ready mosaics from remote sensing data and runs distributed workflows at scale. This notebook uses two RasterFlow capabilities:

run_mosaics_change on the RasterFlow Remote client for computing difference and LOO medoid scores over time
build_zarr_multiscales for producing a visualization-ready Zarr mosaic

Change Detection with RasterFlow Remote

For this run, we call run_mosaics_change on the RasterFlow Remote client. The method takes an input Zarr mosaic and writes out a derived Zarr mosaic containing one or more temporal change scores.

Here we compute both difference-over-time and LOO medoid scores in AEF embedding space, with robust z-score normalization applied across time for each pixel. The snippet below matches the current remote-client API used for this Colorado example.

from affine import Affine
from rasterio.transform import rowcol

# construct a bounding box selector for Colorado
transform = Affine(*ds.attrs["spatial:transform"])
row_min, col_min = rowcol(transform, -109.060253, 41.003444)
row_max, col_max = rowcol(transform, -102.041524, 37.000232)
selector = {"x": [col_min, col_max], "y": [row_min, row_max]}
selector

{'x': [np.int32(789701), np.int32(867833)],
 'y': [np.int32(475138), np.int32(519702)]}

from rasterflow_remote import DistanceMetricEnum, RasterflowClient, TemporalScoreMethodEnum

client = RasterflowClient()

result = client.run_mosaics_change(
    store="s3://us-west-2.opendata.source.coop/tge-labs/aef-mosaic/",
    score_methods=[
        TemporalScoreMethodEnum.DIFFERENCE,
        TemporalScoreMethodEnum.LOO_MEDOID,
    ],
    distance_metric=DistanceMetricEnum.ALPHA_EARTH_COSINE,
    data_var="embeddings",
    robust_time_zscore=True,
    selector=selector,
)

LOO Medoid Scores for Small Samples

To measure how unusual a given annual embedding is relative to its temporal context, we use a leave-one-out medoid (LOO medoid) score in embedding space. This is well suited to the AEF setting here, where we only have nine annual observations and want a reference that is less sensitive to outliers than a simple average.

Given a time series of embeddings $\{e_1, \dots, e_T\}$ , the LOO medoid score at time $t_k$ compares $e_k$ to a representative embedding computed from the remaining observations:

$E_{-k} = \{e_j \mid j \neq k\}.$

We define that representative as the medoid of $E_{-k}$ , i.e. the embedding that minimizes total distance to all others:

$e_{\text{medoid}}^{(-k)} = \arg\min_{e_i \in E_{-k}} \sum_{e_j \in E_{-k}} d(e_i, e_j).$

The LOO medoid score is then: $\mathrm{LOO}_{\mathrm{medoid}}(t_k) = d\left(e_k, e_{\text{medoid}}^{(-k)}\right),$

where $d(\cdot, \cdot)$ is a distance metric (e.g., Euclidean or cosine distance).

Interpretation: $\mathrm{LOO}_{\mathrm{medoid}}(t_k)$ measures how far the embedding at time $t_k$ is from the medoid of the remaining time series. Larger values indicate that the year-specific embedding is less consistent with the rest of the temporal sequence, suggesting a potential change or anomalous state with respect to all other years.

Results

Using the workflow above, we compute two normalized change layers across the nine-year Colorado subset of the AEF mosaic: a difference-over-time score and a LOO medoid score.

The intermediate result below is a Zarr dataset. You may notice:

9 time steps, corresponding to the annual AEF embeddings from 2017 through 2025
2 score bands, one for consecutive difference and one for LOO medoid deviation
~235Gb of data

ds = xr.open_zarr(
    "s3://sandbox-wherobots-mosaics-tmp/rasterflow/rasterflow/development/rbz2h49wvpjt7h4nzq5g/76e986b947cc.zarr"
)
ds

<xarray.Dataset> Size: 253GB
Dimensions:     (time: 9, band: 2, y: 44800, x: 78336)
Coordinates:
  * time        (time) int32 36B 2017 2018 2019 2020 2021 2022 2023 2024 2025
  * band        (band) object 16B 'difference_alpha_earth_cosine_distance_rob...
  * y           (y) float64 358kB 41.0 41.0 41.0 41.0 ... 36.98 36.98 36.98
  * x           (x) float64 627kB -109.1 -109.1 -109.1 ... -102.0 -102.0 -102.0
Data variables:
    embeddings  (time, band, y, x) float32 253GB dask.array<chunksize=(1, 1, 256, 256), meta=np.ndarray>
Attributes: (12/15)
    proj:code:               EPSG:4326
    spatial:dimensions:      ['y', 'x']
    spatial:transform:       [8.983111749910169e-05, 0.0, -180.0, 0.0, -8.983...
    spatial:transform_type:  affine
    spatial:bbox:            [-180.0, -83.36280346631479, 180.2213438735178, ...
    spatial:shape:           [1859584, 4009984]
    ...                      ...
    geoemb:model:            https://developers.google.com/earth-engine/datas...
    geoemb:source_data:      https://source.coop/tge-labs/aef/v1/annual/
    geoemb:data_type:        int8
    geoemb:gsd:              8.983111749910169e-05
    geoemb:quantization:     {'method': 'signed_square', 'original_dtype': 'f...
    zarr_conventions:        [{'uuid': 'f17cb550-5864-4468-aeb7-f3180cfb622f'...

Visualization

To make the result easier to explore, we run the build_zarr_multiscales function to create a multiscale Zarr mosaic optimized for interactive viewing and analysis, then publish it to a public bucket for anyone to explore.

xr.open_datatree("s3://wherobots-examples/rasterflow/mosaics/co-aef-change.zarr")

Using the multiscales convention and building on carbonplan’s zarr-layer, we can visualize the nine-year mosaic.

If the viewer below is not working, try fullscreen mode here.

Examples of Changes in Boulder, Colorado

The video below shows the normalized change scores around Boulder, Colorado. It is a quick qualitative pass rather than a formal validation exercise, but several patterns stand out more clearly once the scores are viewed over time.

Takeaways and Next Steps

This follow-up suggests that leave-one-out medoid scoring, combined with robust time-wise normalization, is more useful for inspecting AEF change signal than raw distance alone. The normalization makes change magnitude more comparable across pixels, while the leave-one-out medoid highlights anomalous signal across the full observation window rather than only between adjacent time steps. This does not make AEF embeddings a complete change-detection product, but it does provide a clearer and more consistent way to identify which years look unusual at each pixel.

That matters because AEF operates at a spatial and temporal resolution where many changes are subtle, cumulative, or mixed with recurring seasonal variation. A normalized score does not solve interpretation by itself, but it makes the first-pass search for interesting change much more practical.

In this notebook, we used run_mosaics_change on the RasterFlow Remote client to compute change scores from the global AEF mosaic, built a multiscales derivative for browser-based exploration, and reviewed several examples around Boulder. Followup work may include using the embeddings to also help filter the object types of interest for taking yet another step closer towards a more complete change-detection product, especially where we have external labels or stronger domain priors.

If you are interested in this approach, reach out to the Wherobots team at support@wherobots.com.

For a complementary approach aggregating embeddings at the field level, see Zonal Statistics on AlphaEarth Embeddings.

Get Started with RasterFlow

START BUILDING

Key takeaways

Raw AlphaEarth Foundations embedding distances do not live on a consistent scale across pixels or land-cover types. This follow-up scores each annual embedding by how much it stands out from the rest of its local time series using leave-one-out (LOO) medoid scores in embedding space, then robust z-score scaling through time per pixel.
The workflow reads the global AEF Zarr mosaic on Source Cooperative: annual 64-band embeddings at 10-meter resolution in EPSG:4326, nine snapshots from 2017 through 2025. The global store is 3.81 PiB with about 1.02 billion chunks; the demo subsets Colorado rather than the full mosaic, and only land-intersecting chunks are present.
RasterFlow's run_mosaics_change computed both difference-over-time and LOO medoid scores with AlphaEarth cosine distance and robust_time_zscore=True. The Colorado result is a ~235 GB Zarr (9 time steps × 2 score bands). build_zarr_multiscales produced a public multiscale mosaic at s3://wherobots-examples/rasterflow/mosaics/co-aef-change.zarr.
The post is a practical workflow, not a definitive change-detection benchmark. LOO medoid plus robust time-wise normalization makes change magnitude more comparable across pixels and highlights anomalous years across the full nine-year window, not only between adjacent years. Examples around Boulder are qualitative.

AlphaEarth Embeddings, Zonal Statistics, and PCA

Posted on May 6, 2026October 4, 2026 by Ben Pruden

Author: Len Strnad

AlphaEarth Foundations embeddings can be useful at the pixel level when the task is tightly local, such as inspecting boundaries, fine-scale heterogeneity, or small spatial disturbances. But many practical workflows need signals that stay meaningful after aggregation, whether to fields, parcels, patches, or other management units. Zonal statistics provide a simple way to test that idea: aggregate embeddings within field polygons and ask whether the field-level vectors still preserve interpretable structure.

In this notebook, we compute mean AlphaEarth Foundations (AEF) embeddings for about 3.5 million Iowa fields across 2023-2025, then project those field-level vectors into RGB using the first three PCA components. The result is a GeoParquet dataset we can inspect spatially and compare across fields.

The main takeaway is straightforward: the aggregated embeddings form visible clusters that may correspond to crop type or management practice. That does not prove a downstream classifier will work, but it is encouraging evidence that zonal AEF summaries could support lightweight field-level analysis with simple models rather than per-pixel inference.

RasterFlow Tools Used

build_and_predict_mosaic_recipe on the RasterFlow Remote client for building a mosaic and running inference using the Fields of the World model.
vectorize_mosaic on the RasterFlow Remote client for converting field predictions from raster mosaics into GeoParquet.
An experimental zonal-statistics workflow for aggregating AEF embeddings over field geometries.
An experimental PCA workflow for appending the first three components for visualization.

Global Zarr AEF Mosaic

For this analysis, we use the global AlphaEarth Foundations Zarr mosaic hosted on Source Cooperative. Thanks to Taylor Geospatial and Source Cooperative for publishing the dataset in a form that can be accessed lazily with standard Python tooling.

The mosaic stores annual AEF embeddings at 10-meter resolution in EPSG:4326, with nine temporal snapshots from 2017 through 2025. Working from one global store keeps the workflow simple: we can subset Iowa, compute field-level summaries, and publish a derived product without assembling local inputs.

import xarray as xr

ds = xr.open_zarr("s3://us-west-2.opendata.source.coop/tge-labs/aef-mosaic/", consolidated=False)
ds

<xarray.Dataset> Size: 4PB
Dimensions:     (time: 9, band: 64, y: 1859584, x: 4009984)
Coordinates:
  * time        (time) int32 36B 2017 2018 2019 2020 2021 2022 2023 2024 2025
  * band        (band) object 512B 'A00' 'A01' 'A02' 'A03' ... 'A61' 'A62' 'A63'
  * y           (y) float64 15MB 83.69 83.69 83.69 ... -83.36 -83.36 -83.36
  * x           (x) float64 32MB -180.0 -180.0 -180.0 ... 180.2 180.2 180.2
Data variables:
    embeddings  (time, band, y, x) int8 4PB dask.array<chunksize=(1, 64, 256, 256), meta=np.ndarray>
Attributes: (12/15)
    proj:code:               EPSG:4326
    spatial:dimensions:      ['y', 'x']
    spatial:transform:       [8.983111749910169e-05, 0.0, -180.0, 0.0, -8.983...
    spatial:transform_type:  affine
    spatial:bbox:            [-180.0, -83.36280346631479, 180.2213438735178, ...
    spatial:shape:           [1859584, 4009984]
    ...                      ...
    geoemb:model:            https://developers.google.com/earth-engine/datas...
    geoemb:source_data:      https://source.coop/tge-labs/aef/v1/annual/
    geoemb:data_type:        int8
    geoemb:gsd:              8.983111749910169e-05
    geoemb:quantization:     {'method': 'signed_square', 'original_dtype': 'f...
    zarr_conventions:        [{'uuid': 'f17cb550-5864-4468-aeb7-f3180cfb622f'...

At first glance, 3.81 PiB and 1,024,049,664 chunks in S3 sound impractically large. For this demo, a few details keep the workflow manageable:

only land-intersecting chunks are present
the dataset is sharded into nine spatial chunks, which lowers the effective object count
we subset Iowa rather than operating on the full global mosaic

This mosaic is the same one used in Change Detection Using AlphaEarth Foundations (Part 2).

Zonal Statistics Experimental Workflow

We use rayzon, a side project I built to compute zonal statistics using Ray.

Without tuning, it computes zonal statistics for about 3 million Iowa fields across three years in roughly five minutes on a Ray cluster with eight workers and 16 cores per worker.

It also supports t-digest-backed quantile estimation for statistics beyond mean, min, and max, though that path is still experimental. If this workflow is relevant to your team, let us know so we can prioritize productizing this type of analysis in RasterFlow.

Zonal Statistics Results

Using that workflow, we compute mean zonal statistics for Iowa fields across three years of the global AEF Zarr mosaic.

The resulting GeoParquet is published at s3://wherobots-examples/rasterflow/vectors/iowa_aef_zonal_stats/ for anyone who wants to inspect it.

import pyarrow.parquet as pq

ds = pq.ParquetDataset("s3://wherobots-examples/rasterflow/vectors/iowa_aef_zonal_stats/")
ds.schema

geometry: binary
score_mean: float
score_std: float
source_store: string
local_x_offset: int64
local_y_offset: int64
global_x_offset: int64
global_y_offset: int64
path: string
feature_id: int64
bbox: struct<xmax: double, xmin: double, ymax: double, ymin: double>
  child 0, xmax: double
  child 1, xmin: double
  child 2, ymax: double
  child 3, ymin: double
mean: list<element: double>
  child 0, element: double
pca_components: list<element: float>
  child 0, element: float
hilbert_key: int64
time: dictionary<values=string, indices=int32, ordered=0>
label: dictionary<values=string, indices=int32, ordered=0>
-- schema metadata --
geo: '{"version": "1.1.0", "primary_column": "geometry", "columns": {"geo' + 1808

The schema exposes three columns of interest:

geometry: the geometry of each field, with one row per field
mean: a 64-element list containing the mean value for each AEF embedding dimension within the field
pca_components: a 3-element list containing the first three PCA components of each mean embedding, used here as RGB channels for visualization

We also partition the dataset using a Hive-style layout on time and label. In this example, label is always field; the model also produces field_boundaries and non_field, but those classes are not the focus here.

# what are the partitioning keys?
ds.partitioning.dictionaries

[<pyarrow.lib.StringArray object at 0x14199d660>
 [
   "2023-01-01",
   "2024-01-01",
   "2025-01-01"
 ],
 <pyarrow.lib.StringArray object at 0x14199d1e0>
 [
   "field"
 ]]

The output contains about 3.5 million Iowa fields across the 2023-2025 period.

sum([fragment.count_rows() for fragment in ds.fragments])

Visualizing Confidences

The dataset also includes score_mean, the mean model confidence over each polygonized field. Because fields are formed from contiguous pixels above the vectorization threshold, this gives us a compact way to inspect how confident the model was within each extracted geometry. Here we use the 0.5 threshold from vectorize_mosaic.

If the viewer below is not working, try fullscreen mode here.

We can see a great deal of variance in the confidence scores between fields, which can provide some critical feedback for modeling efforts and help guide stratefied sampling for well-balaned training sets.

Visualizing Zonal Statistics of Mean Embeddings with PCA

We also visualize each field’s mean AEF embedding by mapping the first three PCA components to RGB in GeoParquet.

If the viewer below is not working, try fullscreen mode here.

The map shows clear color clustering, especially purples, yellow-greens, and blues. To check whether that spatial pattern is reflected numerically, we next plot the first two PCA components in a scatter plot.

import pandas as pd
import matplotlib.pyplot as plt

df = pd.read_parquet(
    "s3://wherobots-examples/rasterflow/vectors/iowa_aef_zonal_stats/",
    columns=["pca_components"],
)

sample = df.sample(n=20000)
sample["pca_1"] = sample["pca_components"].apply(lambda x: x[0])
sample["pca_2"] = sample["pca_components"].apply(lambda x: x[1])
sample["pca_3"] = sample["pca_components"].apply(lambda x: x[2])

sample.plot.scatter(
    x="pca_1", y="pca_2", c=sample["pca_components"].tolist(), alpha=0.5, figsize=(10, 8)
)

# overlay 2d density plot
plt.hexbin(sample["pca_1"], sample["pca_2"], gridsize=50, cmap='Greys', alpha=0.2, marginals=True)

# add normalized density to 
plt.title("PCA Components Colored by (3) PCA Components")
plt.xlabel("PCA Component 1")
plt.ylabel("PCA Component 2")
plt.show()

We see three broad clusters in the 2D PCA scatter plot, and the cluster colors align with the RGB visualization above. The purple and green groups appear tighter than the yellow group, but the overall separation is still encouraging. That suggests a simple linear probe on top of the AEF embeddings may be able to distinguish field types or management patterns without manual feature engineering or complex models.

Beyond classification, AEF may also be useful for stratified sampling, supervised training set design, or in situ analysis of field management practices.

Takeaways and Next Steps

Aggregated AEF embeddings remain informative at field scale, not just pixel scale.
The zonal-statistics workflow scales to millions of fields and hundreds of millions of pixels on a modest Ray cluster.
GeoParquet and a single global Zarr mosaic work well together for this style of analysis.
Field-level confidence scores provide an additional quality signal alongside the embeddings.
The PCA visualization shows clear clustering that may relate to crop type or management practice and is worth testing with labeled data.
For change detection applications using the same embeddings, see Change Detection Using AlphaEarth Foundations (Part 2) and How Agricultural Fields Change in AlphaEarth Foundations.

If you are interested in this approach, reach out to the Wherobots team at support@wherobots.com.

Get Started

Key takeaways

Mean AlphaEarth Foundations embeddings were computed for about 3.5 million Iowa fields (3,516,044 rows) across 2023–2025, then projected to RGB using the first three PCA components. The result is a GeoParquet dataset at s3://wherobots-examples/rasterflow/vectors/iowa_aef_zonal_stats/.
The experimental zonal-statistics workflow (rayzon) computed zonal stats for about 3 million Iowa fields across three years in roughly five minutes on a Ray cluster with eight workers and 16 cores per worker, without tuning.
Each row has a field geometry, a 64-element mean embedding, and a 3-element pca_components list used as RGB. The dataset is Hive-partitioned on time (2023-01-01, 2024-01-01, 2025-01-01) and label (field). Field polygons come from Fields of the World predictions vectorized at a 0.5 threshold.
PCA scatter of 20,000 sampled fields shows three broad clusters (purples, yellow-greens, blues) that may correspond to crop type or management practice. That is encouraging evidence for lightweight field-level analysis with simple models, not proof a downstream classifier will work.

Detecting Objects From Text Prompts with RasterFlow and Segment Anything 3

Posted on May 4, 2026October 3, 2026 by Ryan Avery

RasterFlow now Supports Detections from Text!

RasterFlow now makes it simple to run promptable geospatial vision models across large aerial and satellite imagery collections, removing the need to build bespoke inference pipelines. With our new SAM 3 support, you can prompt for concepts like “roofs”, “roads”, or “shipping containers” and turn those detections into vector outputs ready for analysis in Wherobots, from city scale to country scale.

Here we see SAM 3 localize and capture roofs really well. It also segments roads pretty cleanly, though there is still room for improvement. Results are shown on the inference input, a NAIP 30 centimeter basemap, and include all instance pixels above 50% confidence.

In this post, we put SAM 3 to the test detecting roofs in suburban neighborhoods, shipping containers in crowded loading docks, and tractors across agricultural landscapes. Throughout, we comment on where it succeeds, where it falls short, and how results compare to previous SAM models. Finally, we discuss what’s next for computer vision in Earth observation (EO) and how the community can build on models like SAM to create more promptable, flexible applications that handle the scale and diversity of remote sensing imagery.

How Segment Anything Model 3 Improves on SAM 1 and SAM 2 for Remote Sensing

In 2023, the Segment Anything Model (SAM) set a new paradigm for computer vision. While most models were trained to address one task at a time, SAM demonstrated that a single model could reliably classify, localize, and segment many kinds of objects in complex natural scenes.

SAM set a new standard task for computer vision models that went beyond classification or semantic segmentation: Promptable Visual Segmentation (PVS), which takes spatial prompts such as boxes, points, or masks as input and predicts masks.

However, SAM 1’s performance on out-of-distribution imagery domains had issues. Back in 2023, while I was working on detection projects with 10 meter resolution imagery, I found it often took multiple rounds of prompting to get useful segments, and many results still needed manual cleanup.

For example, in a rural town, SAM 1 can segment many features across different spatial scales, but it still misses all roads.

SAM 1 segments many features across different spatial scales in Dikwa, Nigeria, but still misses roads.

SAM 2 improved accuracy and added video support, but both SAM 1 and SAM 2 were limited to making predictions without associated labels. The masks they produced were not grounded in categories useful for deriving insights.

Fast forward to today: SAM 3 addresses Promptable Concept Segmentation (PCS), a task where a model can accept either spatial prompts (masks, boxes, pixels) or text prompts like “cat”, “dog”, even “roofs”!

The SAM 3 architecture supports Promptable Concept Segmentation, accepting text and image exemplars in addition to spatial prompts to produce masks with labeled concepts. It also benefits from pretraining on overhead imagery. Source: YouTube.

This opens up new possibilities for detecting objects in imagery. Because SAM 3’s training data spans many imaging domains, including overhead aerial imagery, it can succeed in many Earth observation contexts where previous model generations struggled. For example, one can use SAM 3 to predict “roofs” in NAIP imagery simply by asking, with no model training or ad hoc labeling needed, and get pretty stellar results.

With RasterFlow, we can create imagery mosaics and predict roofs simply by prompting SAM 3 with short noun phrases like “roofs”.

To put SAM 3 to the test on Earth observation imagery, we generated many more SAM 3 predictions on top of National Agriculture Program 30 centimeter imagery using RasterFlow, our scalable mosaic building and inference engine. Below details the performance of RasterFlow on this high resolution detection task.

RasterFlow Task	Billed Runtime	Output	Total Output Size	Spatial Scale of Output
build_gti_mosaic	11m 47s	NAIP Zarr mosaic	133.08 GB (123.94 GiB)	~114.15 km × 66.81 km
predict_mosaic_geometries	39m 51s	GeoParquet geometry results	263.85 MB (251.63 MiB)	~111.02 km × 66.81 km

See an AI agent take plain-language question and orchestrating pipelines to end results using RasterFlow and the Wherobots MCP server, check out the upcoming live session.

See it in action

Join us live

How to Run SAM 3 on Aerial and Satellite Imagery with RasterFlow

What used to take a complex mix of imagery ETL, bespoke inference pipelines, and self-provisioned infrastructure for large Earth observation processing now takes under an hour and two simple Python functions. Let’s check it out below.

Build an Aerial Imagery Mosaic from NAIP

First, we will generate a mosaic: a seamless, stitched-together image from many independent remote sensing scenes. RasterFlow handles all the data sourcing, loading, cleaning, and partitioning into a data asset optimized for inference for a particular model.

from rasterflow_remote import RasterflowClient

rf_client = RasterflowClient()
mosaic_output = rf_client.build_gti_mosaic(
    gti="s3://wherobots-examples/rasterflow/indexes/naip_index.parquet",
    aoi="s3://wherobots-examples/rasterflow/aois/marion_county.parquet",
    bands=["red", "green", "blue", "nir"],
    location_field="url",
    crs_epsg=3857,
    time_column="year",
    skip_xy_coords=False,
    xy_chunksize=1024,
    query="res == 0.3 and time >= '2022-01-01' and time <= '2023-01-01'",
    requester_pays=True,
    sort_field="time",
)
print(mosaic_output.uri)

Run the SAM 3 Inference on the Mosaic

With our mosaic built, we can now run inference with SAM 3. Unlike more rigid models which only predict one category, SAM 3 can accept one or more text prompts and detect all matching objects in a single pass. The runtime for a given batch of imagery scales linearly with the number of detections. More prompts tends to mean more detections, so keep in mind that the more prompts you add, the longer you can expect an inference run to take.

from rasterflow_remote.data_models import GeometryActorEnum, MergeModeEnum

model_output = rf_client.predict_mosaic_geometries(
    store="s3://wherobots-examples/rasterflow/mosaics/marion_county.zarr",
    model_path="https://huggingface.co/wherobots/sam3-text-geometry-pt2/resolve/main/full_sam3_pipeline.pt2",
    patch_size=1008,
    clip_size=0,
    device="cuda",
    features=["red", "green", "blue"],
    labels=["roads", "airplanes", "airports", "roofs", "solar panels", "swimming pools", "shipping containers", "tractors"],
    actor=GeometryActorEnum.TEXT_TO_VECTOR_GEOMETRIES,
    max_batch_size=1,
    confidence_threshold=0.5,
    merge_mode=MergeModeEnum.NONE,
    xy_block_multiplier=1,
)

Analyze the SAM 3 Detections with WherobotsDB

We can visualize both the mosaic and detections directly in the notebook with an embedded RasterFlow map. The RGB imagery mosaic and the SAM 3 detections have been web-optimized with RasterFlow for fast and fluid browsing.

If the embed does not load in your environment, open it in a new tab here.

After visually exploring detections, we can load the results into WherobotsDB for quantitative analysis. WherobotsDB lets us post-process geometries, calculate zonal statistics on other rasters, count objects, measure clustering, and much more.

import os
from sedona.spark import *
from pyspark.sql.functions import *
from wherobots import vtiles

config = (
    SedonaContext.builder()
    .getOrCreate()
)
sedona = SedonaContext.create(config)

parquet_path = "s3://wherobots-examples/rasterflow/model-outputs/marion_county_sam3/"
df = sedona.read.format("geoparquet").load(parquet_path)

df.printSchema()

df.show(10)

Here we’ll use ST_Area to estimate the total square km area of all roofs in Marion County.

roofs_df = df.filter(col("label") == "roofs")
source_srid = roofs_df.selectExpr("ST_SRID(geometry) AS srid").first()["srid"]
source_crs = f"EPSG:{source_srid}"
print(source_crs)

# Run the SRID inspection cell above first so source_crs reflects the data's actual CRS.
roof_areas = roofs_df.withColumn(
    "area_sq_m",
    expr(f"ST_AreaSpheroid(ST_Transform(geometry, '{source_crs}', 'EPSG:4326'))"),
).cache()

roof_areas_agg = roof_areas.agg(
    sum("area_sq_m").alias("total_area_sq_m"),
    (sum("area_sq_m") / lit(1_000_000.0)).alias("total_area_sq_km"),
)

roof_areas_agg.show(truncate=False)

roof_size_stats = roof_areas.agg(
    count("*").alias("roof_count"),
    avg("area_sq_m").alias("avg_roof_sq_m"),
    expr("percentile_approx(area_sq_m, 0.5)").alias("median_roof_sq_m"),
    (avg("area_sq_m") * lit(10.7639)).alias("avg_roof_sq_ft"),
)

stats = roof_size_stats.first()

print(f"roof_count: {stats['roof_count']:,}")
print(f"avg_roof_sq_m: {stats['avg_roof_sq_m']:.2f}")
print(f"median_roof_sq_m: {stats['median_roof_sq_m']:.2f}")

label_counts = df.groupBy("label").count().orderBy(col("count").desc(), col("label"))
label_counts.show(truncate=False)

We can also generate web-optimized vectors in PMTiles format. PMTiles are easily shareable and plug directly into geospatial visualization applications, making them a great choice for distributing detection results.

df_tiles = df.withColumn("layer", col("label"))

output_path = 's3://wherobots-examples/rasterflow/model-outputs/marion_county_sam3.pmtiles'
vtiles.generate_pmtiles(df_tiles, output_path)

SAM 3 Accuracy on Earth Observation: Where It Succeeds and Where It Falls Short

After experimenting with SAM 3 on NAIP, I’m impressed with the range of categories that SAM 3 can positively identify. At the same time, the model still has clear precision and recall gaps for more niche semantic categories, like “tractors” or “shipping containers”.

SAM 3 correctly detects tractors in a harvested field.

In an agricultural context, SAM 3 finds some tractors but also misses some.

SAM 3 correctly identifies shipping containers in a storage yard.

False positives: SAM 3 labels rectangular roofs as “shipping containers”.

These examples illustrate the gap between SAM 3’s strengths on common categories and its current limitations on more niche ones. I’m bullish that if SAM 3 were fine-tuned on a larger corpus of mixed high-resolution imagery and high-quality labels, it would perform even better outside of its primary domain of natural imagery.

Even so, I think the potential for SAM 3 for simpler categories like “roofs” is underutilized. SAM 3 could potentially improve many existing datasets we rely on to make decisions, like Overture Buildings and detections of other kinds of structures.

What to Try Next with SAM 3 on Wherobots

Some things we didn’t showcase but I recommend trying in your experiments:

Can you find and count vehicles that are a certain color?
Can you use SAM 3 to segment city streets? What about rural roads?
Try to map forest canopy! SAM 3 can do a pretty great job at segmenting concepts that don’t map cleanly to individual objects, but are more textural.
For a larger-scale example of RasterFlow in production, see how RasterFlow processed 2.3B Overture Maps features for Fields of the World.

If you’re excited to try SAM 3 on Wherobots, sign up for RasterFlow Private Preview. You can also get in touch at ryan@wherobots.com or talk to us here. I’d love to hear about what you’re looking to build and how SAM 3 could fit into your detection workflows.

Try SAM3 on Wherobots

START BUILDING

Key takeaways

RasterFlow now supports SAM 3 promptable concept segmentation: text prompts such as ‘roofs’, ‘roads’, or ‘shipping containers’ produce vector outputs from aerial and satellite mosaics. SAM 3 accepts spatial prompts or text, and it was pretrained on overhead imagery, unlike SAM 1/2.
On Marion County NAIP 30 cm, build_gti_mosaic billed 11 minutes 47 seconds and wrote a 133.08 GB Zarr mosaic (~114.15 km × 66.81 km). predict_mosaic_geometries billed 39 minutes 51 seconds and wrote 263.85 MB of GeoParquet. Eight labels were prompted in one pass at confidence_threshold 0.5.
SAM 3 localizes roofs well on the NAIP basemap and segments roads reasonably, with room for improvement. Niche classes such as tractors and shipping containers show clear precision and recall gaps, including rectangular roofs labeled as shipping containers.
What used to take a mix of imagery ETL, bespoke inference pipelines, and self-provisioned infrastructure now takes under an hour and two Python functions: build_gti_mosaic then predict_mosaic_geometries. Detections load into WherobotsDB as GeoParquet for ST_AreaSpheroid, counts, and PMTiles via vtiles.generate_pmtiles.

How Agricultural Fields Change in AlphaEarth Foundations

Posted on April 7, 2026October 4, 2026 by Ben Pruden

Author: Len Strnad

AlphaEarth Foundations is a geospatial AI model from Google DeepMind. It compresses a year of satellite observations into 64-dimensional embedding vectors for every 10-meter pixel on Earth’s land surfaces and shallow coastal waters annually. In this notebook, we use AEF to examine how agricultural fields evolve across several regions and ask a practical question: what kinds of field-level change become visible when we inspect the embeddings directly rather than relying only on imagery?

We begin by defining a few agricultural areas of interest and building Zarr mosaics over each one from COGs available in Source Cooperative. From there, we compare several views of temporal change: the first three embedding bands rendered as RGB, a Principal Component Analysis (PCA) projection that makes broader structure easier to interpret, and an experimental distance-over-time workflow that highlights where embedding values remain stable or shift over time.

The goal is not to claim a single definitive interpretation of AEF. Instead, this notebook shows a practical workflow for exploring seasonal cycles, repeated interventions, and other field-scale dynamics in embedding space.

RasterFlow Tools Used

RasterFlow is Wherobots’ serverless inference engine for Earth Observation data. It builds inference-ready mosaics from satellite imagery and runs distributed workflows at scale. This notebook uses three RasterFlow capabilities:

build_gti_mosaic for constructing Zarr mosaics from AlphaEarth Foundations data
an experimental PCA workflow for dimensionality reduction
and an experimental distance-over-time workflow for measuring embedding change between timesteps.

Selecting Agricultural AOIs Across Iowa, California, and Japan

We start with three small AOIs that contain agricultural fields in different regions. This gives us a compact but varied set of examples for exploring how AlphaEarth Foundations embeddings change over time.

import geopandas as gpd
from shapely.geometry import box
from folium import Map

iowa1 = [-94.732361, 41.911475, -94.237289, 42.112996]
california1 = [-120.11095, 36.63082, -120.018082, 36.683839]
japan1 = [141.018019, 39.226668, 141.067801, 39.249803]
aois = [iowa1, california1, japan1]

gdf = gpd.GeoDataFrame({"geometry": [box(*aoi) for aoi in aois]}, crs="EPSG:4326")

# start on the first AOI
m = Map(location=gdf.iloc[0].geometry.centroid.coords[0][::-1], zoom_start=8)
gdf.explore(height="120%", width="120%", m=m)

Make this Notebook Trusted to load map: File -> Trust Notebook

# store this out to parameterize which Zarr mosaics we want to build for AEF
gdf.to_parquet("s3://union-sandbox-unionai-wherobots/aois/field_change_aef.parquet")

AlphaEarth Foundations Zarr Mosaics by AOI

The next step is to create an AlphaEarth Foundations embeddings Zarr mosaic for each AOI. This gives us a time-indexed raster store that we can query efficiently for visualization and downstream analysis. The output here is from the build_gti_mosaic RasterFlow workflow.

import geopandas as gpd

mosaic_index = gpd.read_parquet(
    "s3://union-sandbox-unionai-wherobots/ha/wherobots/wherobots-mosaics/development/rnmkrr4sbxsdf24t288p/d8cq2y2k3b1r4p9gr0k5n1ldh/1/7q/fdjxpi5q/61a6e8f637cce7ddf177ff7a6637aceb/"
)
mosaic_index

	geometry	location
0	POLYGON ((-94.23729 41.91148, -94.23729 42.113…	s3://sandbox-wherobots-mosaics-tmp/mosaics/rnm…
1	POLYGON ((-120.01808 36.63082, -120.01808 36.6…	s3://sandbox-wherobots-mosaics-tmp/mosaics/rnm…
2	POLYGON ((141.0678 39.22667, 141.0678 39.2498,…	s3://sandbox-wherobots-mosaics-tmp/mosaics/rnm…

The result is one Zarr mosaic per AOI, which we can load independently for inspection and comparison.

Inspect an Example AEF Zarr Mosaic

import xarray as xr

ds = xr.open_zarr(mosaic_index.iloc[0]["location"])
ds

<xarray.Dataset> Size: 38GB
Dimensions:      (time: 9, band: 64, y: 3020, x: 5512)
Coordinates:
  * time         (time) datetime64[ns] 72B 2017-01-01 2018-01-01 ... 2025-01-01
  * band         (band) object 512B 'band_0' 'band_1' ... 'band_62' 'band_63'
  * y            (y) float64 24kB 5.178e+06 5.178e+06 ... 5.148e+06 5.148e+06
  * x            (x) float64 44kB -1.055e+07 -1.055e+07 ... -1.049e+07
    spatial_ref  int64 8B ...
Data variables:
    variables    (time, band, y, x) float32 38GB dask.array<chunksize=(1, 1, 1024, 1024), meta=np.ndarray>

# pull this into memory once to reuse
subset = ds["variables"][:, :3, :, :].compute()

import imageio.v3 as iio
import numpy as np


def to_uint8_robust(array: xr.DataArray, low: xr.DataArray, high: xr.DataArray) -> np.ndarray:
    arr = array.transpose("y", "x", ...).astype("float32")
    scale = np.where((high - low) == 0, 1.0, (high - low))
    arr = (arr - low) / scale
    arr = np.clip(arr, 0, 1)
    return (arr * 255).astype(np.uint8)


def robust_params(da: xr.DataArray) -> tuple[xr.DataArray, xr.DataArray]:
    q_low, q_high = 1, 99
    lo = da.quantile(q_low / 100, dim=("time", "y", "x"))
    hi = da.quantile(q_high / 100, dim=("time", "y", "x"))
    return lo, hi

low, high = robust_params(subset)
frames = [to_uint8_robust(subset.sel(time=t), low, high) for t in subset.time]
iio.imwrite("assets/rgb_animation.gif", frames, duration=1000, loop=0)

To establish a baseline view, we render the first three embedding bands as an RGB animation over time. The color mapping is not directly semantic, but it provides a quick way to spot recurring seasonal patterns and abrupt changes.

Using PCA to Identify Seasonal Patterns in Embeddings

A more interpretable view comes from projecting the 64-dimensional embeddings into their first three principal components. PCA concentrates the dominant variance into three channels, which often makes field-level change easier to see.

# training sample for PCA using just the first year and a spatial slice
sample = ds["variables"][:, :, :1024, :1024].compute()
sample = sample.transpose("band","time","y","x").data.reshape(64, -1)
sample = sample[:, ~np.isnan(sample).any(axis=0)]

from sklearn.decomposition import PCA

model = PCA(n_components=3)
model.fit(sample.T)

PCA(n_components=3)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

To keep the notebook lightweight, we pull all time steps and bands for a small spatial subset rather than the full mosaic. AlphaEarth Foundations embeddings contain 64 bands, so this smaller window reduces data transfer while still preserving the temporal behavior we want to inspect. The same workflow can be scaled up for larger analyses.

# can take 30+ seconds
subset = ds["variables"][:, :, :1024, :2048].compute()

# Transform with PCA
flattened = subset.transpose("band", "time", "y", "x").data.reshape(64, -1).T
transformed = model.transform(flattened)
img = transformed.T.reshape(3, subset.shape[0], subset.shape[2], subset.shape[3]).transpose(
    1, 0, 2, 3
)

# Use the existing xarray data structure to build a nice home for PCA data.
subset_copy = subset.copy()
array = subset_copy.isel(band=slice(0, 3))

# insert the transformed PCA data
array.data = img
array["band"] = ["pca1", "pca2", "pca3"]

# we chunk only to get a better view of the dataset
array.to_dataset(name='variables')

low, high = robust_params(array)
frames = [to_uint8_robust(array.sel(time=t), low, high) for t in array.time]
iio.imwrite("assets/rgb_animation_pca.gif", frames, duration=1000, loop=0)

The PCA animation makes temporal structure easier to see than the raw RGB view. Field boundaries and repeated change patterns stand out more clearly, while nearby urban areas show a different signature of structural change.

Compared with the first three raw bands, the PCA view surfaces a more coherent pattern of cyclical change. In agricultural areas, the alternating shifts between green and brown tones likely correspond to crop rotations, harvest cycles, irrigation, or other recurring management practices.

PCA is still a simplification, and it comes with tradeoffs:

PCA is a linear transformation, so it may miss more complex relationships in the embedding space.
PCA can be sensitive to outliers and to the sample used to fit the model.
PCA adds computational overhead because fitting and applying the transform requires additional matrix operations.

Another option is to stay in the original embedding space and measure distance between embeddings over time. That gives us a direct way to ask how similar or dissimilar each pixel is from one timestep to the next.

Measure Change Directly in Embedding Space

Here we use an experimental RasterFlow workflow to compute distance over time in the original embedding space.

This workflow is not yet publicly available, but it is useful for exploring whether similarity metrics can reveal patterns that PCA may smooth over. If you are interested in this approach, please reach out to the Wherobots team (support@wherobots.com). For a field-level aggregation approach using the same embeddings, see Zonal Statistics on AlphaEarth Embeddings.

index = gpd.read_parquet(
    "s3://union-sandbox-unionai-wherobots/dw/wherobots/wherobots-rasterflow/development/rfzhkxdjv782gtq9xtzz/a2dwp687c37a3kc1z6593pdax/1/ae/famfddly/0fbb239f05eb805bf4128b8be91316ce/"
)

# we look at an exmple mosaic for now
ex_store = index.iloc[0]["location"]
ds = xr.open_zarr(ex_store).compute()
ds

array = ds["variables"][:, :, :, :]
low = np.nanpercentile(array[0].values, 0)
high = np.nanpercentile(array[0].values, 50)
frames = [to_uint8_robust(array.sel(time=t)[0], low, high)[:, :] for t in array.time]
iio.imwrite("assets/distance_over_time_50p.gif", frames, duration=1000, loop=0)

In the animation above, darker pixels indicate embeddings that remain more similar over time, while lighter pixels indicate larger temporal differences.

Interpretation is still evolving, but the visualization suggests several useful patterns.

Persistent bright pixels may mark locations with repeated or abrupt change. Alternating light-dark cycles may indicate recurring agricultural activity such as planting, harvest, or irrigation.

This does not replace domain knowledge, but it provides another lens for understanding how field-scale practices appear in AlphaEarth Foundations embeddings.

What AlphaEarth Foundations Embeddings Show About Field-Level Dynamics

AlphaEarth Foundations embeddings make it possible to examine agricultural change as a time series in embedding space rather than only as raw imagery. In this notebook, we built Zarr mosaics for several AOIs, visualized temporal variation in the first three embedding bands, projected the embeddings with PCA, and then compared that view with a distance-over-time workflow. Together, these helped draw insight into a few of the many ways AlphaEarth Foundations embeddings (and others) can provide insight into field-level dynamics. Check out Change Detection Using AlphaEarth Foundations (Part 2).

Get Started with RasterFlow

Start Building

Key takeaways

AlphaEarth Foundations (Google DeepMind) compresses a year of satellite observations into 64-dimensional embedding vectors for every 10-meter pixel on Earth's land surfaces and shallow coastal waters, annually. This notebook inspects field-level change in embedding space rather than relying only on imagery.
Three agricultural AOIs are used: Iowa, California, and Japan. RasterFlow build_gti_mosaic builds a Zarr mosaic per AOI from AlphaEarth Foundations COGs on Source Cooperative. An example Iowa mosaic is 38 GB (9 time steps, 64 bands, 2017–2025).
Three views of temporal change: the first three embedding bands as RGB, a PCA projection of all 64 bands into three components, and an experimental cosine-distance-over-time workflow. PCA is described as making field-level seasonal structure easier to see than raw RGB bands.
The goal is not a single definitive interpretation of AEF. It is a practical workflow for exploring seasonal cycles, repeated interventions, and other field-scale dynamics in embedding space, using experimental RasterFlow PCA and distance-over-time tools.

	n_components n_components: int, float or ‘mle’, default=None Number of components to keep. if n_components is not set all components are kept:: n_components == min(n_samples, n_features) If “n_components == ‘mle’“ and “svd_solver == ‘full’“, Minka’s MLE is used to guess the dimension. Use of “n_components == ‘mle’“ will interpret “svd_solver == ‘auto’“ as “svd_solver == ‘full’“. If “0 < n_components < 1“ and “svd_solver == ‘full’“, select the number of components such that the amount of variance that needs to be explained is greater than the percentage specified by n_components. If “svd_solver == ‘arpack’“, the number of components must be strictly less than the minimum of n_features and n_samples. Hence, the None case results in:: n_components == min(n_samples, n_features) – 1	3
	copy copy: bool, default=True If False, data passed to fit are overwritten and running fit(X).transform(X) will not yield the expected results, use fit_transform(X) instead.	True
	whiten whiten: bool, default=False When True (False by default) the `components_` vectors are multiplied by the square root of n_samples and then divided by the singular values to ensure uncorrelated outputs with unit component-wise variances. Whitening will remove some information from the transformed signal (the relative variance scales of the components) but can sometime improve the predictive accuracy of the downstream estimators by making their data respect some hard-wired assumptions.	False
	svd_solver svd_solver: {‘auto’, ‘full’, ‘covariance_eigh’, ‘arpack’, ‘randomized’}, default=’auto’ “auto” : The solver is selected by a default ‘auto’ policy is based on `X.shape` and `n_components`: if the input data has fewer than 1000 features and more than 10 times as many samples, then the “covariance_eigh” solver is used. Otherwise, if the input data is larger than 500×500 and the number of components to extract is lower than 80% of the smallest dimension of the data, then the more efficient “randomized” method is selected. Otherwise the exact “full” SVD is computed and optionally truncated afterwards. “full” : Run exact full SVD calling the standard LAPACK solver via `scipy.linalg.svd` and select the components by postprocessing “covariance_eigh” : Precompute the covariance matrix (on centered data), run a classical eigenvalue decomposition on the covariance matrix typically using LAPACK and select the components by postprocessing. This solver is very efficient for n_samples >> n_features and small n_features. It is, however, not tractable otherwise for large n_features (large memory footprint required to materialize the covariance matrix). Also note that compared to the “full” solver, this solver effectively doubles the condition number and is therefore less numerical stable (e.g. on input data with a large range of singular values). “arpack” : Run SVD truncated to `n_components` calling ARPACK solver via `scipy.sparse.linalg.svds`. It requires strictly `0 < n_components < min(X.shape)` “randomized” : Run randomized SVD by the method of Halko et al. .. versionadded:: 0.18.0 .. versionchanged:: 1.5 Added the ‘covariance_eigh’ solver.	‘auto’
	tol tol: float, default=0.0 Tolerance for singular values computed by svd_solver == ‘arpack’. Must be of range [0.0, infinity). .. versionadded:: 0.18.0	0.0
	iterated_power iterated_power: int or ‘auto’, default=’auto’ Number of iterations for the power method computed by svd_solver == ‘randomized’. Must be of range [0, infinity). .. versionadded:: 0.18.0	‘auto’
	n_oversamples n_oversamples: int, default=10 This parameter is only relevant when `svd_solver=”randomized”`. It corresponds to the additional number of random vectors to sample the range of `X` so as to ensure proper conditioning. See :func:`~sklearn.utils.extmath.randomized_svd` for more details. .. versionadded:: 1.1	10
	power_iteration_normalizer power_iteration_normalizer: {‘auto’, ‘QR’, ‘LU’, ‘none’}, default=’auto’ Power iteration normalizer for randomized SVD solver. Not used by ARPACK. See :func:`~sklearn.utils.extmath.randomized_svd` for more details. .. versionadded:: 1.1	‘auto’
	random_state random_state: int, RandomState instance or None, default=None Used when the ‘arpack’ or ‘randomized’ solvers are used. Pass an int for reproducible results across multiple function calls. See :term:`Glossary `. .. versionadded:: 0.18.0	None