Bring Your Own Rasters (BYOR) with RasterFlow

Private Preview

The following content is a read-only preview of an executable Jupyter notebook.To run this notebook interactively:

Go to Wherobots Cloud.
Start a runtime.
Open the notebook.
In the Jupyter Launcher:
1. Click File > Open Path.
2. Paste the following path to access this notebook: examples/Analyzing_Data/RasterFlow_Bring_Your_Own_Rasters_NAIP.ipynb
3. Click Enter.

This notebook demonstrates how to bring your own rasters (BYOR) into RasterFlow by querying a STAC catalog and creating a GTI (GDAL Raster Tile Index). You will learn how to query the NAIP collection, extract Cloud-Optimized GeoTIFF (COG) URLs, and build a mosaic using RasterFlow’s build_gti_mosaics function.

What you will learn

This notebook will teach you to:

Query a STAC catalog to discover available imagery for an area of interest
Create a GTI (GDAL Raster Tile Index) from STAC items
Build a seamless mosaic from multiple image tiles using RasterFlow’s build_gti_mosaics
Visualize the resulting mosaic and inference outputs
(Optional) Run road detection using the ChesapeakeRSC model

NAIP (National Agriculture Imagery Program)

The National Agriculture Imagery Program (NAIP) acquires aerial imagery during the agricultural growing seasons in the continental United States. Key characteristics:

Resolution: 1 meter ground sample distance
Bands: 4-band imagery (Red, Green, Blue, NIR) in a single COG file
Coverage: Continental United States
Update cycle: Typically every 2-3 years per state

We will use the Element84 Earth Search STAC catalog to query NAIP imagery for Montgomery County, Maryland.

Setup and Imports

import os

import geopandas as gpd
import wkls
from wherobots_gl import Map

Selecting an Area of Interest (AOI)

We will use Montgomery County, Maryland as our AOI. This county is well-covered by NAIP imagery and provides a good example for demonstrating the BYOR workflow.

# Generate a geometry for Montgomery County, MD using Well-Known Locations (https://github.com/wherobots/wkls)
gdf = gpd.read_file(wkls.us.md.montgomerycounty.geojson())

print(f"Bounding box: {gdf.total_bounds}")
print(f"CRS: {gdf.crs}")

# Display the AOI boundary on the Wherobots-GL map
minx, miny, maxx, maxy = gdf.total_bounds
center_lat = (miny + maxy) / 2
center_lon = (minx + maxx) / 2

# Wherobots-GL Map is URL-based, so write the AOI to GeoParquet first
aoi_path = os.getenv("USER_S3_PATH") + "montgomery_county_aoi.parquet"
gdf.to_parquet(aoi_path)

Map(
    layers=[{"type": "geoparquet", "source": aoi_path, "name": "Montgomery County, MD"}],
    view={"zoom": 10, "lat": center_lat, "lng": center_lon},
    basemap="satellite",
)

Loading the NAIP STAC collection

We use Wherobots’ built-in STAC reader to load the NAIP collection from Element84’s Earth Search catalog. This catalog provides direct S3 URLs to the imagery files.

from pyspark.sql import functions as F
from sedona.spark import *

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

# Get the actual polygon geometry as WKT for spatial filter
aoi_wkt = gdf.geometry.iloc[0].wkt

print(f"AOI geometry type: {gdf.geometry.iloc[0].geom_type}")

# Load NAIP collection from Element84 Earth Search
naip_df = (
    sedona.read.format("stac")
    .option("itemsLimitMax", "200")
    .option("itemsLimitPerRequest", "50")
    .load("https://earth-search.aws.element84.com/v1/collections/naip")
)

naip_df.createOrReplaceTempView("naip_raw")

print(f"\nLoaded NAIP STAC collection")
print(f"Schema:")
naip_df.printSchema()

Applying spatial and temporal filters

We filter the NAIP collection to only include tiles that:

Intersect with our AOI (Montgomery County)
Were captured in 2017 (the most recent year with 1-meter resolution NAIP imagery available for Maryland so we can test the ChesapeakeRSC model)

# Filter by:
# - Date range: 2017 (full year)
# - Spatial: Intersects with Montgomery County polygon boundary

filtered_df = sedona.sql(f"""
    SELECT 
        id,
        geometry,
        datetime,
        assets
    FROM naip_raw
    WHERE datetime BETWEEN '2017-01-01' AND '2017-12-31'
      AND ST_Intersects(geometry, ST_GeomFromText('{aoi_wkt}'))
""")

filtered_df.createOrReplaceTempView("naip_filtered")

tile_count = filtered_df.count()
print(f"Found {tile_count} NAIP tiles for Montgomery County, MD in 2017")
filtered_df.select("id", "datetime").show(10, truncate=False)

Exploring the STAC asset structure

Before creating the GTI, we need to understand the structure of the STAC assets. NAIP stores all 4 bands (RGBIR) in a single COG file, so we only need one row per tile in our GTI.

