HERE Technologies and Wherobots announce partnership Learn More

Detecting roads with RasterFlow

This notebook will guide you through detecting roads from aerial imagery, using Wherobots RasterFlow and the ChesapeakeRSC model. You will gain a hands-on understanding of how to run models like ChesapeakeRSC on your selected area of interest, vectorize the model outputs, and work with those vectors in WherobotsDB.

ChesapeakeRSC

The ChesapeakeRSC1 model is an open source segmentation model that can detect roads from high resolution imagery, even when those roads may be occluded by tree cover.

This model predicts 2 classes:

  • background
  • road

It was trained on high resolution imagery from the National Agriculture Imagery Program (NAIP) and labeled land cover data for the state of Maryland in the United States.

Preview: model inputs and outputs

The interactive map linked below shows the input imagery and model outputs for this notebook’s example AOI. Toggle layers in the sidebar to compare the inputs, the raster model output, and the vectorized geometries side-by-side.

Layers (top layer first, as in the map’s layer list):

  • ChesapeakeRSC PM Tiles: Vectorized road polygons in PMTiles format, created from the vector geometries
  • ChesapeakeRSC vector geometries: Vectorized road polygons derived from the raster output
  • ChesapeakeRSC model outputs: Single-band raster output from the model (road class)
  • ChesapeakeRSC input mosaic: the 4-band (red, green, blue, near-infrared) NAIP imagery the model runs on, shown here as RGB

View the interactive map here.

Selecting an Area of Interest (AOI)

We will choose an Area of Interest (AOI) for our analysis. ChesapeakeRSC was trained on imagery and land cover labels from Maryland.

We will try it out in Montgomery County, Maryland.

import os

import shapely
import wkls
from sedona.spark import *

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

# Generate a WKT geometry for Montgomery County, Maryland using Well-Known Locations (https://github.com/wherobots/wkls)
aoi = wkls.us.md.montgomerycounty.wkt()

min_lon, min_lat, max_lon, max_lat = shapely.from_wkt(aoi).bounds
print(f"Bounds: [{min_lon:.4f}, {min_lat:.4f}, {max_lon:.4f}, {max_lat:.4f}]")

Selecting a time range and verifying NAIP coverage

NAIP is flown state by state on a multi-year cycle at mixed resolutions, and the ChesapeakeRSC recipe sources 1m imagery. If you change the AOI or the time range, the 1m imagery the model needs may not exist. The USDA’s NAIP coverage map (PDF, 2002-2025) shows what is available where.

from datetime import datetime

# Date range for imagery to be used by the model
start_date = datetime(2017, 1, 1)
end_date = datetime(2018, 1, 1)

# The ChesapeakeRSC recipe runs on 1m NAIP imagery from the 4-band analytic product, since the model also takes
# near-infrared. Wherobots Open Data indexes every scene of that product with its footprint, resolution (gsd)
# and acquisition time.
MODEL_RES = 1.0
sedona.createDataFrame([(aoi,)], ["wkt"]).selectExpr("ST_GeomFromWKT(wkt) AS geom").createOrReplaceTempView("aoi")

# Scenes at the model's resolution touching the AOI (the table and wkls are both EPSG:4326)
sedona.sql(f"""
    SELECT n.geometry, n.year, n.datetime
    FROM wherobots_open_data.imagery.usda_naip_analytic_stac n
    JOIN aoi a ON ST_Intersects(n.geometry, a.geom)
    WHERE n.gsd = {MODEL_RES}
""").cache().createOrReplaceTempView("nearby")

# The recipe needs the AOI fully inside the union of the scene footprints in the date range
coverage = sedona.sql(f"""
    SELECT count(*) AS scenes, array_sort(collect_set(n.year)) AS years,
           ST_Contains(ST_Union_Agg(n.geometry), first(a.geom)) AS covered
    FROM nearby n CROSS JOIN aoi a
    WHERE n.datetime >= TIMESTAMP '{start_date.isoformat()}Z' AND n.datetime < TIMESTAMP '{end_date.isoformat()}Z'
""").first()

