TerraScope STAC API with leafmap¶

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

Setup: Set environment variables before running:

export TERRASCOPE_USERNAME='your_username'
export TERRASCOPE_PASSWORD='your_password'

Installation¶

Uncomment the following line to install leafmap if needed.

In [ ]:

# %pip install -U leafmap

# %pip install -U leafmap

Import libraries¶

In [ ]:

import leafmap
import leafmap.terrascope as terrascope

import leafmap import leafmap.terrascope as terrascope

Authentication¶

In [ ]:

terrascope.login()

terrascope.login()

Search for NDVI Data¶

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

In [ ]:

# 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}")

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

In [ ]:

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

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¶

In [ ]:

layers = terrascope.create_time_layers(items)

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

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¶

In [ ]:

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}")

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¶

In [ ]:

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

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.

In [ ]:

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

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

Inspect the First Item¶

In [ ]:

items[0]

items[0]

Visualize on Map¶

In [ ]:

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

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.

In [ ]:

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

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