Category: Engineering

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.

-- 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.

Buildings: read the elevation under each footprint

-- 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;

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

-- 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;

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.

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

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.

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.

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

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

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.

-- 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;

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.

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:

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.

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.

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.

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

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.

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.

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.