Example: Low-Cost Roof and Solar Inspection with SAM3

from datetime import datetime
import os

import geopandas as gpd
import wkls

# Area of interest. `or` is a Python keyword, so wkls exposes Oregon by name.
# wkls.us.oregon.cities() and wkls.us.oregon.counties() list what is available,
# and a wrong name raises an error suggesting the closest matches.
AOI = wkls.us.oregon.keizer
AOI_NAME = "keizer"

# The window RasterFlow searches for imagery. Keizer's only 30 cm NAIP is 2022, so
# this window is 2022. Another area needs its own year: step 4 prints the years that work.
START = datetime(2022, 1, 1)
END = datetime(2023, 1, 1)

# What to detect. A list requests several object types in one inference call.
# "roofs", "roof", and "rooftops" are three different prompts and can return
# different detections.
TEXT_PROMPT = ["parking lots", "roads", "roofs", "solar panels"]

# Required, and there is no safe default. At 0.3 primary structures detect reliably
# and small accessory structures are missed.
CONFIDENCE_THRESHOLD = 0.3

# Start with a 30 m² roof floor, but retain small solar arrays for review.
# Step 9 shows what each floor drops before you apply it.
MIN_AREA_M2 = 30
SOLAR_MIN_AREA_M2 = 0

# Where results go. USER_S3_PATH is set for you in a Wherobots notebook.
aoi_path = os.getenv("USER_S3_PATH") + f"{AOI_NAME}.parquet"

import matplotlib.pyplot as plt

gdf = gpd.read_file(AOI.geojson())
gdf.to_parquet(aoi_path)

area_km2 = gdf.to_crs(gdf.estimate_utm_crs()).area.sum() / 1e6
print(f"{AOI_NAME}: {area_km2:,.1f} km2")
print(f"Boundary written to {aoi_path}")

fig, ax = plt.subplots(figsize=(4, 4))
for geom in gdf.geometry:
    # `.geoms` for a MultiPolygon AOI, the geometry itself for a Polygon
    for poly in getattr(geom, "geoms", [geom]):
        x, y = poly.exterior.xy
        ax.plot(x, y, color="black", linewidth=0.8)
ax.set_aspect("equal")
ax.set_title(f"{AOI_NAME}: {area_km2:,.1f} km2")
plt.show()

# One row per NAIP scene, with its footprint, resolution (res), acquisition year and time
naip_index = gpd.read_parquet(
    "s3://wherobots-examples/rasterflow/indexes/naip_index.parquet",
    columns=["geometry", "res", "year", "time"],
)

# Every scene over the AOI, compared in the index's CRS
aoi_geom = gdf.geometry.to_crs(naip_index.crs).union_all()
touching = naip_index.iloc[naip_index.sindex.query(aoi_geom, predicate="intersects")]

# What has been flown here, and at what resolution. Only 0.3 is usable by SAM3.
print("Scenes over this AOI, by acquisition year and resolution:")
print(touching.groupby(["year", "res"]).size())

# The SAM3 recipes run on 30 cm NAIP
MODEL_RES = 0.3

# Scenes at the model's resolution touching the AOI
nearby = touching.query(f"res == {MODEL_RES}")

# Narrowed to the requested date range
covering = nearby.query(f"time >= '{START:%Y-%m-%d}' and time <= '{END:%Y-%m-%d}'")

# The recipe needs the AOI fully inside the matching scenes
if not covering.geometry.union_all().contains(aoi_geom):
    # How much of the AOI is covered, measured in an equal-area CRS
    area_crs = gdf.estimate_utm_crs()
    aoi_area = gdf.geometry.to_crs(area_crs).union_all()
    tiles = covering.geometry.to_crs(area_crs).union_all()
    fraction = max(0.0, 1.0 - aoi_area.difference(tiles).area / aoi_area.area)

    # Years that would work, applying the same full-coverage test as the gate above
    available = sorted(
        year
        for year, scenes in nearby.groupby("year")
        if scenes.geometry.union_all().contains(aoi_geom)
    )
    raise ValueError(
        f"{MODEL_RES:g}m NAIP covers {fraction:.1%} of this AOI between "
        f"{START:%Y-%m-%d} and {END:%Y-%m-%d}; the recipe needs the AOI fully covered. "
        + (
            f"Years with complete {MODEL_RES:g}m coverage over this AOI: {available}. "
            "Set START and END to one of them."
            if available
            else f"No year has complete {MODEL_RES:g}m NAIP coverage over this AOI."
        )
    )

print(f"NAIP coverage is available. Found {len(covering)} NAIP scenes at {MODEL_RES:g}m covering the AOI")