# Examine the assets structure
# NAIP has a single 'image' asset containing all 4 bands (RGBIR)

sample_df = sedona.sql("SELECT id, assets FROM naip_filtered LIMIT 1")
sample = sample_df.first()

if sample:
    print(f"Tile ID: {sample['id']}")
    print(f"\nAvailable asset keys:")
    for key in sorted(sample['assets'].keys()):
        print(f"  - {key}")
    
    # Check structure of the 'image' asset (4-band COG)
    print(f"\n--- Structure of 'image' asset ---")
    image_asset = sample['assets'].get('image')
    if image_asset:
        print(f"Type: {type(image_asset)}")
        print(f"Content: {image_asset}")
else:
    print("No samples found - check filters")

# Show the S3 URL for the 4-band COG

sample_df = sedona.sql("""
    SELECT 
        id,
        assets['image'].href as image_url
    FROM naip_filtered 
    LIMIT 5
""")

sample_df.show(truncate=False)

Creating the GTI (GDAL Raster Tile Index)

The GTI is a spatial index that maps tile geometries to their corresponding COG URLs. For more information on GTI, see the reference documentation. Since NAIP has a single 4-band COG per tile, we create one row per tile.

# Create the GTI
# NAIP has a single 4-band COG per tile, so we just need one row per tile

gti_df = sedona.sql("""
    SELECT 
        id as tile_id,
        ST_AsText(geometry) as geometry_wkt,
        geometry,
        datetime,
        assets['image'].href as url
    FROM naip_filtered
    WHERE assets['image'] IS NOT NULL
""")

# Filter out any rows where URL is null
gti_df = gti_df.filter(F.col("url").isNotNull())

gti_df.createOrReplaceTempView("gti")

row_count = gti_df.count()
print(f"GTI has {row_count} rows (one row per NAIP tile)")
gti_df.select("tile_id", "datetime", "url").show(10, truncate=80)

Saving the GTI as GeoParquet

We convert the Spark DataFrame to a GeoDataFrame and save it as GeoParquet for use with RasterFlow.

from shapely import wkt

# Select columns for the GTI (excluding the Sedona geometry, using WKT instead)
gti_export_df = gti_df.select(
    "tile_id",
    "geometry_wkt",
    "datetime",
    "url"
)

# Convert Spark DataFrame to Pandas
gti_pandas = gti_export_df.toPandas()

# Convert WKT geometry to Shapely geometry and create GeoDataFrame
gti_pandas['geometry'] = gti_pandas['geometry_wkt'].apply(wkt.loads)
gti_gdf = gpd.GeoDataFrame(gti_pandas, geometry='geometry', crs="EPSG:4326")

# Drop the WKT column (no longer needed)
gti_gdf = gti_gdf.drop(columns=['geometry_wkt'])

# Save to S3 as GeoParquet
gti_path = os.getenv("USER_S3_PATH") + "naip_gti.parquet"
gti_gdf.to_parquet(gti_path)

print(f"GTI saved to: {gti_path}")
print(f"\nGTI Schema:")
print(gti_gdf.dtypes)
print(f"\nSample rows:")
gti_gdf.head()

Verifying the GTI

Before building the mosaic, we verify that the GTI was saved correctly and visualize the tile footprints.

# Reload and verify the GTI
gti_verify = gpd.read_parquet(gti_path)

print(f"GTI contains {len(gti_verify)} rows")
print(f"Unique tiles: {gti_verify['tile_id'].nunique()}")
print(f"Date range: {gti_verify['datetime'].min()} to {gti_verify['datetime'].max()}")

print(f"\nSample URLs:")
for i, row in gti_verify.head(3).iterrows():
    print(f"  {row['tile_id']}: {row['url'][:70]}...")

# Visualize the GTI footprints with the AOI boundary on the Wherobots-GL map

# Calculate center from the AOI
minx, miny, maxx, maxy = gdf.total_bounds
center_lat = (miny + maxy) / 2
center_lon = (minx + maxx) / 2

print(f"Showing {len(gti_verify)} NAIP tile footprints over Montgomery County, MD")

Map(
    layers=[
        {"type": "geoparquet", "source": gti_path, "name": "NAIP Tile Footprints"},
        {"type": "geoparquet", "source": aoi_path, "name": "Montgomery County, MD"},
    ],
    view={"zoom": 10, "lat": center_lat, "lng": center_lon},
    basemap="satellite",
)

Building a mosaic with RasterFlow

Now we use RasterFlow’s build_gti_mosaics function to create a seamless 4-band mosaic from the NAIP tiles.

Note: The naip-analytic S3 bucket is requester-pays, so we set requester_pays=True.

Expected runtime: ~5-10 minutes depending on the number of tiles.

from rasterflow_remote import RasterflowClient

client = RasterflowClient()

