Open Earth MonitorLandsatMGM Pipeline
Open Earth Monitor

LandsatMGM Data Processing Pipeline

Monthly, cloud-free Landsat surface reflectance for 1997–2025, built with temporal aggregation and SoftImpute gap filling.

Open Earth MonitorOpenGeoHub

Use ← →, swipe, or the buttons below to move between slides. Click any figure to enlarge it.

Each pass images a swath about 185 km wide.

Landsat

The longest continuous record of Earth's land surface from space.

1972Continuous Earth observation archive
30 mSpatial resolution, with consistent data since 1982
FreeUnrestricted, public data access since 2008
GlobalMonitoring forests, cities, agriculture, and water
Observations are organized into 1° × 1° tiles.

GLAD Landsat ARD

Analysis-ready data covering 1997–2025, the input to this pipeline.

Sensor harmonization and 16-day composites

Minimizes sensor effects across Landsat 5 TM, 7 ETM+, and 8 OLI/TIRS by cross-calibrating surface reflectance against MODIS.

Consistent spectral coverage

23 observations per year across all seven spectral bands.

19,132 standardized tiles

1° × 1° latitude/longitude tiles of 4004 × 4004 pixels, stored as int16. About 600 TB in total.

Artifacts in the raw composites

Blue band, tile 003W_45N. Clouds, sensor stripes, and missing swaths are common in single 16-day composites. Click an image to enlarge it.

Blue band composite with cloud and stripe artifacts
Blue band composite with artifacts
Blue band composite with artifacts

Version 1: cloud-free composites

The first version worked, but had four limitations that motivated a redesign.

Temporal limitations

60 days

Bi-monthly steps are too coarse for many applications.

Data type

8 bit

Values were compressed and smoothed.

Limited quality

Artifacts

Remaining cloud and stripe artifacts.

No quality mask

Uncertainty

No mask showing which values were observed or filled.

Aim

What LandsatMGM had to achieve.

1

Increase temporal frequency

2

Develop a more advanced gap-filling approach

3

Provide better-quality, ready-to-use data and side products

4

Build a workflow that can run at scale on OpenGeoHub infrastructure

Potential reconstruction solutions

Deep learning approaches were considered first.

Spatio-temporal

Spatio-temporal and super-resolution models

Deep learning models (e.g., CFC, super-resolution) that use spatial context and temporal patterns together.

Generative AI

GANs and ER-CGN models

Synthesize missing data with generative adversarial networks or ER-CGN frameworks for realistic texture.

Hardware and feasibility limit

Most high-end generative and spatio-temporal architectures are too slow at this scale. Testing showed they are not feasible on our current CPU infrastructure.

Super-resolution example

Selected area in Brazil, September.

Compare

MODIS: blue band, 500 m. Model output: blue band, 30 m. Google satellite imagery as visual reference.

MODIS 500 m Model output 30 m
↔

Multi-step gap-filling pipeline

A lighter, CPU-friendly alternative in two stages.

Stage 1, temporal

Temporal aggregation

Aggregates observations over short time windows to resolve most cloud gaps and establish a strong baseline.

Stage 2, context and PCA

Neighbor context and SoftImpute

Initial values are estimated from temporal neighbors. SoftImpute, an iterative low-rank matrix completion, then reconstructs the persistent gaps.

LandsatMGM: monthly cloud-free

LandsatMGM extends the consistent, gap-free archive to 1997–2025.

High-precision pipeline

Keeps the original 16-bit data, with multi-step gap filling.

Temporal aggregation and SoftImpute

SoftImpute matrix completion fills the gaps left after aggregation.

Long-term products and quality masks

29-year quartiles and medians, PCA components, and quality masks showing the gap-filling level of every pixel.

Pipeline overview

Step through the pipeline, from GLAD Landsat ARD to the LandsatMGM data cube.

Click an image in the diagram to enlarge it

Stage 1: temporal aggregation

Cloud, haze, and shadow are removed with the GLAD quality mask. Each monthly value is then a weighted mean of up to seven 16-day observations. Click the observations to mark them clear or cloudy.

Rule

  • At least 1 clear observation in the target month M
  • At least 2 clear observations in total
  • Weights: 1 for month M, 0.25 for the neighboring months

Months have 1 to 3 composites. A missing M₃ simply counts as unavailable.

