Module 3

Data Access & APIs

Programmatic access to satellite imagery through Copernicus, USGS, Google Earth Engine, and STAC APIs

Note: The code examples below are illustrative and require external internet access and valid API credentials to execute. They demonstrate the standard workflows for each data access method. Register for free accounts at the respective platforms to run them.

1. Copernicus Data Space Ecosystem (CDSE)

The Copernicus Data Space Ecosystem (dataspace.copernicus.eu) is the official distribution platform for all Sentinel data. It replaced the legacy Copernicus Open Access Hub in October 2023. Access is free after registration.

Key Features

• Full Sentinel-1, 2, 3, 5P archive (PB onwards)
• OData and OpenSearch catalogue APIs
• Sentinel Hub integration for on-the-fly processing
• JupyterHub environment with pre-installed tools
• S3 object storage access for batch processing

Searching & Downloading with sentinelsat

The sentinelsat library provides a Python interface to the Copernicus catalogue. Note: as of 2024, sentinelsat supports the new CDSE endpoints.

Copernicus CDSE Sentinel-2 Search Demo

Python

script.py113 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from datetime import datetime, timedelta

print("=" * 62)
print("  Copernicus CDSE - Sentinel-2 Search Demo")
print("  (Simulated API response - real usage requires credentials)")
print("=" * 62)
print()

# ── Simulated Sentinel-2 product catalogue response ──
# In production: from sentinelsat import SentinelAPI
# api = SentinelAPI('user', 'pass', 'https://catalogue.dataspace.copernicus.eu')
# products = api.query(area=wkt, date=(...), platformname='Sentinel-2', ...)

aoi = "Athens, Greece: POLYGON((23.6 37.9, 23.8 37.9, 23.8 38.1, 23.6 38.1, 23.6 37.9))"
print(f"AOI: {aoi}")
print(f"Date range: 2024-06-01 to 2024-06-30")
print(f"Platform: Sentinel-2, Level-2A, Cloud < 20%")
print()

# Simulate catalogue search results
np.random.seed(42)
n_products = 8
dates = [datetime(2024, 6, 1) + timedelta(days=np.random.randint(0, 30)) for _ in range(n_products)]
dates.sort()
tiles = ['34SGH', '34SGJ', '34SFH', '34SGH', '34SGJ', '34SFH', '34SGH', '34SGJ']
clouds = np.random.uniform(0, 20, n_products)
sizes_mb = np.random.uniform(600, 900, n_products)

print(f"Found {n_products} products:")
print(f"{'Title':<55s} {'Cloud%':>7s} {'Size':>8s}")
print("-" * 72)
for i in range(n_products):
    title = f"S2A_MSIL2A_{dates[i].strftime('%Y%m%d')}T092031_N0510_R093_{tiles[i]}_{dates[i].strftime('%Y%m%d')}T120000"
    print(f"{title:<55s} {clouds[i]:6.1f}% {sizes_mb[i]:7.0f}MB")

best = np.argmin(clouds)
print(f"\nBest scene: cloud cover = {clouds[best]:.1f}%")
print(f"  -> S2A_MSIL2A_{dates[best].strftime('%Y%m%d')}..._{tiles[best]}")

# ── Sentinel-2 band information ──
print("\n── Sentinel-2 MSI Band Configuration ──")
bands = [
    ("B01", 443, 60, "Coastal aerosol"),
    ("B02", 490, 10, "Blue"),
    ("B03", 560, 10, "Green"),
    ("B04", 665, 10, "Red"),
    ("B05", 705, 20, "Red-Edge 1"),
    ("B06", 740, 20, "Red-Edge 2"),
    ("B07", 783, 20, "Red-Edge 3"),
    ("B08", 842, 10, "NIR"),
    ("B8A", 865, 20, "NIR narrow"),
    ("B09", 945, 60, "Water vapour"),
    ("B10", 1375, 60, "Cirrus"),
    ("B11", 1610, 20, "SWIR 1"),
    ("B12", 2190, 20, "SWIR 2"),
]
print(f"{'Band':<6s} {'Centre nm':>10s} {'Res m':>6s} {'Purpose'}")
for b, wl, res, desc in bands:
    print(f"{b:<6s} {wl:>10d} {res:>6d} {desc}")

# ── Plot: cloud cover distribution + timeline ──
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
fig.patch.set_facecolor('#0f172a')
fig.suptitle('Sentinel-2 L2A Catalogue Search: Athens, June 2024', color='#7dd3fc', fontsize=13)