if not coverage["covered"]:
    # Years that would work, applying the same full-coverage test to each year on its own
    available = [
        row.year
        for row in sedona.sql("""
            SELECT n.year FROM nearby n CROSS JOIN aoi a
            GROUP BY n.year
            HAVING ST_Contains(ST_Union_Agg(n.geometry), first(a.geom))
            ORDER BY n.year
        """).collect()
    ]
    raise ValueError(
        f"{MODEL_RES:g}m NAIP does not fully cover this AOI between "
        f"{start_date:%Y-%m-%d} and {end_date:%Y-%m-%d}. "
        + (
            f"Years with complete {MODEL_RES:g}m coverage over this AOI: {available}."
            if available
            else f"No year has complete {MODEL_RES:g}m NAIP coverage over this AOI."
        )
    )

print(f"NAIP coverage is available. Found {coverage['scenes']} NAIP scenes at {MODEL_RES:g}m covering the AOI")
print(f"Acquisition years: {coverage['years']}")

Using the Wherobots Tasks API

Each RasterFlow workflow runs as a task through the Wherobots Python SDK. run_task below starts a task run, waits for it to finish, and returns its output value. Recipes return an object with input_mosaic and result; standalone tasks return a URI. Call wb.tasks.get("<task name>").input_schema to see every input a task accepts.

import time

from wherobots import TaskRunState, Wherobots

# In a Wherobots notebook the SDK uses your Wherobots credentials, elsewhere it reads WHEROBOTS_API_KEY
wb = Wherobots()


def run_task(task_name, inputs, poll_seconds=30):
    """Run a RasterFlow task, wait for it to finish, and return its output value."""
    run = wb.task_runs.get(wb.tasks.run(task_name, inputs, region="aws-us-west-2").task_run_id)
    print(f"Started {task_name} run {run.task_run_id}")
    while run.is_active:
        time.sleep(poll_seconds)
        run = wb.task_runs.get(run.task_run_id)
    if run.state != TaskRunState.SUCCEEDED:
        raise RuntimeError(f"{task_name} run {run.task_run_id} ended {run.state}: {run.error_message}")
    # A successful run can temporarily have no results if the output read fails.
    for _ in range(3):
        if run.results is not None:
            break
        time.sleep(poll_seconds)
        run = wb.task_runs.get(run.task_run_id)
    if run.results is None:
        raise RuntimeError(
            f"Results for successful run {run.task_run_id} are unavailable. "
            "Fetch this run again with wb.task_runs.get(run_id); do not resubmit the task."
        )
    # The single task output can be a URI or a structured recipe result.
    (output,) = run.results.values()
    return output

Running a model

RasterFlow has pre-defined workflows to simplify orchestration of the processing steps for model inference.These steps include:

  • Ingesting imagery for the specified Area of Interest (AOI)
  • Generating a seamless image from multiple image tiles (a mosaic)
  • Running inference with the selected model

The recipe returns input_mosaic for the source imagery and result for the predictions. Each contains a mosaic_index_uri for its GeoParquet index and a first_row_mosaic for its first Zarr store.

Note: This step will take approximately 30 minutes to complete.

recipe_output = run_task(
    "predict_mosaic_recipe",
    {
        # Area of Interest as WKT (EPSG:4326)
        "aoi": aoi,

        # Date range for imagery to be used by the model (set in the coverage check above), as RFC 3339 timestamps
        "start": start_date.isoformat() + "Z",
        "end": end_date.isoformat() + "Z",

        # Coordinate Reference System for the output mosaic
        "target_crs": "EPSG:3857",

        # The model recipe to be used for inference (ChesapeakeRSC in this case)
        "model_recipe": "chesapeakersc",
    },
)

# The recipe returns both the input imagery and the prediction mosaic.
input_mosaic_index_uri = recipe_output["input_mosaic"]["mosaic_index_uri"]
input_mosaic_store = recipe_output["input_mosaic"]["first_row_mosaic"]
print(f"Input imagery: {input_mosaic_store}")

model_output_index_uri = recipe_output["result"]["mosaic_index_uri"]

# Inspect the prediction index (one row per output mosaic location)
model_output_index_df = sedona.read.format("geoparquet").load(model_output_index_uri)
model_output_index_df.show(truncate=60)
# This example AOI has one geometry, so we select the first (only) mosaic location.
model_output_store = model_output_index_df.first()["location"]

model_output_store

(Optional) Build an optimized Zarr for visualization

RasterFlow writes its outputs as Zarr stores at native resolution. To explore them interactively on cloud.wherobots.com/map, you can build an optimized multiscale Zarr. build_zarr_multiscales adds downsampled overview levels (image pyramids) to the store so the map can stream coarse tiles when zoomed out and full-resolution pixels when zoomed in.

