Sentinel-2 & Landsat Cloud Masking Workflows #

In carbon accounting and MRV automation, optical satellite data remains the foundational input for land-use change detection, biomass estimation, and disturbance tracking. Atmospheric interference—particularly cloud cover, cirrus haze, and topographic shadow—introduces systematic bias into spectral indices and temporal composites. The Sentinel-2 & Landsat Cloud Masking Workflows stage sits at the critical intersection of raw ingestion and analytical readiness. When integrated into broader Satellite Imagery Processing for Emissions Tracking architectures, robust masking ensures that downstream algorithms operate on spectrally consistent, artifact-free surfaces. Without rigorous pre-processing, cloud leakage propagates through NDVI/EVI time series, inflates false-positive disturbance signals, and compromises compliance-grade carbon stock reporting.

Sensor-Specific QA Band Architecture #

Cloud masking is not a uniform threshold operation across optical constellations. Sentinel-2 Level-2A products distribute the Scene Classification Layer (SCL) and the QA60 bitmask, while Landsat 8/9 rely on the QA_PIXEL integer bitmask. Each encodes atmospheric classes using distinct bit structures and spatial resolutions:

• Sentinel-2 SCL: 20m native resolution. Classes include 3 (Cloud Shadow), 8 (Cloud Medium Probability), 9 (Cloud High Probability), and 10 (Thin Cirrus).
• Sentinel-2 QA60: 60m resolution. Bit 10 = opaque clouds, Bit 11 = cirrus.
• Landsat 8/9 QA_PIXEL: 30m resolution. Bits 3–4 encode cloud confidence, Bit 5 encodes cloud shadow, Bit 1 encodes cirrus.

Production pipelines must parse these structures using bitwise operations rather than naive value comparisons. Misinterpreting bit positions or ignoring confidence tiers leads to either over-masking (discarding usable pixels in humid tropics) or under-masking (retaining cloud-contaminated pixels that skew carbon flux calculations). Official specifications from ESA and USGS provide the definitive bitmaps required for accurate decoding: Sentinel-2 MSI Processing Levels and Landsat Collection 2 Quality Assessment Bands.

Production Implementation #

The following implementation demonstrates a production-ready cloud masking routine designed for MRV pipelines. It handles multi-sensor QA parsing, spatial alignment, error handling, and fallback routing when metadata is incomplete or corrupted. The workflow leverages rasterio for I/O and geospatial alignment, numpy for vectorized bitwise operations, xarray/dask for chunked memory management, and prefect for orchestration and audit logging.

import logging
import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.warp import calculate_default_transform, reproject
import xarray as xr
import dask.array as da
from prefect import flow, task
from dataclasses import dataclass
from typing import Optional, Tuple, Dict
from pathlib import Path

# Structured logging configuration for audit trails
logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s | %(levelname)s | %(name)s | %(message)s",
    datefmt="%Y-%m-%dT%H:%M:%S"
)
logger = logging.getLogger("mrv.cloud_masking")

@dataclass
class MaskConfig:
    cloud_threshold: float = 0.30  # Max acceptable cloud fraction
    shadow_include: bool = True
    fallback_skip: bool = False
    target_crs: str = "EPSG:4326"
    target_res: float = 10.0
    confidence_tiers: Tuple[int, ...] = (2, 3)  # Medium/High confidence for Landsat

def _parse_sentinel2_scl(scl_path: Path) -> np.ndarray:
    """Return boolean mask: True = valid, False = cloud/shadow/cirrus."""
    with rasterio.open(scl_path) as src:
        scl = src.read(1)
    # SCL classes: 3=shadow, 8=med_cloud, 9=high_cloud, 10=cirrus
    invalid = np.isin(scl, [3, 8, 9, 10])
    return ~invalid

def _parse_landsat_qa(qa_path: Path, config: MaskConfig) -> np.ndarray:
    """Parse Landsat QA_PIXEL using bitwise operations."""
    with rasterio.open(qa_path) as src:
        qa = src.read(1)
    # Bits: 1=cirrus, 3-4=cloud_confidence, 5=cloud_shadow
    cloud_conf = (qa >> 3) & 0b11
    shadow = (qa >> 5) & 0b1
    cirrus = (qa >> 1) & 0b1
    
    cloud_mask = np.isin(cloud_conf, config.confidence_tiers)
    shadow_mask = shadow == 1 if config.shadow_include else np.zeros_like(qa, dtype=bool)
    cirrus_mask = cirrus == 1
    
    return ~(cloud_mask | shadow_mask | cirrus_mask)

def _align_and_resample(mask: np.ndarray, src_profile: Dict, config: MaskConfig) -> np.ndarray:
    """Reproject mask to target CRS/resolution using nearest-neighbor."""
    transform, width, height = calculate_default_transform(
        src_profile["crs"], config.target_crs,
        src_profile["width"], src_profile["height"],
        *src_profile["bounds"], resolution=config.target_res
    )
    out_shape = (height, width)
    out_mask = np.empty(out_shape, dtype=np.uint8)
    
    reproject(
        source=mask.astype(np.uint8),
        destination=out_mask,
        src_transform=src_profile["transform"],
        src_crs=src_profile["crs"],
        dst_transform=transform,
        dst_crs=config.target_crs,
        resampling=Resampling.nearest
    )
    return out_mask.astype(bool)