import pandas as pd

# Flight dates of the scenes behind the mosaic
flights = pd.to_datetime(covering["time"])
if flights.dt.tz is not None:
    flights = flights.dt.tz_localize(None)

newest, oldest = flights.max(), flights.min()
today = pd.Timestamp.today().normalize()

print(f"{len(covering)} scenes, flown {oldest.date()} to {newest.date()}")
print(f"Newest imagery: {(today - newest).days / 365.25:.1f} years old")
print(f"Oldest imagery: {(today - oldest).days / 365.25:.1f} years old")
print()
print("Scenes per flight date:")
print(flights.dt.date.value_counts().sort_index())

from rasterflow_remote import RasterflowClient
from rasterflow_remote.data_models import DatasetEnum, GeometryModelRecipes

rf_client = RasterflowClient()

# Preview the same source imagery, AOI, date window, and mosaic CRS
mosaic_index = rf_client.build_mosaics(
    datasets=[DatasetEnum.NAIP_30CM],
    aoi=aoi_path,
    start=START,
    end=END,
    target_crs=3857,
)

# The URI of the store for the first mosaic location, which is what an AOI of a single
# geometry produces. `wkls` returns one geometry per city or county.
mosaic_store = mosaic_index.first_row_mosaic
if mosaic_store is None:
    raise RuntimeError(
        "build_mosaics returned no mosaic store. The index it wrote is "
        f"{mosaic_index.mosaic_index_uri}; check the run in Workload History."
    )

print(f"Mosaic index: {mosaic_index.mosaic_index_uri}")
print(f"Mosaic store: {mosaic_store}")

# Reads the mosaic and writes a second store carrying overview levels and per-band
# histograms. The source mosaic is read, never modified.
optimized_store = rf_client.build_zarr_multiscales(source_store=mosaic_store)

print(f"Source mosaic:   {mosaic_store}")
print(f"Optimized store: {optimized_store.uri}")

from wherobots_gl import Map

# A dark basemap makes the mosaic's own edges visible. Over a satellite basemap a
# gap in the mosaic still looks like imagery.
Map(
    layers=[{"type": "zarr", "source": optimized_store.uri, "name": "NAIP preview mosaic"}],
    basemap="dark",
)

model_output = rf_client.predict_mosaic_geometries_recipe(
    # Path to our AOI in GeoParquet or GeoJSON format
    aoi=aoi_path,

    # Date range for imagery to be used by the model (verified above)
    start=START,
    end=END,

    # CRS for the internal mosaic; SAM3 geometry output is EPSG:4326
    target_crs=3857,

    # The model recipe and text prompt for object detection
    model_recipe=GeometryModelRecipes.SAM3_TEXT_GEOMETRY,
    text_prompt=TEXT_PROMPT,
    confidence_threshold=CONFIDENCE_THRESHOLD,
)

# The URI is None when nothing cleared the confidence threshold
if model_output.uri is None:
    raise RuntimeError(
        "No detections. Lower CONFIDENCE_THRESHOLD, or reword TEXT_PROMPT: "
        "'roofs', 'roof', and 'rooftops' are three different prompts."
    )

detections_gdf = gpd.read_parquet(model_output.uri)
print(f"{len(detections_gdf):,} detections in {detections_gdf.crs.name} (EPSG:{detections_gdf.crs.to_epsg()})")
detections_gdf.head()

print(f"Total detections: {len(detections_gdf)}")
print(f"Columns: {list(detections_gdf.columns)}")
print()
print(detections_gdf.groupby("label").size())
print()
print("Confidence score stats:")
print(detections_gdf.groupby("label")["bbox_score"].describe()[["count", "min", "50%", "max"]].round(3))
print()
print(f"`time` values on the detections: {sorted(pd.to_datetime(detections_gdf['time']).dt.date.unique())}")
print("Flight dates from the scene index:", sorted(flights.dt.date.unique()))

from sedona.spark import *

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

(
    sedona.read.format("geoparquet")
    .load(model_output.uri)
    .createOrReplaceTempView("sam3_detections")
)

# ST_AreaSpheroid needs lon/lat. The recipe returns EPSG:4326, so this is a no-op,
# and it keeps working if a run ever comes back projected.
epsg = detections_gdf.crs.to_epsg() if detections_gdf.crs else 4326
geom = "geometry" if epsg == 4326 else f"ST_Transform(geometry, 'EPSG:{epsg}', 'EPSG:4326')"
floor_sql = f"CASE WHEN label = 'solar panels' THEN {SOLAR_MIN_AREA_M2} ELSE {MIN_AREA_M2} END"
print(f"Detections are EPSG:{epsg}; measuring area on `{geom}`")

