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

import rasterio

# the pixels, anonymously
with rasterio.open(sentinel_item.assets["red"].href) as ds:            # works
    print("Sentinel-2:", ds.width, "x", ds.height, ds.dtypes[0])
try:
    with rasterio.Env(AWS_NO_SIGN_REQUEST="YES"):
        rasterio.open(landsat_item.assets["red"].href)                 # AccessDenied
except rasterio.errors.RasterioIOError as e:
    print(f"Landsat: {type(e).__name__}: {e}")

Sentinel-2: 10980 x 10980 uint16
Landsat: RasterioIOError: AccessDenied: Anonymous users cannot invoke requests against Requester Pays buckets. Please authenticate.

Time span, feature size, two different archives

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 about 5.6 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:

for name, item in (("Sentinel-2 L2A", sentinel_item), ("Landsat C2 L2", landsat_item)):
    declared = item.assets["red"].extra_fields["raster:bands"][0]
    print(f"{name:16} scale {declared['scale']:<9g} offset {declared['offset']}")

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 between −0.090 and −0.098 in every band (Example 4), 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) & ~shadow.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.

One loader for both sensors

The data used on this page

Every Sentinel-2 L2A and Landsat Collection 2 Level-2 item that Earth Search lists for 2025 over a 10 km window on Snowdon, 4.16–4.00° W, 53.03–53.12° N. Sentinel-2 pixels are read from Earth Search's open COGs. Landsat pixels cannot be read anonymously from the USGS bucket (step 5), so they come from the Microsoft Planetary Computer copy of the same scenes, whose URLs its client signs for free.

import planetary_computer
from pystac_client import Client

AOI = (-4.16, 53.03, -4.00, 53.12)
earth = Client.open("https://earth-search.aws.element84.com/v1")
s2_items = list(earth.search(collections=["sentinel-2-l2a"], bbox=list(AOI),
                             datetime="2025-01-01/2025-12-31").items())
ls_items = list(earth.search(collections=["landsat-c2-l2"], bbox=list(AOI),
                             datetime="2025-01-01/2025-12-31").items())
cloud = lambda i: i.properties["eo:cloud_cover"]
summer = [i for i in s2_items if "2025-06" <= i.properties["datetime"][:7] <= "2025-09"]
sentinel_item = min(summer, key=cloud)                     # least cloudy, June-September
landsat_item = min(ls_items, key=cloud)                    # least cloudy of the year

mirror = Client.open("https://planetarycomputer.microsoft.com/api/stac/v1",
                     modifier=planetary_computer.sign_inplace)
landsat_mirror = lambda item: mirror.get_collection("landsat-c2-l2").get_item(item.id)
print(len(s2_items), "Sentinel-2 items;", len(ls_items), "Landsat items")
print("Sentinel-2:", sentinel_item.id, f"{cloud(sentinel_item):.2f}% cloud")
print("Landsat:   ", landsat_item.id, f"{cloud(landsat_item):.2f}% cloud")

201 Sentinel-2 items; 79 Landsat items
Sentinel-2: S2C_30UVD_20250712_0_L2A 0.06% cloud
Landsat:    LC09_L2SP_204023_20250410_02_T1 0.03% cloud

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

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

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
        shadow = (qa >> 4) & 1                      # "clear" means no cloud, not no shadow
        clear = (qa >> 6) & 1
        return clear.astype(bool) & ~shadow.astype(bool) & ~fill.astype(bool)
    raise ValueError(f"unknown sensor {sensor!r}")


def window_shape(ds, aoi):
    """The AOI's bounds in the dataset's CRS, and its shape on the dataset's grid."""
    bounds = transform_bounds("EPSG:4326", ds.crs, *aoi)
    win = rasterio.windows.from_bounds(*bounds, transform=ds.transform).round_offsets().round_lengths()
    return bounds, (int(win.height), int(win.width))


