Skip to content

image image image

TerraScope STAC API with leafmap

This notebook demonstrates how to access TerraScope data using the leafmap.terrascope module.

Setup: Set environment variables before running: 

1
2
export TERRASCOPE_USERNAME='your_username'
export TERRASCOPE_PASSWORD='your_password'

Installation

Uncomment the following line to install leafmap if needed.

1
# %pip install -U leafmap

Import libraries

1
2
import leafmap
import leafmap.terrascope as terrascope

Authentication

1
terrascope.login()

Search for NDVI Data

Search data from the terrascope-s2-ndvi-v2 collection.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
# Brussels area: [west, south, east, north]
bbox = [4.2194, 50.7564, 4.4790, 50.9207]
items = terrascope.search_ndvi(
    bbox=bbox,
    start="2025-05-01",
    end="2025-06-01",
    max_cloud_cover=10,
)

print(f"Found {len(items)} scenes:")
for date in terrascope.get_item_dates(items):
    print(f"  {date}")

Single Layer Visualization

Visualize a single NDVI scene. Raw pixel values (0-250) are scaled to NDVI using value * 0.004 - 0.08, giving a range of [-0.08, 0.92].

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
center = [(bbox[1] + bbox[3]) / 2, (bbox[0] + bbox[2]) / 2]

m = leafmap.Map(center=center, zoom=14)
m.add_basemap("Esri.WorldImagery")
m.add_raster(
    items[0].assets["NDVI"].href,
    layer_name=f"NDVI {items[0].datetime.date()}",
    colormap="RdYlGn",
    vmin=0,
    vmax=250,  # scale = 0.004, offset = -0.08
)
m.add_colormap(cmap="RdYlGn", label="NDVI", vmin=-0.08, vmax=0.92)
m

Time Slider Animation

1
2
3
4
5
layers = terrascope.create_time_layers(items)

m = leafmap.Map(center=center, zoom=14)
m.add_time_slider(layers, time_interval=1)
m

Data Analysis with rioxarray

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
import rioxarray
import numpy as np

# Clean up stale tile servers before analysis
terrascope.cleanup_tile_servers()

first_item = items[0]
print(f"Analyzing: {first_item.datetime.date()}")

with rioxarray.open_rasterio(first_item.assets["NDVI"].href, mask_and_scale=True) as ds:
    clipped = ds.rio.clip_box(*bbox, crs="EPSG:4326")
    data = clipped.sel(band=1).values

    print(f"\nNDVI Statistics:")
    print(f"  Min: {np.nanmin(data):.2f}")
    print(f"  Max: {np.nanmax(data):.2f}")
    print(f"  Mean: {np.nanmean(data):.2f}")

Explore Available Collections

1
2
3
4
collections = terrascope.list_collections()
print(f"Available collections ({len(collections)}):")
for c in sorted(collections):
    print(f"  {c}")

ESA WorldCover Visualization

Search and visualize the ESA WorldCover 10m 2021 v2 global land cover map.

1
2
3
4
items = terrascope.search(
    collection="esa-worldcover-map-10m-2021-v2",
    bbox=bbox,
)

Inspect the First Item

1
items[0]

Visualize on Map

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
m = leafmap.Map(center=center, zoom=14)
m.add_basemap("Esri.WorldImagery")
m.add_raster(
    items[0].assets["Map"].href,
    layer_name="Land cover",
)
m.add_legend(
    title="Land Cover Type",
    builtin_legend="ESA_WorldCover",
)
m

Clean Up

Optionally log out to clear GDAL authentication credentials from the environment.

1
2
# Optional: logout when done
# terrascope.logout()