Sentinel-2 or Landsat? Choosing a Satellite Imagery Source #

Problem statement #

Both are free, both are optical multispectral, both have red, near-infrared and shortwave-infrared bands, and both are indexed in the same STAC catalogues. Choosing between them looks like a preference.

It is not. The two differ in resolution, revisit rate, band scaling, mask format, archive depth and — the one that surprises people mid-project — how you are allowed to read the pixels.

Measured over one 10 km window in Snowdonia for calendar year 2025:

Sentinel-2 L2A items    201
Landsat C2 L2 items      79

Two and a half times as many looks, at three times the spatial resolution. That settles a lot of cases before any other consideration.

Quick answer #

# the pixels, anonymously
rasterio.open(sentinel_item.assets["red"].href)     # works

with rasterio.Env(AWS_NO_SIGN_REQUEST="YES"):
    rasterio.open(landsat_item.assets["red"].href)  # AccessDenied

RasterioIOError: AccessDenied: Anonymous users cannot invoke requests
against Requester Pays buckets.

Step-by-step solution #

1. Start with the question's time span #

This is usually decisive and it needs no analysis.

• Anything before 2015 is Landsat. Sentinel-2A did not launch until June 2015 and full two-satellite coverage came later. Landsat's usable multispectral record runs back to the 1980s.
• Recent monitoring favours Sentinel-2, on both resolution and revisit — 201 opportunities against 79 in the window measured here.
• A long series ending recently needs both, harmonised, which is a real piece of work: different band widths, different resolutions, different atmospheric correction.

2. Check the resolution against your features #

Sentinel-2's 10 m red and near-infrared give 100 m² pixels; Landsat's 30 m give 900 m². A 0.5 ha field is 50 Sentinel-2 pixels and 5 Landsat pixels.

Below about 1 ha, Landsat gives you a handful of pixels, most of them mixed with the boundary. Above about 10 ha the difference stops mattering for a mean, and Landsat's longer archive starts to win.

3. Handle the two scalings separately #

Both declare a scale and an offset, and they are different numbers:

Sentinel-2 L2A   scale 0.0001    offset -0.1
Landsat C2 L2    scale 2.75e-05  offset -0.2

Neither should be applied without checking. For the Sentinel-2 COGs measured here, applying the declared offset gives water a reflectance of −0.09 in every band, which means it was already applied — see Radiometric levels explained.

Landsat Collection 2 surface reflectance genuinely does need DN * 2.75e-5 - 0.2. Run the same water check on whichever copy you have, because the answer is a property of the archive, not of the mission.

4. Write two different mask functions #

The mask formats are not variations on a theme:

# Sentinel-2: SCL holds a class code per pixel
usable = ~np.isin(scl, (0, 1, 3, 8, 9, 10))

# Landsat: QA_PIXEL packs flags into bits
fill        = (qa >> 0) & 1
dilated     = (qa >> 1) & 1
cirrus      = (qa >> 2) & 1
cloud       = (qa >> 3) & 1
shadow      = (qa >> 4) & 1
snow        = (qa >> 5) & 1
clear       = (qa >> 6) & 1
water       = (qa >> 7) & 1
usable = clear.astype(bool) & ~fill.astype(bool)

np.isin on a bit field returns nonsense that looks like a mask. This is the single most common bug when a Sentinel-2 pipeline is pointed at Landsat.

5. Check the access model before you plan the compute #

The Sentinel-2 L2A COGs on AWS are openly readable — a windowed read from a laptop needs no credentials. The USGS Landsat Collection 2 bucket is requester-pays, and anonymous reads fail outright, as shown above.

That is not a blocker; it is a configuration and a cost line. Set credentials and AWS_REQUEST_PAYER=requester, or use a mirror such as a commercial or institutional STAC that hosts open copies. But discover it now rather than when the batch job fails at three in the morning.

6. Use Landsat when you need the thermal band #

Landsat carries a thermal infrared instrument; Sentinel-2 does not. Anything about land surface temperature, evapotranspiration, urban heat or thermal anomalies is Landsat, full stop — regardless of what the resolution comparison says.

Code examples #

Example 1 — one loader for both, with the differences isolated #

import numpy as np
import rasterio
from rasterio.enums import Resampling