# Build a 4-band mosaic (Red, Green, Blue, NIR)
mosaic_index = client.build_gti_mosaics(
    # Path to the GTI GeoParquet file
    gti=os.getenv("USER_S3_PATH") + "naip_gti.parquet",
    
    # Area of Interest (AOI)
    aoi=gdf,
    
    # Output bands (NAIP has 4 bands: r, g, b, ir)
    # Note: We use 'ir' to match the ChesapeakeRSC model's expected input bands
    bands=["r", "g", "b", "ir"],
    
    # Column in the GTI containing the COG URLs
    location_field="url",
    # Group tiles for time dimension by any column in the gti index
    time_column="datetime",
    
    # Output CRS (UTM Zone 18N for Maryland)
    target_crs=32618,
    
    # NAIP bucket is requester-pays
    requester_pays=True,
    
    # NAIP uses 0 as nodata value
    nodata=0.0,
)

# This example AOI has one geometry, so we select the first (only) mosaic location.
mosaic_store = mosaic_index.first_row_mosaic

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

(Optional) Build an optimized Zarr for visualization

RasterFlow writes its outputs as Zarr stores at native resolution. To explore them smoothly with the wherobots_gl.Map() widget, 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.

# optimized_store = client.build_zarr_multiscales(source_store=mosaic_store)
# optimized_store

Visualizing outputs

You can visualize the mosaic and inference outputs from this notebook inline with the wherobots_gl.Map() widget — it renders Zarr, GeoParquet, and COG directly in the notebook (as shown above for the AOI and tile footprints), so you can inspect either the mosaic output or the road-probability output without leaving the notebook.

(Optional) Running inference with ChesapeakeRSC

Now that we have built the NAIP mosaic, we can optionally run a semantic segmentation model to detect roads. The ChesapeakeRSC model was trained on NAIP imagery from Maryland and can detect roads even when occluded by tree cover. The model outputs two classes:

background
road

Expected runtime: ~10-20 minutes depending on the mosaic size.

# # Uncomment to run road detection on the NAIP mosaic
# 
# from rasterflow_remote import MosaicToMosaicActorEnum
# from rasterflow_remote.data_models import MergeModeEnum
# 
# # ChesapeakeRSC model - direct URL (no HuggingFace auth required)
# MODEL_PATH = "https://huggingface.co/wherobots/chesapeakersc-pt2/resolve/main/chesapeakersc-ep.pt2"
# 
# prediction_store = client.predict_mosaic(
#     mosaics=mosaic_store,
#     model_path=MODEL_PATH,
#     patch_size=512,
#     clip_size=64,
#     device="cuda",
#     features=["r", "g", "b", "ir"],  # Must match the bands in our mosaic
#     labels=["background", "road"],
#     actor=MosaicToMosaicActorEnum.SEMANTIC_SEGMENTATION_PYTORCH,
#     max_batch_size=256,
#     merge_mode=MergeModeEnum.WEIGHTED_AVERAGE,
# )
# 
# print(f"Prediction index URI: {prediction_store.uri}")

You can visualize the road-probability output from the inference workflow inline as well — point a wherobots_gl.Map() at the prediction Zarr store (e.g. prediction_store.uri), the same way we mapped the AOI and tile footprints above.

Getting Started

Spatial Queries

WherobotsDB

Data Connections

RasterFlow

Advanced Topics

GeoStats, Map Matching & RasterFlow

Bring Your Own Rasters (BYOR) with RasterFlow

What you will learn

NAIP (National Agriculture Imagery Program)

Setup and Imports

Selecting an Area of Interest (AOI)

Loading the NAIP STAC collection

Applying spatial and temporal filters

Exploring the STAC asset structure

Creating the GTI (GDAL Raster Tile Index)

Saving the GTI as GeoParquet

Verifying the GTI

Building a mosaic with RasterFlow

(Optional) Build an optimized Zarr for visualization

Visualizing outputs

(Optional) Running inference with ChesapeakeRSC

​What you will learn

​NAIP (National Agriculture Imagery Program)

​Setup and Imports

​Selecting an Area of Interest (AOI)

​Loading the NAIP STAC collection

​Applying spatial and temporal filters

​Exploring the STAC asset structure

​Creating the GTI (GDAL Raster Tile Index)

​Saving the GTI as GeoParquet

​Verifying the GTI

​Building a mosaic with RasterFlow

​(Optional) Build an optimized Zarr for visualization

​Visualizing outputs

​(Optional) Running inference with ChesapeakeRSC

What you will learn

NAIP (National Agriculture Imagery Program)

Setup and Imports

Selecting an Area of Interest (AOI)

Loading the NAIP STAC collection

Applying spatial and temporal filters

Exploring the STAC asset structure

Creating the GTI (GDAL Raster Tile Index)

Saving the GTI as GeoParquet

Verifying the GTI

Building a mosaic with RasterFlow

(Optional) Build an optimized Zarr for visualization

Visualizing outputs

(Optional) Running inference with ChesapeakeRSC