# What the floors would drop, before applying them
sedona.sql(f"""
    SELECT
      label,
      COUNT(*)                                                  AS detections,
      ROUND(PERCENTILE_APPROX(ST_AreaSpheroid({geom}), 0.5), 1) AS median_area_m2,
      ROUND(PERCENTILE_APPROX(ST_AreaSpheroid({geom}), 0.9), 1) AS p90_area_m2,
      SUM(CASE WHEN ST_AreaSpheroid({geom}) <= {floor_sql} THEN 1 ELSE 0 END)
                                                                AS below_floor
    FROM sam3_detections
    GROUP BY label
    ORDER BY detections DESC
""").show(truncate=False)

sedona.sql("CREATE DATABASE IF NOT EXISTS examples_temp.sam3_db")

table_name = f"examples_temp.sam3_db.sam3_{AOI_NAME}"

sedona.sql(f"""
    WITH measured AS (
      SELECT
        {geom}                    AS geometry,
        label                     AS layer,
        bbox_score,
        ST_AreaSpheroid({geom})   AS area_m2,
        {floor_sql}              AS min_area_m2
      FROM sam3_detections
    )
    SELECT geometry, layer, bbox_score, area_m2
    FROM measured WHERE area_m2 > min_area_m2
""").writeTo(table_name).createOrReplace()

df_filtered = sedona.table(table_name)
kept, total = df_filtered.count(), sedona.table("sam3_detections").count()

if kept == 0:
    raise RuntimeError(
        f"The per-prompt area floors removed all {total:,} detections. Check the area summary "
        "above: if every median is near zero, the geometries are not lon/lat and "
        "ST_AreaSpheroid cannot measure them."
    )

print(f"{table_name}: kept {kept:,} of {total:,} detections after per-prompt area floors")
df_filtered.groupBy("layer").count().orderBy("layer").show()

from wherobots_gl import Map

# The widget loads layers by URL, so write the filtered detections out before mapping them
filtered_path = os.getenv("USER_S3_PATH") + f"sam3_{AOI_NAME}_filtered.parquet"
df_filtered.write.format("geoparquet").mode("overwrite").save(filtered_path)

Map(
    layers=[{"type": "geoparquet", "source": filtered_path, "name": f"SAM3 detections: {AOI_NAME}"}],
    basemap="satellite",
)

from wherobots import vtiles

full_tiles_path = os.getenv("USER_S3_PATH") + f"sam3_{AOI_NAME}_tiles.pmtiles"
vtiles.generate_pmtiles(df_filtered, full_tiles_path)
print(full_tiles_path)

vtiles.show_pmtiles(full_tiles_path)

# Paste each URI into Layer URL at https://cloud.wherobots.com/map, then Add Layer
print("Layers to add")
print(f"  NAIP preview mosaic (zarr): {optimized_store.uri}")
print(f"  SAM3 PM Tiles (pmtiles):  {full_tiles_path}")

# map_url opens the mosaic on its own; it is None when no viewer URL is available
if optimized_store.map_url:
    print(f"\nMosaic in the map viewer: {optimized_store.map_url}")

(
    sedona.read.format("geoparquet")
    .load("s3://your-bucket/portfolio/policies.parquet")  # your own policy records
    .createOrReplaceTempView("my_portfolio")
)

# The flight date from step 5. The oldest scene is the conservative choice when a
# mosaic spans more than one flight.
imagery_date = oldest.date()

sedona.sql(f"""
    CREATE OR REPLACE TEMP VIEW scored_roofs AS
    SELECT
      p.policy_id,
      r.geometry,
      r.area_m2                                          AS roof_area_m2,
      DATE '{imagery_date}'                              AS imagery_date,
      -- How stale is the evidence?
      months_between(current_date(), DATE '{imagery_date}') / 12
                                                         AS imagery_age_years,
      -- Straight-line consumed life as of today, not as of the flyover.
      -- Replace 25.0 with the expected service life of the covering on that risk.
      greatest(0.0, least(1.0,
        (year(current_date()) - p.roof_installed_year) / 25.0
      ))                                                 AS consumed_life_fraction
    FROM examples_temp.sam3_db.sam3_{AOI_NAME} r
    JOIN my_portfolio p
      ON ST_Intersects(r.geometry, p.geometry)
    WHERE r.layer = 'roofs'
""")

sedona.sql("""
    SELECT * FROM scored_roofs
    ORDER BY consumed_life_fraction DESC, imagery_age_years DESC
""").show()