# Cloud cover histogram
ax1.bar(range(n_products), clouds, color='#0ea5e9', edgecolor='#0369a1', width=0.6)
ax1.axhline(y=20, color='#ef4444', linestyle='--', alpha=0.7, label='Cloud threshold (20%)')
ax1.bar(best, clouds[best], color='#22c55e', edgecolor='#166534', width=0.6)
ax1.set_xlabel('Product index', color='white')
ax1.set_ylabel('Cloud cover (%)', color='white')
ax1.set_title('Cloud Cover per Scene', color='#bae6fd')
ax1.legend(fontsize=9, facecolor='#1e293b', edgecolor='#334155', labelcolor='white')
ax1.set_facecolor('#0f172a')
ax1.tick_params(colors='white')
for s in ax1.spines.values(): s.set_color('#334155')

# Band wavelength plot
wls = [b[1] for b in bands]
ress = [b[2] for b in bands]
colors = ['#ef4444' if r==60 else '#f59e0b' if r==20 else '#22c55e' for r in ress]
ax2.scatter(wls, [10]*len(wls), c=colors, s=80, zorder=5, edgecolors='white', linewidths=0.5)
for b, wl, res, desc in bands:
    ax2.annotate(b, (wl, 10), textcoords="offset points", xytext=(0, 12),
                fontsize=7, color='#7dd3fc', ha='center')
ax2.axvspan(400, 700, alpha=0.1, color='#22c55e', label='Visible')
ax2.axvspan(700, 1000, alpha=0.1, color='#ef4444', label='NIR')
ax2.axvspan(1000, 2500, alpha=0.1, color='#f59e0b', label='SWIR')
ax2.set_xlabel('Wavelength (nm)', color='white')
ax2.set_title('Sentinel-2 Band Positions', color='#bae6fd')
ax2.set_xlim(350, 2400)
ax2.set_ylim(8, 14)
ax2.set_yticks([])
ax2.legend(fontsize=8, facecolor='#1e293b', edgecolor='#334155', labelcolor='white', loc='upper right')
ax2.set_facecolor('#0f172a')
ax2.tick_params(colors='white')
for s in ax2.spines.values(): s.set_color('#334155')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

print("\n── How to run with real credentials ──")
print("from sentinelsat import SentinelAPI")
print("api = SentinelAPI('your_email', 'your_password',")
print("                  'https://catalogue.dataspace.copernicus.eu')")
print("products = api.query(area=wkt, date=('20240601','20240630'),")
print("                     platformname='Sentinel-2', processinglevel='Level-2A')")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

2. USGS EarthExplorer & Landsat

The USGS EarthExplorer (earthexplorer.usgs.gov) hosts the complete Landsat archive plus many other datasets. Thelandsatxplore library provides programmatic access via the USGS M2M API.

Available Datasets

Landsat Collection 2

L1-L9, Level-1 & Level-2 (1972-present)

Sentinel-2

Mirror of Copernicus archive

ASTER

14-band VNIR/SWIR/TIR, 15-90 m

SRTM / NASADEM

Global DEM at 30 m / 1 arc-second

Searching Landsat with landsatxplore

USGS EarthExplorer Landsat Search Demo

Python

script.py76 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

print("=" * 62)
print("  USGS EarthExplorer - Landsat Search Demo")
print("  (Simulated - real usage requires USGS EROS credentials)")
print("=" * 62)
print()

# In production:
# from landsatxplore.api import API
# api = API('username', 'password')
# scenes = api.search(dataset='landsat_ot_c2_l2', latitude=37.97, ...)

# ── Simulated Landsat catalogue response ──
np.random.seed(123)
scenes = []
for i in range(6):
    day = 3 + i * 5
    sat = np.random.choice([8, 9])
    path, row = 183, 33
    cloud = np.random.uniform(0, 15)
    scenes.append({
        'display_id': f'LC0{sat}_L2SP_{path:03d}{row:03d}_20240{6}{day:02d}_02_T1',
        'date': f'2024-06-{day:02d}',
        'cloud': cloud,
        'satellite': sat,
        'path': path, 'row': row,
    })

print(f"AOI: Athens, Greece (37.97N, 23.73E)")
print(f"Dataset: landsat_ot_c2_l2 (Landsat 8/9 Collection 2 Level 2)")
print(f"Found {len(scenes)} scenes:\n")
for s in scenes:
    print(f"  {s['display_id']}")
    print(f"    Date: {s['date']}  Cloud: {s['cloud']:.1f}%  Landsat {s['satellite']}")