def read_window(href, bounds, shape, resampling):
    with rasterio.open(href) as ds:
        win = rasterio.windows.from_bounds(*bounds, transform=ds.transform)
        return ds.read(1, window=win, out_shape=shape, resampling=resampling)


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:
        bounds, shape = window_shape(ds, aoi)

    qa = read_window(assets[cfg["qa"]], bounds, shape, Resampling.nearest)

    bands = {}
    for name, asset in cfg["bands"].items():
        raw = read_window(assets[asset], bounds, shape, 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


s2_bands, s2_mask = load(sentinel_item, "sentinel-2-l2a", AOI)
ls_bands, ls_mask = load(landsat_mirror(landsat_item), "landsat-c2-l2", AOI)

  sentinel-2-l2a S2C_30UVD_20250712_0_L2A: (1017, 1087), 100.0% usable
  landsat-c2-l2 LC09_L2SP_204023_20250410_02_T1: (339, 362), 100.0% usable

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")


revisit_report(earth, AOI, 2025)

  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 difference is not hypothetical. Landsat 8 on 2 April 2025 and Sentinel-2 on 3 April, averaged to the same 30 m grid, gave mean NDVIs of 0.516 and 0.481 over 122,700 clear pixels of the window, with a correlation of only 0.812 (Example 4) — a day of change, a different view angle and sub-pixel misregistration as well as band widths. That one pair's fit is not a calibration, which is the point below.

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.

Example 4 — the measurements behind the explanation

import re
from collections import Counter

# 1. Where the looks come from
orbit = lambda i: re.search(r"_R(\d{3})_", i.properties["s2:product_uri"]).group(1)
print("Sentinel-2:", len(s2_items), "items on", len({i.properties["datetime"][:10] for i in s2_items}), "days;",
      "tiles", dict(Counter(i.properties["grid:code"] for i in s2_items)), "relative orbits", dict(Counter(map(orbit, s2_items))),
      "platforms", dict(Counter(i.properties["platform"] for i in s2_items)))
print("Landsat:   ", len(ls_items), "items on", len({i.properties["datetime"][:10] for i in ls_items}), "days;",
      "WRS-2 path/row", dict(Counter(i.id.split("_")[2] for i in ls_items)),
      "platforms", dict(Counter(i.properties["platform"] for i in ls_items)))

# 2. Is Sentinel-2's declared offset already applied? Look at water.
water = s2_mask & (read_window(sentinel_item.assets["scl"].href, *window_shape(
    rasterio.open(sentinel_item.assets["red"].href), AOI), Resampling.nearest) == 6)
print(f"Sentinel-2 water pixels: {water.sum():,}")
for band in ("blue", "green", "red", "nir"):
    dn_scaled = np.nanmedian(s2_bands[band][water])
    print(f"  {band:5} DN x 1e-4 = {dn_scaled:.4f}   with the declared -0.1 = {dn_scaled - 0.1:.4f}")

# 3. What np.isin does to a Landsat bit field, on a cloudy scene
cloudy = min(ls_items, key=lambda i: abs(cloud(i) - 50))
signed = landsat_mirror(cloudy)
with rasterio.open(signed.assets["qa_pixel"].href) as ds:
    bounds, shape = window_shape(ds, AOI)
qa = read_window(signed.assets["qa_pixel"].href, bounds, shape, Resampling.nearest)
print(f"{cloudy.id}, {cloud(cloudy):.1f}% scene cloud; most common QA_PIXEL values in the window:")
for value, n in Counter(qa.ravel().tolist()).most_common(4):
    flags = [name for bit, name in ((0, "fill"), (1, "dilated"), (2, "cirrus"), (3, "cloud"), (4, "shadow"),
                                     (5, "snow"), (6, "clear"), (7, "water")) if value >> bit & 1]
    print(f"  {value:6d} = {value:016b}  {n / qa.size:6.1%}  {flags}")
wrong = ~np.isin(qa, (0, 1, 3, 8, 9, 10))
right = usable_mask("landsat-c2-l2", qa)
cloudy_px = ((qa >> 3) & 1).astype(bool)
clear_only = ((qa >> 6) & 1).astype(bool) & ~(qa & 1).astype(bool)
print(f"  Sentinel-2 rule keeps {wrong.mean():.1%}, including {(wrong & cloudy_px).sum() / cloudy_px.sum():.1%} "
      f"of the cloud pixels; clear-bit alone keeps {clear_only.mean():.1%}; clear without shadow {right.mean():.1%}")

# 4. The same ground one day apart: Landsat 2 April, Sentinel-2 3 April
ls_day = next(i for i in ls_items if i.properties["datetime"][:10] == "2025-04-02")
s2_day = min((i for i in s2_items if i.properties["datetime"][:10] == "2025-04-03"), key=cloud)
lb, lm = load(landsat_mirror(ls_day), "landsat-c2-l2", AOI)
with rasterio.open(landsat_mirror(ls_day).assets["red"].href) as ds:
    bounds, shape = window_shape(ds, AOI)
s2_on_30m = {b: read_window(s2_day.assets[b].href, bounds, shape, Resampling.average).astype("float32") * 1e-4
             for b in ("red", "nir")}
ndvi_ls = (lb["nir"] - lb["red"]) / (lb["nir"] + lb["red"])
ndvi_s2 = (s2_on_30m["nir"] - s2_on_30m["red"]) / (s2_on_30m["nir"] + s2_on_30m["red"])
ok = lm & np.isfinite(ndvi_ls) & np.isfinite(ndvi_s2)
b, a = np.polyfit(ndvi_ls[ok], ndvi_s2[ok], 1)
print(f"NDVI on {ok.sum():,} shared 30 m pixels: Landsat mean {ndvi_ls[ok].mean():.3f}, Sentinel-2 mean "
      f"{ndvi_s2[ok].mean():.3f}, r = {np.corrcoef(ndvi_ls[ok], ndvi_s2[ok])[0, 1]:.3f}, "
      f"fit S2 = {a:.3f} + {b:.3f} x Landsat")

Sentinel-2: 201 items on 172 days; tiles {'MGRS-30UVD': 201} relative orbits {'037': 100, '080': 101} platforms {'sentinel-2a': 62, 'sentinel-2c': 67, 'sentinel-2b': 72}
Landsat:    79 items on 79 days; WRS-2 path/row {'204023': 39, '205023': 40} platforms {'landsat-8': 40, 'landsat-9': 39}
Sentinel-2 water pixels: 23,380
  blue  DN x 1e-4 = 0.0096   with the declared -0.1 = -0.0904
  green DN x 1e-4 = 0.0092   with the declared -0.1 = -0.0908
  red   DN x 1e-4 = 0.0039   with the declared -0.1 = -0.0961
  nir   DN x 1e-4 = 0.0018   with the declared -0.1 = -0.0982
LC09_L2SP_204023_20250528_02_T1, 47.6% scene cloud; most common QA_PIXEL values in the window:
   22280 = 0101011100001000   87.1%  ['cloud']
   23888 = 0101110101010000    3.8%  ['shadow', 'clear']
   23826 = 0101110100010010    3.7%  ['dilated', 'shadow']
   21824 = 0101010101000000    2.0%  ['clear']
  Sentinel-2 rule keeps 100.0%, including 100.0% of the cloud pixels; clear-bit alone keeps 6.4%; clear without shadow 2.2%
  landsat-c2-l2 LC08_L2SP_204023_20250402_02_T1: (339, 362), 100.0% usable
NDVI on 122,700 shared 30 m pixels: Landsat mean 0.516, Sentinel-2 mean 0.481, r = 0.812, fit S2 = -0.229 + 1.377 x Landsat

Revisit and resolution in practice

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. A two-satellite Sentinel-2 constellation repeats every 5 days. That predicts about 1.6:1; the measurement was 201 to 79, about 2.5:1.

Swath overlap is not the difference, because both missions have it here. The window lies in two Sentinel-2 relative orbits (037 and 080, with 100 and 101 items) and in two Landsat paths (204 and 205, with 39 and 40), so each sees it twice per repeat cycle. The difference is a third satellite: in 2025 Sentinel-2A, 2B and 2C were all imaging, with 62, 72 and 67 items. Three satellites on a 10-day orbit is one pass per 3.3 days per relative orbit, which predicts 2.4:1 against Landsat — and on 29 days two of them imaged the window from different orbits, so the 201 items fell on 172 days.

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 — nor is 22280, cloud. On a Landsat 9 scene of 28 May 2025 with 47.6% scene cloud, the Sentinel-2 rule kept 100.0% of the window, including every cloud pixel, with no error raised. The bit rule kept 2.2%.

The bit field has its own trap. Bit 6, "clear", means no cloud and no dilated cloud — not no shadow. 23888 is a cloud-shadow pixel with the clear bit set, and it was 3.8% of that window, so testing the clear bit alone kept 6.4%. Test the shadow bit too, as usable_mask does, or the Landsat mask keeps shadow that the Sentinel-2 one removes.

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.

Landsat 7's failure and band naming

• 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
• Global Forest Change explained — the Landsat-based forest loss record
• Landsat thermal data explained — the thermal record, sensors and drift

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 in 2025, measured here as 201 against 79 items over one 10 km window at 53° N, with three Sentinel-2 satellites imaging.

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.