Change Detection with RasterFlow

Public 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_ChangeDetection.ipynb
3. Click Enter.

This notebook demonstrates RasterFlow’s change detection inference workflow using the Fields of the World (FTW)¹ model as an example. Related notebooks:

RasterFlow_FTW.ipynb — Complete FTW field boundary detection tutorial
RasterFlow_Bring_Your_Own_Model.ipynb — Export custom PyTorch models

Why Fields of the World (FTW)?

FTW predicts 3 classes (non_field_background, field, field_boundaries). While FTW is a field boundary segmentation model (not a true change detection model), it takes bi-temporal input (two seasons), making it useful for demonstrating change detection.

Model signature

Change detection models receive multi-temporal imagery stacked along the channel dimension:

	Shape
Input	`(Batch, Channels × TimeSteps, H, W)` — e.g. `(N, 8, 256, 256)` for 4 bands × 2 seasons
Output	`(Batch, Classes, H, W)` — e.g. `(N, 3, 256, 256)` for 3 classes

class ChangeDetectionModel(nn.Module):
    def forward(self, x):  # x: (N, Channels*TimeSteps, H, W)
        return logits      # (N, Classes, H, W)

We’ll use a particular Actor to run change detection inference using Rasterflow: semantic_segmentation_change_detection_pytorch.

Selecting an Area of Interest (AOI)

To start, we will choose an Area of Interest (AOI) for our analysis. The area around Haskell County, Kansas has some interesting crop field patterns so we will try out the model there.

import os

import shapely
import wkls
from sedona.spark import *

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

# Generate a WKT geometry for Haskell County, Kansas using Well-Known Locations (https://github.com/wherobots/wkls)
aoi = wkls.us.ks.haskellcounty.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}]")

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 the URI of its output. Call wb.tasks.get("<task name>").input_schema to see every input a task accepts.

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

Building two annual median composites

Use s2_median with median_compositing_interval="annual" to create the two input time steps needed by this model. The UTC request [2023-01-01, 2025-01-01) produces medians for [2023-01-01, 2024-01-01) and [2024-01-01, 2025-01-01), labelled 2023-01-01 and 2024-01-01. Without explicit annual grouping, s2_median produces one median over the entire request window. That output would have only one time step and would not meet this example’s two-time-step input contract. RasterFlow applies quality filtering and cloud masking before computing per-pixel medians and writing the Zarr store. Reflectance bands use the s2_median prefix.

mosaic_index_uri = run_task(
    "build_mosaics",
    {
        "datasets": ["s2_median"],
        # Area of interest as WKT in EPSG:4326
        "aoi": aoi,
        "start": "2023-01-01T00:00:00Z",
        "end": "2025-01-01T00:00:00Z",
        "median_compositing_interval": "annual",
        # UTM zone 14N, which covers Haskell County
        "target_crs": "EPSG:32614",
    },
)

# Inspect the mosaic index (one row per output mosaic location)
mosaic_index_df = sedona.read.format("geoparquet").load(mosaic_index_uri)
mosaic_index_df.show(truncate=60)

# This example AOI has one geometry, so we select the first (only) mosaic location.
input_store = mosaic_index_df.first()["location"]

input_store

Inspect the Zarr store with xarray. The annual grouping creates two time steps, labelled 2023-01-01 and 2024-01-01, for this change-detection example.

import xarray as xr
import s3fs
import zarr

fs = s3fs.S3FileSystem(profile="default", asynchronous=True)
zstore = zarr.storage.FsspecStore(fs, path=input_store[5:])
mosaic = xr.open_zarr(zstore)

mosaic['variables']

We can also see what bands are included in the mosaic.

mosaic['variables']['band']

Mosaic Inference

Now that we have built the mosaic and understand the features that were included, we can run the run_mosaic_inference task. We’ll pass in the band names we want to predict, as well as other inference configs that specify how to run inference.