This step is optional and can take a few minutes for large outputs, so the code below is commented out by default — uncomment it to run it.

Pass aoi to select a chunk-aligned bounding window from a single Zarr source. It accepts polygon or multipolygon WKT in EPSG:4326, or an S3 GeoJSON or GeoParquet URI with a CRS. Pixels outside the polygon remain within that rectangular window. Omit aoi to process the full store.

# optimized_store = run_task("build_zarr_multiscales", {"source_store": model_output_store, "aoi": aoi})
# optimized_store

Visualizing outputs

You can visualize the Zarr, GeoParquet, and other geospatial outputs using cloud.wherobots.com/map.

Vectorize the raster model outputs

The output for the ChesapeakeRSC model is a raster with two classes: background and road.

We can run a separate flow to convert the roads into vector geometries, based on the confidence threshold.Converting these results to geometries allows us to more easily post process the results or join the results with other vector data.

import xarray as xr
import s3fs
import zarr
# Determine the classes that are predicted by the model
fs = s3fs.S3FileSystem(profile="default", asynchronous=True)
zstore = zarr.storage.FsspecStore(fs, path=model_output_store[5:])
ds = xr.open_zarr(zstore)
model_features = ds['band'].data.tolist()

# Remove the 'background' class for vectorization
vector_features = [f for f in model_features if f != 'background']  
# Note: this should take about 5 minutes to complete
vectorized_results_uri = run_task(
    "vectorize_mosaic",
    {
        "mosaic_uri": model_output_store,
        "features": vector_features,
        "threshold": 0.5,
        "stats": True,
        "medial_skeletonize": False,
    },
)

print(vectorized_results_uri)

Save the vectorized results to the catalog

We can store these vectorized outputs in the catalog by using WherobotsDB to persist the GeoParquet results.

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

df = sedona.read.format("geoparquet").load(vectorized_results_uri)
df = df.withColumnRenamed("label", "layer")
df.writeTo("examples_temp.chesapeake_db.chesapeake_vectorized").createOrReplace()

Visualize the vectorized results

To visualize the vectorized results, we will choose a small subset of the outputs that intersect with a bounding box around the center of the AOI.We will also filter out results with a score lower than 0.7.

from sedona.spark.sql.st_constructors import ST_MakeEnvelope
from sedona.spark.sql.st_predicates import ST_Intersects
import pyspark.sql.functions as f

# Compute a small bounding box around the center of the AOI (bounds from the AOI cell above)
center_lat = (min_lat + max_lat) / 2
center_lon = (min_lon + max_lon) / 2
delta = 0.01

min_lon = center_lon - delta
max_lon = center_lon + delta
min_lat = center_lat - delta
max_lat = center_lat + delta

df = df.filter(
    ST_Intersects(f.col("geometry"), ST_MakeEnvelope(min_lon, min_lat, max_lon, max_lat))
).filter("score_mean > 0.7")
from wherobots_gl import Map

# Wherobots-GL Map is URL-based, so write the filtered subset to GeoParquet first.
vectorized_subset_path = os.getenv("USER_S3_PATH") + "chesapeake_vectorized_subset.parquet"
df.write.format("geoparquet").mode("overwrite").save(vectorized_subset_path)

Map(
    layers=[{"type": "geoparquet", "source": vectorized_subset_path, "name": "Vectorized results"}],
    view={"zoom": 12, "lat": center_lat, "lng": center_lon},
    basemap="dark",
)

Generate PM Tiles for visualization

To improve visualization performance of a large number of geometries, we can use Wherobots built-in high performance PM tile generator.

from wherobots import vtiles

full_tiles_path = os.getenv("USER_S3_PATH") + "chesapeake_tiles.pmtiles"
vtiles.generate_pmtiles(df, full_tiles_path)
vtiles.show_pmtiles(full_tiles_path)
full_tiles_path

References

  1. Robinson, C., Corley, I., Ortiz, A., Dodhia, R., Lavista Ferres, J. M., & Najafirad, P. (2024). Seeing the roads through the trees: A benchmark for modeling spatial dependencies with aerial imagery. arXiv preprint arXiv:2401.06762. https://arxiv.org/abs/2401.06762
Related notebooks
VIEW MORE

RasterFlow SAM3

Detecting Objects From Text Prompts with RasterFlow This notebook will guide you through detecting objects...
LEARN MORE