STAC-Based Exploration and Visualization of Sentinel-2 and Sentinel-1 Data
Learn how to access Sentinel-2 and Sentinel-1 images via the EOPF STAC Catalog using pystac-client
Table of Contents¶
Introduction¶
This notebook provides a comprehensive walkthrough for accessing, filtering, and visualizing Sentinel-1 and Sentinel-2 satellite imagery using a STAC (SpatioTemporal Asset Catalog) interface. Leveraging the pystac-client, xarray, and matplotlib libraries, the notebook demonstrates how to:
Connect to a public EOPF STAC catalog hosted at https://
stac .core .eopf .eodc .eu. Perform structured searches across Sentinel-2 L2A and Sentinel-1 GRD collections, filtering by spatial extent, date range, cloud cover, and orbit characteristics.
Retrieve and load data assets directly into xarray for interactive analysis.
Visualize Sentinel-2 RGB composites along with pixel-level cloud coverage masks.
Render Sentinel-1 backscatter data (e.g., VH polarization) for selected acquisitions.
Setup¶
Start importing the necessary libraries
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
import numpy as np
import pystac_client
import xarray as xr
from pystac_client import CollectionSearch
from matplotlib.gridspec import GridSpecSTAC Collection Discovery¶
Using CollectionSearch to list available STAC collections
# Initialize the collection search
search = CollectionSearch(
url="https://stac.core.eopf.eodc.eu/collections", # STAC /collections endpoint
)
# Retrieve all matching collections (as dictionaries)
for collection_dict in search.collections_as_dicts():
print(collection_dict["id"])sentinel-2-l2a
sentinel-3-slstr-l1-rbt
sentinel-3-olci-l2-lfr
sentinel-2-l1c
sentinel-3-slstr-l2-lst
sentinel-1-l1-slc
sentinel-3-olci-l1-efr
sentinel-3-olci-l1-err
sentinel-1-l2-ocn
sentinel-1-l1-grd
sentinel-3-olci-l2-lrr
/home/sdhinakaran/micromamba/envs/new_eopf/lib/python3.11/site-packages/pystac_client/collection_search.py:121: PystacClientWarning: Unable to parse extent from collection=sentinel-3-olci-l2-lrr
warnings.warn(
Sentinel-2 Item Search¶
Querying the Sentinel-2 L2A collection by bounding box, date range, and cloud cover
catalog = pystac_client.Client.open("https://stac.core.eopf.eodc.eu")
# Search with cloud cover filter
items = list(
catalog.search(
collections=["sentinel-2-l2a"],
bbox=[7.2, 44.5, 7.4, 44.7],
datetime=["2025-01-30", "2025-05-01"],
query={"eo:cloud_cover": {"lt": 20}}, # Cloud cover less than 20%
).items()
)
print(items)[<Item id=S2B_MSIL2A_20250430T101559_N0511_R065_T32TLQ_20250430T131328>, <Item id=S2C_MSIL2A_20250425T102041_N0511_R065_T32TLQ_20250425T155812>, <Item id=S2C_MSIL2A_20250418T103041_N0511_R108_T32TLQ_20250418T160655>, <Item id=S2C_MSIL2A_20250405T102041_N0511_R065_T32TLQ_20250405T175414>]
Quicklook Visualization for Sentinel-2¶
We can use the RGB quicklook layer contained in the Sentinel-2 EOPF Zarr product for a visualization of the content:
item = items[0] # extracting the first item
ds = xr.open_dataset(
item.assets["product"].href,
**item.assets["product"].extra_fields["xarray:open_datatree_kwargs"],
) # The engine="eopf-zarr" is already embedded in the STAC metadata
ds.quality_l2a_quicklook_r60m_tci.plot.imshow(rgb="quality_l2a_quicklook_r60m_band")
We can view side by side the RGB quicklook and the SCL (Scene Classification Layer) for all the returned STAC Items:
def visualize_items_table(items):
"""
Visualize multiple Sentinel-2 items with scene IDs centered above each row,
followed by two columns (RGB and SCL).
Parameters:
items (list): List of STAC items
"""
num_items = len(items)
row_height = 7 # Approximate height per item
fig_height = min(row_height * num_items, 24)
fig = plt.figure(figsize=(14, fig_height))
# Create height ratios: 0.3 for title row, 1.7 for image row (per item)
height_ratios = []
for _ in range(num_items):
height_ratios.extend([0.3, 1.7]) # first is title, second is images
gs = GridSpec(
num_items * 2, 2, height_ratios=height_ratios, hspace=0.1, wspace=0.05
)
# Create colorbar axis (positioned below all plots)
scl_cbar_ax = fig.add_axes([0.3, 0.03, 0.4, 0.015]) # [left, bottom, width, height]
for i, item in enumerate(items):
try:
# Add centered title for the row
title_ax = fig.add_subplot(gs[i * 2, :]) # Span both columns
title_text = f"{item.id} (Cloud: {item.properties['eo:cloud_cover']:.1f}%)"
title_ax.text(
0.5,
0.5,
title_text,
ha="center",
va="center",
fontsize=10,
fontweight="bold",
)
title_ax.axis("off")
# Load the dataset
ds = xr.open_dataset(
item.assets["product"].href,
**item.assets["product"].extra_fields["xarray:open_datatree_kwargs"],
)
# RGB Column
ax1 = fig.add_subplot(gs[i * 2 + 1, 0])
try:
if "quality_l2a_quicklook_r60m_tci" in ds:
ds.quality_l2a_quicklook_r60m_tci.plot.imshow(ax=ax1)
ax1.set_title("True Color (TCI)", pad=5, fontsize=9)
else:
ax1.text(
0.5,
0.5,
"No RGB quicklook",
ha="center",
va="center",
fontsize=8,
)
except Exception:
ax1.text(0.5, 0.5, "RGB Error", ha="center", va="center", fontsize=8)
ax1.axis("off")
# SCL Column
ax2 = fig.add_subplot(gs[i * 2 + 1, 1])
try:
if "conditions_mask_l2a_classification_r60m_scl" in ds:
scl = ds.conditions_mask_l2a_classification_r60m_scl.compute()
# Apply classification attributes
scl.attrs.update(
{
"flag_values": [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11],
"flag_meanings": "no_data sat_or_defect_pixel topo_casted_shadows cloud_shadows vegetation not_vegetation water unclassified cloud_medium_prob cloud_high_prob thin_cirrus snow_or_ice",
"flag_colors": "#000000 #ff0000 #2f2f2f #643200 #00a000 #ffe65a #0000ff #808080 #c0c0c0 #ffffff #64c8ff #ff96ff",
}
)
cmap = mcolors.ListedColormap(scl.attrs["flag_colors"].split(" "))
norm = mcolors.BoundaryNorm(
boundaries=np.arange(13) - 0.5, ncolors=12 # For 12 classes
)
img_scl = scl.plot.imshow(
ax=ax2, cmap=cmap, norm=norm, add_colorbar=False
)
ax2.set_title("Scene Classification", pad=5, fontsize=9)
# Create colorbar once
if i == 0:
cbar = fig.colorbar(
img_scl,
cax=scl_cbar_ax,
orientation="horizontal",
ticks=np.arange(12),
)
class_labels = [
label.replace("_", " ").title()
for label in scl.attrs["flag_meanings"].split()
]
cbar.ax.set_xticklabels(
class_labels, rotation=45, ha="right", fontsize=8
)
cbar.set_label("Land Cover Classes", fontsize=9)
else:
ax2.text(
0.5, 0.5, "No SCL layer", ha="center", va="center", fontsize=8
)
except Exception:
ax2.text(0.5, 0.5, "SCL Error", ha="center", va="center", fontsize=8)
ax2.axis("off")
except Exception as e:
print(f"Error processing item {item.id}: {str(e)}")
continue
plt.tight_layout()
plt.subplots_adjust(bottom=0.07, top=0.95) # Adjust for colorbar and titles
plt.show()
visualize_items_table(items)/tmp/ipykernel_358028/1482541263.py:122: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
plt.tight_layout()
Sentinel-1 GRD Item Search¶
Querying Sentinel-1 L1 GRD collection given bounding box, orbit direction and orbit track number.
# Load the STAC catalog
catalog = pystac_client.Client.open("https://stac.core.eopf.eodc.eu")
# Search for Sentinel-1 GRD items matching the metadata characteristics
items = list(
catalog.search(
collections=["sentinel-1-l1-grd"],
bbox=[29.922651, -9.648995, 32.569355, -7.391839],
datetime=[
"2025-06-09T16:19:07Z",
"2025-06-09T16:19:40Z",
],
query={
"sat:orbit_state": {"eq": "ascending"},
"sat:relative_orbit": {"eq": 174},
},
).items()
)
items[<Item id=S1A_IW_GRDH_1SDV_20250609T161936_20250609T162001_059571_07655C_4E9C>,
<Item id=S1A_IW_GRDH_1SDV_20250609T161907_20250609T161936_059571_07655C_0AB1>]Sentinel-1 GRD Visualization¶
Opening the first STAC Item and plot the VH polarization
item = items[1]
ds = xr.open_dataset(
item.assets["product"].href,
**item.assets["product"].extra_fields["xarray:open_datatree_kwargs"],
variables="VH_measurements*",
)
plt.imshow(ds["VH_measurements_grd"][::10, ::10].to_numpy(), vmin=0, vmax=200)
import matplotlib.pyplot as plt
import xarray as xr
def visualize_sentinel1_items(items, downsample_factor=10, vmin=0, vmax=200):
num_items = len(items)
if num_items == 0:
print("No items to visualize")
return
ncols = 2
nrows = (num_items + 1) // 2
fig = plt.figure(figsize=(12, 4 * nrows + 1))
gs = GridSpec(nrows + 1, ncols, height_ratios=[1] * nrows + [0.15], hspace=0.4)
images = []
for idx, item in enumerate(items):
row = idx // 2
if num_items % 2 == 1 and idx == num_items - 1: # Center last plot if odd
ax = fig.add_subplot(gs[row, :])
else:
col = idx % 2
ax = fig.add_subplot(gs[row, col])
ds = xr.open_dataset(
item.assets["product"].href,
**item.assets["product"].extra_fields["xarray:open_datatree_kwargs"],
variables="VH_measurements*",
)
img = ax.imshow(
ds["VH_measurements_grd"][
::downsample_factor, ::downsample_factor
].to_numpy(),
vmin=vmin,
vmax=vmax,
)
images.append(img)
ax.set_title(
f"{item.id}\nOrbit {item.properties.get('sat:relative_orbit', 'N/A')} - {item.properties.get('sat:orbit_state', 'N/A').capitalize()}",
fontsize=8,
)
ax.axis("off")
cbar_ax = fig.add_subplot(gs[-1, :])
fig.colorbar(
images[0], cax=cbar_ax, orientation="horizontal", label="VH Backscatter (dB)"
)
cbar_ax.xaxis.set_label_position("top")
cbar_ax.xaxis.set_ticks_position("top")
fig.suptitle(
f"Sentinel-1 GRD VH Measurements ({num_items} scenes)", y=1.02, fontsize=13
)
plt.show()
visualize_sentinel1_items(items)
The Sentinel-1 visualization is off due to the EOPF CPM version used for the data conversion. Once that the following issue will be solved we will update the STAC Catalog and this notebook to reflect the updates.