mosaic_uri is the mosaic we just created for sentinel-2
model_path can be a path to a Pytorch 2 Archive file on s3 or Huggingface. In this case we’ll use the Fields of the World model from the Wherobots’ Huggingface Collection.
patch_size controls the XY size of the array input to the model
clip_size in conjunction with merge_mode, controls how to run overlapping windowed inference to reduce edge effects relative to non-overlapping inference
device specifies what device to use for running inference. Our runtimes are GPU only at this time, but you could run on the CPU device if you wanted to!
features selects what bands and in what order to run inference on
labels sets the labels on the models output. Order matters here based on what your model returns!
actor determines the kind of inference handler to use. The correct actor to select depends on the task your model is solving and the kind of input it expects. For change detection, we support models that take in an input where time and channels dimensions are stacked together. So the channel dimension is ordered like so in in this example:

        "<time_1>s2_median:B04",
        "<time_1>s2_median:B03",
        "<time_1>s2_median:B02",
        "<time_1>s2_median:B08",
        "<time_2>s2_median:B04",
        "<time_2>s2_median:B03",
        "<time_2>s2_median:B02",
        "<time_2>s2_median:B08"

max_batch_size: Maximum number of patches fed through the model in a single forward pass.
- Higher values = better GPU utilization but more GPU memory required
- Lower values = less memory but potentially slower inference

model_output_index_uri = run_task(
    "run_mosaic_inference",
    {
        "mosaic_uri": input_store,
        "model_path": "https://huggingface.co/wherobots/ftw-v1.1-pt2/resolve/b20b5238c72e2db0954bd836523efc588d259940/ftw-v1.1-ep.pt2",
        "patch_size": 256,
        "clip_size": 32,
        "device": "cuda",
        "features": [
            "s2_median:B04",
            "s2_median:B03",
            "s2_median:B02",
            "s2_median:B08",
        ],
        "labels": [
            "non_field_background",
            "field",
            "field_boundaries",
        ],
        "actor": "semantic_segmentation_change_detection_pytorch",
        "max_batch_size": 128,
        "merge_mode": "weighted_average",
    },
)

# Inspect the mosaic 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

Let’s check out our mosaic result with xarray to see what a Zarr change detection store looks like.

import s3fs
import zarr

fs = s3fs.S3FileSystem(profile="default", asynchronous=True)
zstore = zarr.storage.FsspecStore(fs, path=model_output_store[5:])
result_mosaic = xr.open_zarr(zstore)

The prediction result mosaic is shaped similar to the input, height and width are the same. Rather than have 4 bands, it has 3 bands for the prediction result, representing the categories. The FTW model outputs 3 categories, whereas other models may output a binary classification and there will be only one prediction band. Instead of 2 time steps, there is only 1, representing the change interval. We index the time coordinate of the Zarr store using the start date used to build the input mosaic. We also store a time delta coordinate to represent the interval of time between the start_date and end_date.

result_mosaic['variables']

start_time = result_mosaic['variables'].coords['time'].values[0]
print("start time:", start_time)
end_time = (result_mosaic['variables'].coords['time'] + result_mosaic['variables'].coords['time_delta']).values[0]
print("end time", end_time)

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 FTW model is a raster with three classes as bands: field, field_boundaries, and non_field_background. We will run a separate flow to convert the fields and field boundaries into vector geometries. Converting these results to geometries allows us to more easily post process the results or join the results with other vector data.

# Determine the classes that are predicted by the model
model_features = result_mosaic['band'].data.tolist()

