Skip to main content
Private Preview
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.
  2. Select the specific notebook you wish to run.
  3. Click Run Model in Notebook.
This notebook demonstrates how to generate a harvest season Sentinel-2 median mosaic using Wherobots RasterFlow. You will learn how to define an Area of Interest, build a seasonal composite mosaic, and visualize the results as RGB imagery.

​
Sentinel-2 Harvest and Growing Season Mosaics

RasterFlow provides pre-built seasonal mosaic datasets that use latitude-based heuristics to determine the optimal imaging window for agricultural analysis. S2_MED_PLANTING: Sentinel-2 median composite dataset for planting season. S2_MED_HARVEST: Sentinel-2 median composite dataset for harvest season. In this example, we will generate a mosaic for the harvest season. The mosaic contains 5 bands:
  • s2med_harvest:B02 - Blue (10m)
  • s2med_harvest:B03 - Green (10m)
  • s2med_harvest:B04 - Red (10m)
  • s2med_harvest:B08 - NIR (10m)
  • s2med_harvest:N_VALID_PIXELS - Count of valid pixels used in the median calculation

​
Setup and Imports

import os
from datetime import datetime

import geopandas as gpd
import hvplot.xarray
import s3fs
import wkls
import xarray as xr
import zarr
from pyproj import Transformer
from rasterflow_remote import RasterflowClient
from rasterflow_remote.data_models import DatasetEnum

​
Initializing the RasterFlow client

rf_client = RasterflowClient()

​
Selecting an Area of Interest (AOI)

We will use Haskell County, Kansas as our AOI. This is an agricultural region in the US Great Plains, ideal for demonstrating harvest-season imagery. 
# Generate a geometry for Haskell County, Kansas using Well-Known Locations (https://github.com/wherobots/wkls)
gdf = gpd.read_file(wkls["us"]["ks"]["Haskell County"].geojson())

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

# Store bounds for visualization later
min_lon, min_lat, max_lon, max_lat = gdf.total_bounds

print(f"AOI saved to: {aoi_path}")
print(f"Bounds: [{min_lon:.4f}, {min_lat:.4f}, {max_lon:.4f}, {max_lat:.4f}]")

​
Building the harvest season mosaic

This step builds a Sentinel-2 median composite mosaic for the harvest season only. The workflow will:
  1. Query Sentinel-2 imagery for the AOI and date range
  2. Filter to harvest season dates based on latitude
  3. Apply cloud filtering (< 75% cloud cover)
  4. Compute a pixel-wise median composite
  5. Output the result as a Zarr store
Note: This step will take approximately 10-15 minutes to complete. If you want to skip this build_mosaic step, you can uncomment the next cell and start with a pre-generated mosaic. 
# Use the pre-generated mosaic instead
# mosaic_output = "s3://wherobots-examples/rasterflow/mosaics/haskell_harvest.zarr"

mosaic_output = rf_client.build_mosaic(
    # Dataset type for harvest season median mosaic
    datasets=[DatasetEnum.S2_MED_HARVEST],
    
    # Path to our AOI in GeoParquet format
    aoi=aoi_path,
    
    # Date range for imagery to be used
    start=datetime(2024, 1, 1),
    end=datetime(2025, 1, 1),
    
    # Output CRS (Web Mercator)
    crs_epsg=3857,
)

print(f"Mosaic saved to: {mosaic_output}")

​
Visualizing the mosaic

We load the Zarr store and visualize the harvest mosaic as an RGB composite. 
# Open the Zarr store from S3
fs = s3fs.S3FileSystem(profile="default", asynchronous=True)
zstore = zarr.storage.FsspecStore(fs, path=mosaic_output.removeprefix('s3://'))  # FsspecStore expects a path without the 's3://' scheme prefix
ds = xr.open_zarr(zstore)

# Create a transformer to convert from lat/lon to meters
transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)

# Transform bounding box coordinates from lat/lon to meters
(min_x, max_x), (min_y, max_y) = transformer.transform(
    [min_lon, max_lon],
    [min_lat, max_lat]
)

# Subset the data to our AOI bounds
# Determine y-coordinate order (ascending vs descending) to slice correctly
if ds.y[0] < ds.y[-1]:
    ds_subset = ds.sel(x=slice(min_x, max_x), y=slice(min_y, max_y))
else:
    ds_subset = ds.sel(x=slice(min_x, max_x), y=slice(max_y, min_y))

# Select RGB bands and access the variables array
# Band names are prefixed with the dataset name (e.g., s2med_harvest:B04)
rgb = ds_subset.sel(band=["s2med_harvest:B04", "s2med_harvest:B03", "s2med_harvest:B02"])["variables"]

# Handle time dimension if present
if "time" in rgb.dims:
    rgb = rgb.isel(time=0)

# Normalize values for display (Sentinel-2 reflectance * 10000)
rgb_normalized = rgb / 3000
rgb_normalized = rgb_normalized.clip(0, 1)

# Visualize as RGB composite
rgb_normalized.hvplot.rgb(
    x="x",
    y="y",
    bands="band",
    data_aspect=1,
    xaxis=False,
    yaxis=False,
    title="Harvest Season Mosaic - Haskell County, KS (2024)",
    frame_width=600,
)

​
Summary

This notebook demonstrated how to:
  1. Use RasterFlow to build a harvest-season Sentinel-2 median mosaic
  2. Load and visualize the resulting Zarr store as an RGB composite

​
Next steps

  • Compare with a planting season mosaic using DatasetEnum.S2_MED_PLANTING
  • Run the FTW model for field boundary detection (see RasterFlow_FTW.ipynb)
  • Use the mosaic as input for model inference using predict_mosaic