@task(retries=2, retry_delay_seconds=10)
def generate_cloud_mask(qa_path: Path, sensor: str, config: MaskConfig) -> xr.DataArray:
    """Orchestrate mask generation, alignment, and xarray wrapping."""
    try:
        with rasterio.open(qa_path) as src:
            profile = src.profile.copy()
            profile["bounds"] = src.bounds  # profile dict has no bounds key
            if sensor.lower() == "sentinel2":
                mask = _parse_sentinel2_scl(qa_path)
            elif sensor.lower() in ("landsat8", "landsat9"):
                mask = _parse_landsat_qa(qa_path, config)
            else:
                raise ValueError(f"Unsupported sensor: {sensor}")
            
            aligned_mask = _align_and_resample(mask, profile, config)
            
            # Wrap in xarray with dask chunking for memory efficiency
            return xr.DataArray(
                da.from_array(aligned_mask, chunks=(1024, 1024)),
                dims=["y", "x"],
                coords={
                    "y": np.arange(aligned_mask.shape[0]) * config.target_res,
                    "x": np.arange(aligned_mask.shape[1]) * config.target_res
                }
            )
    except Exception as e:
        logger.error(f"Mask generation failed for {qa_path}: {e}")
        if config.fallback_skip:
            logger.warning("Fallback: Returning full-coverage mask (conservative)")
            return xr.DataArray(da.ones((1024, 1024), dtype=bool, chunks=(1024, 1024)), dims=["y", "x"])
        raise

@flow(name="mrva_cloud_masking_pipeline")
def run_cloud_masking_pipeline(
    qa_paths: list[Path],
    sensors: list[str],
    config: MaskConfig
) -> Dict[str, xr.DataArray]:
    """Execute multi-sensor masking workflow with audit logging."""
    masks = {}
    for qa_path, sensor in zip(qa_paths, sensors):
        logger.info(f"Processing {sensor} QA: {qa_path.name}")
        mask = generate_cloud_mask(qa_path, sensor, config)
        cloud_fraction = 1.0 - mask.mean().compute()
        
        if cloud_fraction > config.cloud_threshold:
            logger.warning(
                f"Cloud fraction {cloud_fraction:.2%} exceeds threshold {config.cloud_threshold:.2%} for {qa_path.name}"
            )
        
        masks[qa_path.stem] = mask
        
    logger.info(f"Pipeline complete. Generated {len(masks)} masks.")
    return masks

Debugging, Validation & Compliance Mapping #

Mask quality directly dictates the uncertainty bounds reported in carbon inventories. Validation must occur at three levels:

1. Statistical QA: Compute per-tile cloud fraction histograms. Sudden spikes in cloud_fraction > 0.30 often indicate misaligned QA bands or corrupted metadata. Implement automated flagging that routes suspect tiles to manual review queues.
2. Spectral Cross-Check: Apply the mask to reflectance bands and verify that NDVI/EVI distributions shift toward expected biome baselines. Persistent high NDVI in masked regions typically reveals cirrus leakage, which scatters near-infrared radiation.
3. Compliance Alignment: MRV frameworks such as Verra VM0042 and ISO 14064-2 require explicit documentation of data quality controls. The structured logging and threshold routing in the pipeline above generate immutable audit trails. Each mask’s cloud_fraction, confidence tier application, and fallback status must be persisted alongside the processed imagery to satisfy third-party verification audits.

When masking is inconsistent across acquisition dates, temporal compositing algorithms produce artificial step-changes. Properly masked inputs ensure that Temporal Aggregation for Land-Use Change routines compute genuine phenological or disturbance signals rather than atmospheric artifacts.

Pipeline Integration & Orchestration #

Cloud masking is rarely an isolated step. In production MRV systems, it feeds directly into change detection engines. Masked reflectance stacks are ingested by Deforestation Alert Generation Pipelines where spectral divergence thresholds trigger near-real-time notifications. To maintain throughput across continental-scale footprints, masking tasks should be distributed via dask.distributed or orchestrated through cloud-native schedulers.

For teams scaling ingestion across heterogeneous archives, STAC-compliant catalog traversal eliminates manual path resolution. The companion guide, Automating Sentinel-2 Cloud Masking with STAC and Rasterio, details how to bind STAC item metadata directly to the MaskConfig schema, enabling dynamic CRS inference and automated QA band discovery.

When integrating with multi-sensor fusion layers, ensure that mask alignment precedes spectral harmonization. Misaligned masks propagate through atmospheric correction steps, violating the assumption of pixel-level correspondence required for Advanced Multi-Sensor Data Fusion Techniques. By enforcing strict bitwise parsing, explicit CRS resampling, and structured fallback routing, engineering teams can guarantee that carbon accounting pipelines remain spectrally rigorous, computationally efficient, and audit-ready.