> ## Documentation Index
> Fetch the complete documentation index at: https://docs.wherobots.com/llms.txt
> Use this file to discover all available pages before exploring further.

# Detecting sidewalks with RasterFlow

<Badge color="purple">Private Preview</Badge>

<Tip>
  The following content is a read-only preview of an executable Jupyter notebook.

  To run this notebook interactively:

  1. Go to the [**Wherobots Model Hub**](https://cloud.wherobots.com/model-hub).
  2. Select the specific notebook you wish to run.
  3. Click **Run Model in Notebook**.
</Tip>

<img src="https://mintcdn.com/wherobots/XTZef6cJhZI1PP0b/tutorials/example-notebooks/images/rasterflow-tile2net-tile2net-banner.jpg?fit=max&auto=format&n=XTZef6cJhZI1PP0b&q=85&s=2576b29de06508a2595e13342be055f4" alt="tile2net banner" width="1292" height="231" data-path="tutorials/example-notebooks/images/rasterflow-tile2net-tile2net-banner.jpg" />

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

### Tile2Net

The [Tile2Net](https://github.com/VIDA-NYU/tile2net)<sup>1</sup> model is an open source segmentation model that can detect sidewalks and other pathways from very high resolution imagery.

This model predicts 4 classes:

* background
* road
* sidewalk
* crosswalk

It was trained on \~19cm aerial imagery from USGS for Cambridge, Manhattan, Brooklyn, and Washington DC.  We will demonstrate results using 30cm resolution data from the [National Agriculture Imagery Program (NAIP)](https://www.usgs.gov/centers/eros/science/usgs-eros-archive-aerial-photography-national-agriculture-imagery-program-naip).

## 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:**

* *Tile2Net input mosaic*: RGB high-resolution aerial imagery the model runs on
* *Tile2Net model outputs*: raster output from the model, with bands for sidewalk / road / crosswalk classes
* *Tile2Net vector geometries*: vectorized sidewalk / road / crosswalk polygons derived from the raster output
* *Tile2Net PM Tiles*: same vectors delivered as PMTiles for fast rendering at scale

View the interactive map [here](https://viz.wherobots.com/?z=14.24\&lat=38.99775\&lng=-76.93743\&l0.id=4b363cf8-20cf-464f-bbaf-820ae32d85de\&l0.src=s3%3A%2F%2Fwherobots-examples%2Frasterflow%2Fmosaics%2Fcollegepark_optimized.zarr\&l0.name=Tile2Net+input+mosaic\&l0.clim=24.672528253600106%2C200.6443167352033\&l0.mode=rgb\&l0.rgb=r%2Cg%2Cb%2Cband\&l0.var=variables\&l1.id=11862f01-e388-4788-bcec-0ade21eec794\&l1.src=s3%3A%2F%2Fwherobots-examples%2Frasterflow%2Fmodel-outputs%2Ftile2net_optimized.zarr\&l1.name=Tile2Net+model+outputs\&l1.clim=0.00009313203892032774%2C0.8193664079376207\&l1.mode=single\&l1.rgb=sidewalk%2Croad%2Ccrosswalk%2Cband\&l1.var=variables\&l2.id=78ecd25a-79c8-49c7-996e-6774a45cc05a\&l2.src=s3%3A%2F%2Fwherobots-examples%2Frasterflow%2Fvectors%2Ftile2net_collegepark.parquet%2F\&l2.name=Tile2Net+vector+geometries\&l2.type=geoparquet\&l2.fill=%233b82f6\&l2.line=%232563eb\&l2.lw=2\&l2.cc=%233b82f6\&l2.cr=5\&l3.id=7ba01f8a-40b6-41df-9766-5553491ec152\&l3.src=s3%3A%2F%2Fwherobots-examples%2Frasterflow%2Fmodel-outputs%2Ftile2net.pmtiles\&l3.name=Tile2Net+PM+Tiles\&l3.type=pmtiles\&l3.fill=%23ef4444\&l3.line=%23dc2626\&l3.lw=2\&l3.cc=%23ef4444\&l3.cr=5\&title=Tile2Net+College+Park+MD).

## Selecting an Area of Interest (AOI)

To start, we will choose an Area of Interest (AOI) for our analysis. Tile2Net was trained on select geographies in the Northeastern United States. So we will pick a location in that geographic region that wasn't in the training data and where these is recent 30cm resolution NAIP data: College Park, Maryland.

```python theme={"system"}
import wkls
import geopandas as gpd
import os

# Generate a geometry for College Park, Maryland using Well-Known Locations (https://github.com/wherobots/wkls)
# NOTE: later in the notebook we fix an area to inspect assuming you ran on College Park. If you change the location here
# you will need to update the visualization area in "Visualize a subset of the model outputs"
gdf = gpd.read_file(wkls.us.md.collegepark.geojson())

# Save the geometry to a parquet file in the user's S3 path
aoi_path = os.getenv("USER_S3_PATH") + "collegepark.parquet"
gdf.to_parquet(aoi_path)
```

## Selecting a time range and verifying NAIP coverage

Because NAIP collects state-level imagery at mixed resolutions, Tile2Net's required 30cm data may not exist for your selected AOI or timeframe. The USDA's [NAIP coverage map](https://www.arcgis.com/home/item.html?id=cdaa8c24cf0844abba74fb9d71432fd4#overview) (PDF, 2002-2025) shows what is available where.

```python theme={"system"}
from datetime import datetime

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

# The Tile2Net recipe runs on 30cm NAIP imagery. The index holds one row per NAIP scene, with its
# footprint, resolution (res) and acquisition time.
MODEL_RES = 0.3
naip_index = gpd.read_parquet(
    "s3://wherobots-examples/rasterflow/indexes/naip_index.parquet",
    columns=["geometry", "res", "year", "time"],
)

# Scenes at the model's resolution touching the AOI, compared in the index's CRS
aoi_geom = gdf.geometry.to_crs(naip_index.crs).union_all()
at_res = naip_index.query(f"res == {MODEL_RES}")
nearby = at_res.iloc[at_res.sindex.query(aoi_geom, predicate="intersects")]

# Narrowed to the requested date range
covering = nearby.query(f"time >= '{start_date:%Y-%m-%d}' and time <= '{end_date:%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_date:%Y-%m-%d} and {end_date:%Y-%m-%d}; the recipe needs the AOI fully covered. "
        + (
            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 {len(covering)} NAIP scenes at {MODEL_RES:g}m covering the AOI")
print(f"Acquisition years: {sorted(covering['year'].unique().tolist())}")
```

## Initializing the RasterFlow client

```python theme={"system"}
from rasterflow_remote import RasterflowClient

from rasterflow_remote.data_models import (
    ModelRecipes, 
    VectorizeMethodEnum
)

rf_client = RasterflowClient()
```

## 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 output is a Zarr store of the model outputs.

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

```python theme={"system"}
model_output_index = rf_client.predict_mosaic_recipe(
    # Path to our AOI in GeoParquet or GeoJSON format
    aoi = aoi_path,

    # Date range for imagery to be used by the model (set in the coverage check above)
    start = start_date,
    end = end_date,

    # Coordinate Reference System EPSG code for the output mosaic   
    target_crs = 3857,

    # The model recipe to be used for inference (Tile2Net in this case)
    model_recipe = ModelRecipes.TILE_2_NET,
)

model_output_index.mosaic_index_gdf
```

```python theme={"system"}
# This example AOI has one geometry, so we select the first (only) mosaic location.
model_output_store = model_output_index.first_row_mosaic

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](https://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.

```python theme={"system"}
# optimized_store = rf_client.build_zarr_multiscales(source_store=model_output_store)
# optimized_store
```

## Visualizing outputs

If RasterFlow is enabled for your organization, you can visualize the Zarr, GeoParquet, and other geospatial outputs using [cloud.wherobots.com/map](https://cloud.wherobots.com/map).

## Vectorize the raster model outputs

The output for the Tile2Net model is a raster with four classes: background, road, sidewalk, crosswalk.

We can run a seperate flow to convert the roads, sidewalks and crosswalks 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 resuilts with other vector data.

```python theme={"system"}
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()

# Only vectorize the 'sidewalk' and 'crosswalk' classes 
relevant_features = ['sidewalk', 'crosswalk']   
vector_features = [f for f in model_features if f in relevant_features]
```

```python theme={"system"}
# Note: this should take about 5 minutes to complete
vectorized_results = rf_client.vectorize_mosaic(
        mosaic = model_output_store,
        features = vector_features,
        threshold = 0.05,
        vectorize_method = VectorizeMethodEnum.SEMANTIC_SEGMENTATION_RASTERIO,
        vectorize_config={"stats": True, "medial_skeletonize": False}
    )

print(vectorized_results)
```

## Save the vectorized results to the catalog

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

```python theme={"system"}
from sedona.spark import *
import pyspark.sql.functions as f
from pyspark.sql.functions import expr

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

```python theme={"system"}
sedona.sql("CREATE DATABASE IF NOT EXISTS examples_temp.tile2net_db")

df = sedona.read.format("geoparquet").load(vectorized_results.uri)
df = df.withColumnRenamed("label", "layer")
df.writeTo("examples_temp.tile2net_db.tile2net_vectorized").createOrReplace()
```

## Visualize the vectorized results

To visualize the vectorized results, we will filter out results with a score lower than 0.2. This threshold was determined through observation to strike a balance: it eliminates obvious noise without being overly aggressive, ensuring that we don't accidentally filter out too many relevant results.

Gaps in prediction for Tile2net are likely due to it being used with lower resolution imagery (30cm resolution) than what it was trained on (19cm resolution), as well as changes in geographic context, and occlusion from trees over the pathways.

```python theme={"system"}
df = df.filter("score_mean > 0.2")
df_filtered = df.withColumn(
    "area_m2",
    expr("ST_AreaSpheroid(geometry)")
).filter("area_m2 > 1000")
```

```python theme={"system"}
from wherobots_gl import Map

# Wherobots-GL Map is URL-based, so write the derived (filtered) vectors to GeoParquet first.
vectors_viz_path = os.getenv("USER_S3_PATH") + "tile2net_vectorized.parquet"
df_filtered.write.format("geoparquet").mode("overwrite").save(vectors_viz_path)

Map(layers=[{"type": "geoparquet", "source": vectors_viz_path, "name": "Vectorized results"}])
```

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

```python theme={"system"}
from wherobots import vtiles

full_tiles_path = os.getenv("USER_S3_PATH") + "tiles.pmtiles"
vtiles.generate_pmtiles(df_filtered, full_tiles_path)
```

```python theme={"system"}
vtiles.show_pmtiles(full_tiles_path)
```

### References

1. **Hosseini, M., Sevtsuk, A., Miranda, F., Cesar Jr, R. M., & Silva, C. T. (2023).** Mapping the walk: A scalable computer vision approach for generating sidewalk network datasets from aerial imagery. *Computers, Environment and Urban Systems*, *101*, 101950.