SENSORS = {
    "sentinel-2-l2a": {
        "bands": {"blue": "blue", "green": "green", "red": "red", "nir": "nir",
                  "swir16": "swir16"},
        "qa": "scl", "scale": 1e-4, "offset": 0.0,     # verified against water
        "reference": "red",
    },
    "landsat-c2-l2": {
        "bands": {"blue": "blue", "green": "green", "red": "red",
                  "nir": "nir08", "swir16": "swir16"},
        "qa": "qa_pixel", "scale": 2.75e-5, "offset": -0.2,
        "reference": "red",
    },
}


def usable_mask(sensor, qa):
    """The one function that genuinely differs between the two."""
    if sensor == "sentinel-2-l2a":
        return ~np.isin(qa, (0, 1, 3, 8, 9, 10))
    if sensor == "landsat-c2-l2":
        fill = (qa >> 0) & 1
        clear = (qa >> 6) & 1
        return clear.astype(bool) & ~fill.astype(bool)
    raise ValueError(f"unknown sensor {sensor!r}")


def load(item, sensor, aoi):
    cfg = SENSORS[sensor]
    assets = {k: v.href for k, v in item.assets.items()}

    with rasterio.open(assets[cfg["bands"][cfg["reference"]]]) as ds:
        shape = window_shape(ds, aoi)

    with rasterio.open(assets[cfg["qa"]]) as ds:
        qa = ds.read(1, out_shape=shape, resampling=Resampling.nearest)

    bands = {}
    for name, asset in cfg["bands"].items():
        with rasterio.open(assets[asset]) as ds:
            raw = ds.read(1, out_shape=shape, resampling=Resampling.bilinear)
        bands[name] = raw.astype("float32") * cfg["scale"] + cfg["offset"]

    mask = usable_mask(sensor, qa)
    for name in bands:
        bands[name] = np.where(mask, bands[name], np.nan)

    print(f"  {sensor} {item.id}: {shape}, {mask.mean():.1%} usable")
    return bands, mask

Putting the scale, offset, band aliases and QA rule in one dictionary per sensor is what keeps a two-sensor pipeline honest. The moment those live inline in the processing code, one of them will be applied to the other sensor.

Example 2 — checking the revisit you will actually get #

from collections import Counter


def revisit_report(client, aoi, year, collections=("sentinel-2-l2a",
                                                   "landsat-c2-l2")):
    """Items, and items you could actually use, per collection."""
    for collection in collections:
        items = list(client.search(collections=[collection], bbox=list(aoi),
                                   datetime=f"{year}-01-01/{year}-12-31").items())
        clouds = [i.properties.get("eo:cloud_cover", 100) for i in items]
        months = Counter(i.properties["datetime"][5:7] for i in items)
        under20 = sum(1 for c in clouds if c < 20)
        print(f"  {collection:18} {len(items):4d} items, "
              f"{under20:3d} under 20% scene cloud, "
              f"{len(months)} months represented")

  sentinel-2-l2a      201 items,  21 under 20% scene cloud, 12 months represented
  landsat-c2-l2        79 items,  11 under 20% scene cloud, 11 months represented

Item count is not the whole story. Run this for your own area before committing — in a persistently cloudy region the number that matters is not overpasses but usable overpasses, and that ratio differs by climate rather than by satellite.

Example 3 — harmonising the two for one time series #

import numpy as np


def harmonise_ndvi(value, sensor, coefficients=None):
    """Adjust one sensor's NDVI onto the other's scale.

    Band widths and centre wavelengths differ, so the same ground gives
    slightly different NDVI. Published cross-calibrations are linear.
    """
    coefficients = coefficients or {
        # NDVI_s2_equivalent = a + b * NDVI_landsat  (illustrative form only —
        # take the constants from a published calibration for your bands)
        "landsat-c2-l2": {"a": 0.0, "b": 1.0},
        "sentinel-2-l2a": {"a": 0.0, "b": 1.0},
    }
    c = coefficients[sensor]
    return c["a"] + c["b"] * np.asarray(value)

The important part of this function is not the arithmetic, which is trivial. It is that the sensor is an explicit argument, so a merged series cannot silently mix two scales. Take the constants from a published cross-calibration for the exact band pair — do not fit them on your own cloudy overlap of a handful of dates.

Explanation #

Why the revisit difference is larger than the orbit difference #

Landsat 8 and 9 each repeat every 16 days, staggered to give 8 days combined. Sentinel-2A, 2B and 2C give 5 days at the equator and better at high latitudes where adjacent swaths overlap.

