Satellite Imagery Analysis — Deep Dive

Moving from notebook-level NDVI maps to production satellite pipelines requires mastering data cube construction, atmospheric correction, multi-temporal compositing, machine-learning classification, and cloud-native I/O. This deep dive covers the patterns that scale.

Building analysis-ready data cubes

Raw satellite scenes arrive as individual granules with varying footprints, CRS, and cloud cover. The first step is assembling them into a regularly gridded data cube.

stackstac: STAC items → xarray

import stackstac
from pystac_client import Client

catalog = Client.open("https://earth-search.aws.element84.com/v1")
search = catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=[-74.05, 40.65, -73.85, 40.85],
    datetime="2025-06-01/2025-08-31",
    query={"eo:cloud_cover": {"lt": 20}},
)
items = list(search.items())

stack = stackstac.stack(
    items,
    assets=["B04", "B08", "SCL"],  # Red, NIR, Scene Classification
    resolution=10,
    epsg=32618,
    bounds_latlon=[-74.05, 40.65, -73.85, 40.85],
)
print(stack)  # <xarray.DataArray (time, band, y, x)>

stackstac lazily constructs a Dask-backed xarray DataArray. No data is downloaded until you call .compute() on a subset. This lets you define processing on a cube spanning hundreds of scenes without running out of memory.

Atmospheric correction

Level-1 (L1C) data measures top-of-atmosphere reflectance — what the satellite sensor recorded. Level-2 (L2A) applies atmospheric correction to estimate bottom-of-atmosphere reflectance — what the ground actually reflects.

For Sentinel-2, the ESA provides L2A products processed with Sen2Cor. For Landsat, USGS provides Collection 2 Level-2 products. When L2A is unavailable, Python options include:

# py6s wraps the 6S radiative transfer model
from Py6S import SixS, AtmosProfile, Wavelength

s = SixS()
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer)
s.wavelength = Wavelength(0.640, 0.670)  # Red band
s.run()
correction_factor = s.outputs.atmos_corrected_reflectance

For most use cases, start with pre-corrected L2A products and only calibrate manually when working with sensors that lack standard processing.

Cloud masking: beyond SCL

The default Scene Classification Layer misses thin cirrus and cloud edges. Production pipelines use multiple strategies:

S2Cloudless (ML-based)

from s2cloudless import S2PixelCloudDetector

detector = S2PixelCloudDetector(threshold=0.4, all_bands=True)
cloud_probs = detector.get_cloud_probability_maps(band_stack)
cloud_mask = cloud_probs > 0.5

Temporal compositing (median pixel)

Instead of masking individual scenes, compute the median across a time window. Clouds are transient, so the median eliminates them:

# stack shape: (time, band, y, x) — Dask array
ndvi_stack = (stack.sel(band="B08") - stack.sel(band="B04")) / \
             (stack.sel(band="B08") + stack.sel(band="B04") + 1e-10)

# Mask cloudy pixels with NaN, then take temporal median
scl = stack.sel(band="SCL")
cloud_free = ndvi_stack.where(~scl.isin([3, 8, 9]))
median_ndvi = cloud_free.median(dim="time").compute()

Median composites are the workhorse of large-area land-cover mapping.

Time-series analysis

Satellite time-series reveal seasonal patterns and anomalies:

import xarray as xr

# Monthly NDVI means
monthly = ndvi_stack.resample(time="1ME").mean()

# Detect anomalies (deviation from long-term monthly mean)
climatology = monthly.groupby("time.month").mean(dim="time")
anomalies = monthly.groupby("time.month") - climatology

Agricultural monitoring uses these anomalies to detect drought stress weeks before it becomes visible at ground level.

Phenology extraction

from scipy.signal import savgol_filter

# Smooth the NDVI time-series for a single pixel
pixel_ts = ndvi_stack.sel(x=-73.95, y=40.72, method="nearest").values
smoothed = savgol_filter(pixel_ts, window_length=7, polyorder=2)

# Green-up date: first day NDVI exceeds 0.4
green_up_idx = np.argmax(smoothed > 0.4)

Machine learning classification

Pixel-based classification

from sklearn.ensemble import RandomForestClassifier
import numpy as np

# Training data: labeled pixels with band values
# X shape: (n_samples, n_bands), y shape: (n_samples,)
X_train = np.column_stack([red_train, green_train, nir_train, swir_train])
y_train = labels  # 0=water, 1=forest, 2=urban, 3=agriculture

clf = RandomForestClassifier(n_estimators=200, n_jobs=-1)
clf.fit(X_train, y_train)

# Predict full image
bands = np.stack([red, green, nir, swir], axis=-1)
flat = bands.reshape(-1, 4)
prediction = clf.predict(flat).reshape(red.shape)

Deep learning with torchgeo

For state-of-the-art accuracy, use pretrained models:

from torchgeo.models import ResNet18_Weights

weights = ResNet18_Weights.SENTINEL2_ALL_MOCO
model = timm.create_model("resnet18", in_chans=13, num_classes=10)
model.load_state_dict(weights.get_state_dict())

torchgeo provides pretrained weights for Sentinel-2, Landsat, and other sensors, significantly reducing training data requirements.

Zonal statistics

Aggregate raster values within vector boundaries — for example, average NDVI per farm field:

from rasterstats import zonal_stats

stats = zonal_stats(
    "fields.geojson",
    "ndvi.tif",
    stats=["mean", "std", "min", "max", "count"],
)
# [{'mean': 0.72, 'std': 0.08, 'min': 0.3, 'max': 0.88, 'count': 4521}, ...]

This bridges raster analysis and vector reporting.

Scalable processing patterns

Tiled processing with Dask

import dask.array as da

# Process 1TB of imagery in 256×256 tiles
result = da.map_blocks(
    compute_ndvi,
    stack.sel(band="B04").data,
    stack.sel(band="B08").data,
    dtype=np.float32,
)
result.to_zarr("ndvi_cube.zarr")

Cloud-native pipeline architecture

STAC Search → stackstac (lazy cube) → Dask processing → COG/Zarr output → STAC catalog
     ↓                                      ↓
  pystac-client                       Coiled/Dask Gateway cluster

This architecture processes continental-scale datasets by distributing tile reads and computations across cloud workers. Data never leaves the cloud — reads go from S3 COGs to compute nodes, and outputs write back to S3.

Common pitfalls

Pitfall	Impact	Solution
Mixing L1C and L2A data	Inconsistent reflectance values	Use L2A exclusively or apply consistent correction
Ignoring nodata values	Corrupted statistics	Always mask with src.read(masked=True) or explicit nodata checks
Wrong band ordering	NDVI inverted or nonsensical	Verify band assignments against product documentation
Computing area in degrees	Meaningless pixel areas	Reproject to UTM before area calculations
Not clipping to AOI	Processing unnecessary data	Use windowed reads or bbox filtering
Single-date analysis	Clouds, shadows, seasonality	Use multi-temporal composites

Real-world applications

• Precision agriculture: NDVI time-series per field → variable-rate fertilizer maps
• Deforestation monitoring: Monthly change detection on Sentinel-2 → alert system (Global Forest Watch uses this approach)
• Urban heat islands: Landsat thermal band → surface temperature maps → city planning
• Flood mapping: SAR (Sentinel-1) backscatter change → water extent during events
• Carbon estimation: Biomass indices + LiDAR calibration → forest carbon stock maps

The one thing to remember: Production satellite analysis is not about individual scenes — it is about building lazy, cloud-native data cubes from STAC catalogs, applying temporal composites to eliminate noise, and scaling computation with Dask so that the same code works for one field or one continent.