---
title: "RasterFlow Bring Your Own Rasters NAIP"
url: "https://wherobots.com/solution-gallery/rasterflow-bring-your-own-rasters-naip/"
published: "2026-01-23"
modified: "2026-11-02"
---

# RasterFlow Bring Your Own Rasters NAIP

## Bring Your Own Rasters (BYOR) with RasterFlow

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)](https://www.usgs.gov/centers/eros/science/usgs-eros-archive-aerial-photography-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 shapely
import wkls
from pyproj import CRS
import pyspark.sql.functions as f
from sedona.spark import *
from wherobots_gl import Map

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

## 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 WKT geometry for Montgomery County, MD 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"Bounding box: [{min_lon:.4f}, {min_lat:.4f}, {max_lon:.4f}, {max_lat:.4f}]")
center_lat = (min_lat + max_lat) / 2
center_lon = (min_lon + max_lon) / 2

# Wherobots-GL Map is URL-based, so write the AOI to GeoParquet for the map layer
aoi_map_path = os.getenv("USER_S3_PATH") + "montgomery_county_aoi.parquet"
aoi_df = sedona.createDataFrame([(aoi,)], ["wkt"]).select(ST_GeomFromWKT("wkt").alias("geometry"))
aoi_df.write.format("geoparquet").mode("overwrite").save(aoi_map_path)

Map(
    layers=[{"type": "geoparquet", "source": aoi_map_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.

```
print(f"AOI geometry type: {shapely.from_wkt(aoi).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}'))
""")

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](https://gdal.org/en/stable/drivers/raster/gti.html).

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,
        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 write the WherobotsDB DataFrame to GeoParquet for use with RasterFlow, tagging the footprints with their EPSG:4326 CRS.

```
# Select the GTI columns and write them to S3 as GeoParquet
gti_export_df = gti_df.select("tile_id", "geometry", "datetime", "url")

gti_path = os.getenv("USER_S3_PATH") + "naip_gti.parquet"
(
    # Tag the footprints as EPSG:4326; RasterFlow reads the whole directory Spark writes as one dataset
    gti_export_df.write.format("geoparquet")
    .option("geoparquet.crs", CRS.from_epsg(4326).to_json())
    .mode("overwrite")
    .save(gti_path)
)

print(f"GTI saved to: {gti_path}")
print("nGTI Schema:")
gti_export_df.printSchema()
print("nSample rows:")
gti_export_df.show(3, truncate=60)
```

## 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 = sedona.read.format("geoparquet").load(gti_path)

print(f"GTI contains {gti_verify.count()} rows")
print(f"Unique tiles: {gti_verify.select('tile_id').distinct().count()}")
earliest, latest = gti_verify.agg(f.min("datetime"), f.max("datetime")).first()
print(f"Date range: {earliest} to {latest}")

print("nSample rows:")
gti_verify.select("tile_id", "datetime", "url").show(3, truncate=70)
```

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

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

Map(
    layers=[
        {"type": "geoparquet", "source": gti_path, "name": "NAIP Tile Footprints"},
        {"type": "geoparquet", "source": aoi_map_path, "name": "Montgomery County, MD"},
    ],
    # Centered on the AOI, using the bounds computed in the AOI cell above
    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.

```
import time

import shapely

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 URI."""
    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."
        )
    # Each RasterFlow task has one output, named by the task (e.g. mosaic_index_uri)
    (output_uri,) = run.results.values()
    return output_uri
```

```
# Build a 4-band mosaic (Red, Green, Blue, NIR)
mosaic_index_uri = run_task(
    "build_gti_mosaics",
    {
        # Path to the GTI GeoParquet file
        "gti_uri": gti_path,

        # Area of interest as WKT in EPSG:4326
        "aoi": aoi,

        # 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": "EPSG: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 = sedona.read.format("geoparquet").load(mosaic_index_uri).first()["location"]

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.

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": mosaic_store, "aoi": aoi})
# 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](https://github.com/isaaccorley/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
#
# prediction_index_uri = run_task(
#     "run_mosaic_inference",
#     {
#         "mosaic_uri": mosaic_store,
#         # ChesapeakeRSC model - direct URL (no HuggingFace auth required)
#         "model_path": "https://huggingface.co/wherobots/chesapeakersc-pt2/resolve/b1729ac4be74040f66d0caf5c0f83d0477d9e091/chesapeakersc-ep.pt2",
#         "patch_size": 512,
#         "clip_size": 64,
#         "device": "cuda",
#         "features": ["r", "g", "b", "ir"],  # Must match the bands in our mosaic
#         "labels": ["background", "road"],
#         "actor": "semantic_segmentation_pytorch",
#         "max_batch_size": 256,
#         "merge_mode": "weighted_average",
#     },
# )
#
# print(f"Prediction index URI: {prediction_index_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 (the `location` column of the prediction index), the same way we mapped the AOI and tile footprints above.