That predicts roughly a 2:1 ratio; the measurement was 201 to 79, about 2.5:1. The extra comes from swath overlap at 53° north, where a single location falls in more than one Sentinel-2 relative orbit.

In a cloudy climate the revisit rate matters more than anything else, because usable observations are a small fraction of overpasses. Two and a half times the looks is two and a half times the chance of catching a clear day.

Why 10 m is not always three times better #

Sentinel-2's 10 m applies to four bands: blue, green, red and the broad near-infrared. The red-edge bands and both shortwave-infrared bands are 20 m.

So an analysis using SWIR — burn severity, moisture, snow, built-up indices — works at 20 m on Sentinel-2 against 30 m on Landsat. That is a 2.25× difference in pixel area rather than 9×, and Landsat's longer archive often outweighs it.

Why the mask difference causes silent bugs #

An SCL band holds small integers 0–11. A QA_PIXEL band holds packed flags, with values in the tens of thousands.

Point Sentinel-2 code at Landsat and np.isin(qa, (0, 1, 3, 8, 9, 10)) will match almost nothing, because a QA_PIXEL value of 21824 (a clear land pixel) is not in the list. The result is a mask that keeps everything, including all the cloud, with no error raised.

Symptom: a suspiciously high usable fraction and a composite that looks like the cloud composite in How to build a cloud-free composite. Always print the usable fraction after masking.

Why requester-pays changes architecture #

An openly readable COG means a windowed read from anywhere with no account, which is what makes the "read only the tiles you need" pattern so effective for Sentinel-2.

Requester-pays means every read is billed to your AWS account, including the reads that turn out to be cloud. That pushes the design towards pre-filtering on metadata, colocating compute in the same region, and caching aggressively — all things worth deciding at the start rather than discovering from a bill.

Edge cases or notes #

• Landsat 7 has the scan-line corrector failure (post-2003): usable, striped, and needs gap handling.
• Landsat 8 and 9 are interchangeable for most purposes; Landsat 5 and 7 need different band indices — the same physical band has a different number.
• Sentinel-2's nir (B08) and nir08 (B8A) are different bands at 10 m and 20 m. Landsat's nir08 is closer to B8A than to B08.
• Both declare a scale and offset; verify both against a dark target.
• Landsat thermal has no Sentinel-2 equivalent.
• Sentinel-2 tiles overlap, so one date can produce two items for one location.
• Landsat's WRS-2 path/row and Sentinel-2's MGRS tiles are different grids; a mosaic needs reprojection either way.
• Harmonised Landsat–Sentinel products exist and are usually a better idea than harmonising them yourself.

Internal links #

• Spectral bands explained — why band widths make the two disagree
• Radiometric levels explained — verifying each archive's scale and offset
• Cloud masking explained — the SCL side of the mask difference
• How to mask clouds in Sentinel-2 imagery with Python — the implementation to adapt
• How to search and download satellite imagery with STAC in Python — querying both catalogues
• How to load Sentinel-2 bands into Python as an analysis-ready array — the loader to generalise
• How to read spatial data from S3 and other object storage — configuring requester-pays access
• STAC explained: how satellite imagery catalogues work — where all this metadata lives

FAQ #

Should I use Sentinel-2 or Landsat? #

Sentinel-2 for anything since 2015 needing resolution or frequency — it gave 201 items against Landsat's 79 over the same area in 2025. Landsat for anything before 2015, anything thermal, and anything at scales where 30 m is enough.

How much more often does Sentinel-2 pass over? #

About 2.5 times as often at mid-latitudes, measured here as 201 against 79 items over one 10 km window in a year.

Do both need a scale and offset applied? #

Both declare one. Landsat Collection 2 genuinely needs DN * 2.75e-5 - 0.2. Whether Sentinel-2 needs its declared −0.1 depends on the archive — check against water before applying either.

Why does my Landsat cloud mask keep everything? #

Because QA_PIXEL is a bit field, not class codes. np.isin matches nothing; use bitwise shifts and tests instead.

Why can I not read Landsat pixels without credentials? #

The USGS Landsat Collection 2 bucket is requester-pays. Anonymous reads return AccessDenied. Set AWS credentials and AWS_REQUEST_PAYER=requester, or use an open mirror.

Can I combine Landsat and Sentinel-2 in one time series? #

Only with a cross-calibration. Band widths and centre wavelengths differ, so the same ground gives slightly different index values. Prefer a published harmonised product.

Does Sentinel-2 have a thermal band? #

No. Anything about surface temperature needs Landsat.