# ── Landsat 8/9 band configuration ──
print("\n── Landsat 8/9 OLI/TIRS Band Configuration ──")
l8_bands = [
    ("B1",  433, 453, 30, "Coastal/Aerosol"),
    ("B2",  450, 515, 30, "Blue"),
    ("B3",  525, 600, 30, "Green"),
    ("B4",  630, 680, 30, "Red"),
    ("B5",  845, 885, 30, "NIR"),
    ("B6",  1560, 1660, 30, "SWIR-1"),
    ("B7",  2100, 2300, 30, "SWIR-2"),
    ("B8",  500, 680, 15, "Pan"),
    ("B9",  1360, 1390, 30, "Cirrus"),
    ("B10", 10600, 11200, 100, "TIR-1"),
    ("B11", 11500, 12500, 100, "TIR-2"),
]
print(f"{'Band':<5s} {'Range (nm)':>12s} {'Res':>5s} {'Purpose'}")
for b, lo, hi, res, desc in l8_bands:
    print(f"{b:<5s} {lo:>5d}-{hi:<5d} {res:>4d}m {desc}")

# Plot
fig, ax = plt.subplots(figsize=(10, 5))
fig.patch.set_facecolor('#0f172a')
clouds = [s['cloud'] for s in scenes]
labels = [s['date'] for s in scenes]
colors = ['#22c55e' if c < 5 else '#f59e0b' if c < 10 else '#ef4444' for c in clouds]
ax.barh(labels, clouds, color=colors, edgecolor='#1e293b', height=0.5)
ax.axvline(x=15, color='#ef4444', linestyle='--', alpha=0.5, label='Cloud threshold')
ax.set_xlabel('Cloud Cover (%)', color='white')
ax.set_title('Landsat Scenes: Athens, June 2024', color='#fde68a')
ax.set_facecolor('#0f172a')
ax.tick_params(colors='white')
ax.legend(fontsize=9, facecolor='#1e293b', edgecolor='#334155', labelcolor='white')
for s in ax.spines.values(): s.set_color('#334155')
plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

Click Run to execute the Python code

Code will be executed with Python 3 on the server

3. Google Earth Engine (GEE)

Google Earth Engine is a cloud-based geospatial analysis platform with a multi-petabyte catalog of satellite imagery and geospatial datasets. Processing runs on Google's infrastructure, eliminating the need to download data locally.

Key Advantages

• Server-side computation: analyze terabytes without downloading
• Massive catalog: Landsat, Sentinel, MODIS, ERA5, and 900+ datasets
• Temporal compositing: easily create mosaics, time series, anomalies
• Python API (ee) + JavaScript Code Editor + geemap for visualization
• Free for research, education, and non-commercial use

NDVI Time Series with GEE Python API

Google Earth Engine NDVI Demo

Python

script.py75 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

print("=" * 62)
print("  Google Earth Engine - NDVI Time Series Demo")
print("  (Simulated - real GEE requires authentication)")
print("=" * 62)
print()

# In production:
# import ee; ee.Initialize(project='your-project')
# s2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
#       .filterBounds(aoi).filterDate(...).filter(cloud < 20)

# ── Simulated monthly NDVI for Nile Delta ──
months = np.arange(1, 13)
month_names = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']

# Nile Delta agriculture: high NDVI in winter (wheat), low in summer (fallow)
ndvi_nile = 0.45 + 0.25 * np.cos(2*np.pi*(months - 2)/12) + np.random.normal(0, 0.02, 12)
ndvi_nile = np.clip(ndvi_nile, 0.1, 0.8)

print("AOI: Nile Delta, Egypt (30.5-31.5E, 30.8-31.4N)")
print("Dataset: COPERNICUS/S2_SR_HARMONIZED (2023)")
print("Filter: cloud < 20%")
print()
print("Monthly Mean NDVI:")
print("-" * 40)
for i, (m, v) in enumerate(zip(month_names, ndvi_nile)):
    bar = '#' * int(v * 40)
    print(f"  {m}: NDVI = {v:.3f} {bar}")

# Plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
fig.patch.set_facecolor('#0f172a')
fig.suptitle('Google Earth Engine: NDVI Time Series (Nile Delta 2023)', color='#7dd3fc', fontsize=13)