# Remove the 'non_field_background' class for vectorization
vector_features = [f for f in model_features if f != 'non_field_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.ftw_db")

vectorize_result_df = sedona.read.format("geoparquet").load(vectorized_results_uri)
vectorize_result_df = vectorize_result_df.withColumnRenamed("label", "layer")
vectorize_result_df.writeTo("examples_temp.ftw_db.ftw_vectorized").createOrReplace()

Visualize the vectorized results

To visualize the vectorized results, we will show the fields around Plymell, Kansas and filter out results with a score lower than 0.5. 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.

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

# Bounding box for Plymell, Kansas
plymell_wkt = "POLYGON ((-100.98 37.83, -100.98 37.68, -100.8 37.68, -100.8 37.83, -100.98 37.83))"
vectorize_result_filtered_df = vectorize_result_df.filter("layer == 'field'").filter("score_mean > 0.5") \
    .filter(ST_Intersects(f.col("geometry"), ST_GeomFromText(lit(plymell_wkt))))

We will apply some WherobotsDB vector post processing operations to refine the model results.

tolerance = 0.00001
decimal_places = 6  # roughly 10cm in lat/lon
df = vectorize_result_filtered_df.withColumn("geometry", f.expr("ST_MakeValid(geometry)"))
df = df.withColumn("geometry", f.expr("ST_SetSRID(geometry, 4326)")) # vectorize_mosaic defaults to returning geometries in the CRS EPSG:4326
df = df.withColumn("geometry", f.expr(f"ST_ReducePrecision(geometry, {decimal_places})"))
df = df.withColumn("geometry", f.expr(f"ST_SimplifyPreserveTopology(geometry, {tolerance})"))
df = df.repartition(200)
df.show()

from wherobots_gl import Map

vectorized_map_path = os.getenv("USER_S3_PATH") + "vectorized_results.parquet"
df.write.format("geoparquet").mode("overwrite").save(vectorized_map_path)
Map(layers=[{"type": "geoparquet", "source": vectorized_map_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. The FTW model has a tendency to create extremely large boundary geometries, which doesn’t play nicely with PMTiles. To avoid this we subdivide the boundary geometries.

fields_df = df.where("layer = 'field'")
# 1266 is the vertex count of the second-largest boundary geometry.
# We use this instead of the absolute largest boundary (an extreme outlier)
# so the subdivision threshold is high enough to keep most geometries intact
# while still forcing that outlier-sized boundary to be subdivided for
# better PMTiles rendering performance.
boundaries_df = df.where("layer = 'field_boundaries'").withColumn("geometry", f.expr("ST_SubDivideExplode(geometry, 1266)"))
tile_features_df = fields_df.withColumn("layer", f.lit("fields")).unionByName(
    boundaries_df.withColumn("layer", f.lit("boundaries"))
)

from wherobots import vtiles
full_tiles_path = os.getenv("USER_S3_PATH") + "tiles.pmtiles"
vtiles.generate_pmtiles(tile_features_df, full_tiles_path)

vtiles.show_pmtiles(full_tiles_path)

Viewing the PMTiles on the Wherobots Cloud map

Open the PMTiles archive on cloud.wherobots.com/map with the link printed below. The map reads the archive from your storage with your Wherobots credentials, so you don’t need a pre-signed URL.

from urllib.parse import quote

print(f"View on the map: https://cloud.wherobots.com/map?l0.src={quote(full_tiles_path, safe='')}")

References

Kerner, H., Chaudhari, S., Ghosh, A., Robinson, C., Ahmad, A., Choi, E., Jacobs, N., Holmes, C., Mohr, M., et al. (2024). Fields of The World: A Machine Learning Benchmark Dataset For Global Agricultural Field Boundary Segmentation. arXiv preprint arXiv:2409.16252. Accepted at AAAI-2025 Artificial Intelligence for Social Impact (AISI) track.
ESA. (2015). Sentinel-2 User Handbook (Issue 1, Rev. 2). European Space Agency. https://sentinels.copernicus.eu/documents/247904/685211/Sentinel-2_User_Handbook

​Why Fields of the World (FTW)?

​Model signature

​Selecting an Area of Interest (AOI)

​Using the Wherobots Tasks API

​Building two annual median composites

​Mosaic Inference

​(Optional) Build an optimized Zarr for visualization

​Visualizing outputs

​Vectorize the raster model outputs

​Save the vectorized results to the catalog

​Visualize the vectorized results

​Generate PM Tiles for visualization

​Viewing the PMTiles on the Wherobots Cloud map

​References