Try a scenario

    Reflectance values are illustrative.

    Real example: aggregation vs gap filling

    NIR band, February 2007, 10.00°S–9.49°S, 58.00°W–57.39°W.

    After stage 1, pixels that fail the rule remain gaps (white). Stage 2 fills them with SoftImpute.

    Drag across the image to compare.

    NIR after gap filling, February 2007 NIR after temporal aggregation only, with gaps Aggregation only After gap filling
    ↔

    Initial values for the gaps

    Before SoftImpute runs, every gap needs a starting value. It comes from the same pixel, band, and month in the surrounding years. Click the years to mark them as observed or missing.

    Same month, nearby years

    Take the 5 closest years that have a valid value for the same month and band.

    Closer years count more

    The weight halves with each year of distance from the target year. The initial value is the weighted mean.

    Fallback: long-term median

    If that month has no valid value in any year, the gap-filled long-term median is used instead.

    Then SoftImpute

    The initial value is only a starting point. SoftImpute re-estimates it from the other months of the same year.

    Target: February, NIR band

    observedmissingused, with weighttarget gap

    Reflectance values are illustrative.

    SoftImpute: from samples to saved model

    Step through sample collection, model calibration, the saved model, and gap filling.

    Click an image in the diagram to enlarge it

    Stage 2: training SoftImpute

    Each pixel is a vector of 99 features. SoftImpute repeatedly fits a low-rank model and shrinks every component’s weight. Weak components fall below a fixed threshold and are removed. Move the slider through the iterations.

    StartIteration
    99Components kept
    –Change in filled values

    Component strength

    Blue: above the fixed threshold. Orange: below it, removed at the next iteration. Grey: how much each weight has shrunk so far.

    One pixel, blue band, 12 months

    Filled dots are observed. Open dots are gaps being reconstructed.

    Illustration of the procedure. The final model after 10 iterations keeps 32 components; intermediate counts and values are schematic.

    Input features

    99 features per sample. Hover over or tap a cell, or click a row in the table to highlight a group.

     
    #FeaturesDescription
    0–8384 spectral7 bands × 12 months (red, NIR, blue, green, SWIR1, SWIR2, thermal)
    84–9512 geo-temperatureMonthly temperature
    96SlopeTerrain slope in degrees
    97TWITopographic wetness index
    98HillshadeHillshade

    Optimization

    Iterations: 10

    Error levels off after about 10 iterations.

    Sigma: 0.01

    Shrinkage threshold chosen from the validation curves.

    Multi-band validation

    RMSE assessed across 7 bands: red, NIR, blue, green, SWIR1, SWIR2, thermal.

    Training and validation curves

    Validation

    80/20 train/test split, plus two extra sample sets from new tiles: pixels with land-cover change and stable pixels. Standardized RMSE per band, compared with the bi-monthly version 1.

    Lower is better

    Green: lower RMSE than bi-monthly. Red: higher.

    Result: NIR band, 2021

    From the global mosaic down to individual fields: zooming in on Levy County, Florida (tile 029N_082W), at the full 30 m resolution.

    Global

    Result: NDVI in northern Italy

    LandsatMGM NDVI for Lombardy, every month from January 1997 to December 2025. Every month is gap-free, so seasonal cycles in the Po Valley farmland and the Alps show without interruption.

    1997-01

    Speed

    Low NDVIHigh NDVI

    Result: USA

    Tile 108W_44N. Around Lovell and Byron. Version 1 leaves most of this area without data, while LandsatMGM is complete.

    Drag across the image to compare

    Version 1: bi-monthly, 8-bit. LandsatMGM: monthly, 16-bit, gap-filled. Google satellite imagery as visual reference.

    Version 1 LandsatMGM
    ↔

    Result: Peru

    Tile 074W_12S. Mountain valleys around Campo Armiño. Version 1 shows stripe artifacts and bright patches; LandsatMGM has far fewer.

    Drag across the image to compare

    Version 1: bi-monthly, 8-bit. LandsatMGM: monthly, 16-bit, gap-filled. Google satellite imagery as visual reference.

    Version 1 LandsatMGM
    ↔

    Result: Russia

    Tile 097E_70N. Tundra with rivers and lakes. Version 1 shows Landsat 7 scan-line stripes and snow; LandsatMGM is stripe-free.

    Drag across the image to compare

    Version 1: bi-monthly, 8-bit. LandsatMGM: monthly, 16-bit, gap-filled. Google satellite imagery as visual reference.

    Version 1 LandsatMGM
    ↔

    Result: deforestation in Brazil

    NDVI from LandsatMGM and its quality mask, tile 063W_11S, around Seringueiras (63.20°W–63.00°W, 12.00°S–11.80°S). One frame per year, 1997–2025.

    View

    Layer

    1997
    19972025
    00.51
    NDVI
    Scroll or pinch to zoom, drag to pan, double-click to zoom in

    NDVI time series

    LandsatMGM NDVI compared with HLS and the bi-monthly version.

    4 land covers

    Short vegetation, tall trees, cropland, built-up

    2014–2024

    Continuous evaluation period

    Baselines

    HLS L30, HLS S30, and bi-monthly ARD2

    NDVI time series for four land covers compared with HLS and bi-monthly data

    Representation and data storage

    How the archive is stored so it stays small, fast to read, and simple to use for training.

    One COG per tile, without overviews

    Each 1° × 1° tile is a Cloud Optimized GeoTIFF with no overview levels. Leaving out the overviews keeps files smaller and as simple as possible, which suits training, where whole tiles are read at full resolution.

    int16 instead of float32

    Values are stored as 16-bit integers, which halves the size of every pixel.

    ZSTD compression

    We tested several compression algorithms, and ZSTD gave the best balance for our data. Stronger settings, such as higher ZSTD levels or other algorithms, saved only a little more space but made reading and writing slower, so we avoided them.

    How much data type and compression matter

    7 spectral bands, 4004 × 4004 pixels per tile.

    Data type

    Amount of data

    –uncompressed
    –as float32
    –stored: int16 + ZSTD
    int16
    float32
    int16 + ZSTD

    Each stored file (1 band, 1 month, int16 with ZSTD) is about 17 MB, compared with 32 MB uncompressed and 64 MB as float32. The full archive is about 745 TB.

    Choosing the compression

    The same data compressed with ZSTD and DEFLATE at every level. Higher levels shrink files a little, but time grows much faster. Hover over or tap a point.

    Show all values

    One mosaic, no duplicate data

    A VRT joins all tiles into one global layer without copying them. Its overviews (.vrt.ovr) are stored as a global COG, so a viewer reads the overview when zoomed out and the original tiles when zoomed in.

    GlobalField

    What gets stored

    Bar length shows the total amount of data stored.

    Tiles (COG)Overviews

    Overview share assumes each level halves the resolution, as usual for overviews: ¼ + ¹⁄₁₆ + … ≈ ⅓.

    STAC catalog and OpenLandMap

    A STAC catalog describes every layer, so the data is ready to ingest into the OpenLandMap platform, where anyone can explore it. Click a step to read more.