# Monthly NDVI
ax1.fill_between(months, ndvi_nile, alpha=0.3, color='#22c55e')
ax1.plot(months, ndvi_nile, 'o-', color='#22c55e', linewidth=2, markersize=6)
ax1.axhline(y=0.4, color='#f59e0b', linestyle='--', alpha=0.5, label='Vegetation threshold')
ax1.set_xticks(months)
ax1.set_xticklabels(month_names, fontsize=8)
ax1.set_ylabel('NDVI', color='white')
ax1.set_title('Monthly Mean NDVI', color='#86efac')
ax1.legend(fontsize=9, facecolor='#1e293b', edgecolor='#334155', labelcolor='white')
ax1.set_facecolor('#0f172a')
ax1.tick_params(colors='white')
for s in ax1.spines.values(): s.set_color('#334155')

# Simulated NDVI map
np.random.seed(42)
ndvi_map = np.random.uniform(0.2, 0.7, (50, 60))
ndvi_map[10:40, 15:50] += 0.15  # agricultural area
ndvi_map[:8, :] = np.random.uniform(0.0, 0.15, (8, 60))  # water
ndvi_map = np.clip(ndvi_map, 0, 0.8)
im = ax2.imshow(ndvi_map, cmap='RdYlGn', vmin=0, vmax=0.8, aspect='auto')
ax2.set_title('NDVI Composite (median)', color='#86efac')
ax2.set_xlabel('Longitude pixels', color='white')
ax2.set_ylabel('Latitude pixels', color='white')
ax2.set_facecolor('#0f172a')
ax2.tick_params(colors='white')
plt.colorbar(im, ax=ax2, label='NDVI', shrink=0.8)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

print("\n── GEE Code Pattern (requires auth) ──")
print("import ee; ee.Initialize(project='...')")
print("s2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')")
print("ndvi = s2.map(lambda img: img.normalizedDifference(['B8','B4']))")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

4. STAC API (SpatioTemporal Asset Catalog)

STAC is an open specification for describing geospatial data assets. It provides a unified, cloud-native interface for discovering and accessing satellite imagery across different providers.

STAC Concepts

Catalog

Top-level container linking to collections

Collection

Group of related items (e.g., Sentinel-2 L2A)

Item

Single spatiotemporal data entity (one scene)

Asset

Actual file (COG, NetCDF, etc.) linked to an item

Major STAC Endpoints

Provider	URL	Key Data
Element 84	earth-search.aws.element84.com/v1	S2 L2A, Landsat C2, NAIP
Microsoft PC	planetarycomputer.microsoft.com/api/stac/v1	S1, S2, Landsat, DEM, ERA5
CDSE	catalogue.dataspace.copernicus.eu/stac	All Sentinel data

Cloud-Native Workflow with pystac-client + odc-stac

STAC API Cloud-Native Data Access Demo

Python

script.py99 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import json

print("=" * 62)
print("  STAC API - Cloud-Native Data Access Demo")
print("  (Simulated - real STAC requires internet access)")
print("=" * 62)
print()

# In production:
# from pystac_client import Client
# catalog = Client.open('https://earth-search.aws.element84.com/v1')
# search = catalog.search(collections=['sentinel-2-l2a'], bbox=[...], ...)

# ── Simulated STAC catalogue ──
stac_endpoints = [
    ("Element 84 Earth Search", "earth-search.aws.element84.com/v1", "Sentinel-2, Landsat, NAIP"),
    ("Microsoft Planetary Computer", "planetarycomputer.microsoft.com/api/stac/v1", "S2, L8/9, MODIS, ERA5"),
    ("Google Earth Engine STAC", "earthengine-stac.storage.googleapis.com", "All GEE collections"),
    ("Copernicus CDSE", "catalogue.dataspace.copernicus.eu/stac", "All Sentinel missions"),
]

print("Major STAC Endpoints:")
for name, url, desc in stac_endpoints:
    print(f"  {name}")
    print(f"    URL: https://{url}")
    print(f"    Collections: {desc}")
    print()

# Simulated STAC search results
np.random.seed(55)
items = []
for i in range(6):
    day = 2 + i * 5
    cloud = np.random.uniform(0, 18)
    items.append({
        "id": f"S2A_34SGH_20240{6}{day:02d}_0_L2A",
        "datetime": f"2024-06-{day:02d}T09:20:31Z",
        "cloud_cover": round(cloud, 1),
        "assets": ["B02","B03","B04","B05","B06","B07","B08","B8A","B11","B12","SCL","visual"],
        "href": f"s3://sentinel-cogs/sentinel-s2-l2a-cogs/34/S/GH/2024/6/S2A_34SGH_20240{6}{day:02d}_0_L2A/"
    })

print(f"Search: sentinel-2-l2a, bbox=[23.6,37.9,23.8,38.1], 2024-06")
print(f"Found {len(items)} items:\n")

for item in items:
    print(f"  {item['id']}")
    print(f"    Date: {item['datetime']}  Cloud: {item['cloud_cover']}%")
    print(f"    Assets: {item['assets'][:6]}...")
    print(f"    COG: {item['href']}B04.tif")
    print()

# ── STAC Item GeoJSON structure ──
stac_item = {
    "type": "Feature",
    "stac_version": "1.0.0",
    "id": items[0]["id"],
    "geometry": {"type": "Polygon", "coordinates": [[[23.6,37.9],[23.8,37.9],[23.8,38.1],[23.6,38.1],[23.6,37.9]]]},
    "properties": {
        "datetime": items[0]["datetime"],
        "eo:cloud_cover": items[0]["cloud_cover"],
        "platform": "sentinel-2a",
        "instruments": ["msi"],
    },
    "assets": {k: {"href": f"{items[0]['href']}{k}.tif", "type": "image/tiff; application=geotiff; profile=cloud-optimized"} for k in items[0]["assets"][:4]}
}
print("STAC Item GeoJSON (truncated):")
print(json.dumps(stac_item, indent=2)[:600] + "\n  ...")

# Plot: STAC workflow diagram
fig, ax = plt.subplots(figsize=(10, 5))
fig.patch.set_facecolor('#0f172a')
steps = ["STAC\nCatalogue", "Search\n(bbox, date)", "Filter\n(cloud<20%)", "Get\nAsset URLs", "Read COG\n(HTTP range)", "xarray\nDataset"]
x = np.arange(len(steps))
for i, step in enumerate(steps):
    color = '#0ea5e9' if i < 4 else '#22c55e'
    ax.add_patch(plt.Rectangle((i-0.35, 0.2), 0.7, 0.6, facecolor=color+'33', edgecolor=color, linewidth=2, zorder=2))
    ax.text(i, 0.5, step, ha='center', va='center', fontsize=10, color='white', fontweight='bold', zorder=3)
    if i < len(steps)-1:
        ax.annotate('', xy=(i+0.4, 0.5), xytext=(i+0.6, 0.5), arrowprops=dict(arrowstyle='->', color='#64748b', lw=2))
ax.set_xlim(-0.6, len(steps)-0.4)
ax.set_ylim(0, 1)
ax.axis('off')
ax.set_title('STAC + COG Cloud-Native Workflow', color='#7dd3fc', fontsize=14, fontweight='bold')
ax.set_facecolor('#0f172a')
plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

print("\nKey advantages of STAC + COG:")
print("  1. No bulk downloads - read only the pixels you need via HTTP range requests")
print("  2. Lazy loading with xarray/dask for terabyte-scale analysis")
print("  3. Provider-agnostic: same code works with any STAC endpoint")
print("  4. Metadata standardized (GeoJSON + extensions)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Choosing the Right Data Access Method

Method	Best For	Pros	Cons
CDSE / sentinelsat	Sentinel full scenes	Official source, complete metadata	Requires download, large files
USGS / landsatxplore	Landsat, ASTER, DEM	50+ year Landsat archive	Download quotas, slower API
Google Earth Engine	Large-area analysis, time series	No downloads, massive compute	Proprietary API, export limits
STAC + COG	Cloud-native, custom pipelines	Standard API, partial reads	Requires cloud compute for speed

Recommendation: Start with Google Earth Engine for exploration and prototyping. Move to STAC + cloud-native workflows for production pipelines. Use CDSE/USGS when you need the full original data products with complete metadata.

Common EO Data Formats

GeoTIFF / COG

Standard raster format. Cloud-Optimized GeoTIFFs (COGs) support partial reads via HTTP range requests.

NetCDF / HDF5

Multi-dimensional scientific data. Used by MODIS, ERA5, ocean products. Supports metadata and compression.

SAFE (Sentinel)

Standard Archive Format for Europe. ZIP containing GeoTIFF bands, metadata XML, and quality masks.

Zarr

Cloud-native chunked array format. Ideal for large time series on S3/GCS. Used in Analysis-Ready Data cubes.

← M2: Physical Principles M4: Earthquake Monitoring →