How to Load Sentinel-2 Bands into Python as an Analysis-Ready Array

Problem statement

A Sentinel-2 scene is not an image file. It is a set of separate single-band GeoTIFFs at three different resolutions, each an unsigned integer array with a scale factor described somewhere else, plus a classification band you need before any of the numbers can be trusted.

"Load the scene" therefore means five decisions:

1. Which bands, and which are at 10 m, 20 m and 60 m?
2. What window — the whole 110 km tile, or an area of interest?
3. Which grid do the coarse bands get resampled onto, and how?
4. What conversion turns the integers into reflectance?
5. What mask marks the pixels that are not ground?

Get any of them wrong and the array still loads. It is just wrong.

Quick answer

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

AOI = (-4.16, 53.03, -4.00, 53.12)          # lon/lat


def load_bands(assets, names, aoi, reference="red", scale=1e-4):
    """Every band on the reference band's grid, as reflectance, plus the mask."""
    with rasterio.open(assets[reference]) as ds:
        bounds = transform_bounds("EPSG:4326", ds.crs, *aoi)
        window = rasterio.windows.from_bounds(
            *bounds, transform=ds.transform).round_offsets().round_lengths()
        shape = (int(window.height), int(window.width))
        profile = ds.profile | {
            "width": shape[1], "height": shape[0],
            "transform": ds.window_transform(window),
            "dtype": "float32", "count": 1, "nodata": np.nan,
        }

    out = {}
    for name in names:
        with rasterio.open(assets[name]) as ds:
            win = rasterio.windows.from_bounds(
                *transform_bounds("EPSG:4326", ds.crs, *aoi),
                transform=ds.transform).round_offsets().round_lengths()
            resampling = (Resampling.nearest if name == "scl"
                          else Resampling.bilinear)
            band = ds.read(1, window=win, out_shape=shape, resampling=resampling)
        out[name] = band if name == "scl" else band.astype("float32") * scale

    return out, profile

bands, profile = load_bands(assets, ["red", "nir", "swir16", "scl"], AOI)
print({k: v.shape for k, v in bands.items()})

{'red': (1017, 1087), 'nir': (1017, 1087),
 'swir16': (1017, 1087), 'scl': (1017, 1087)}

STAC hrefs, a window first, resampling on read

1. Get the asset hrefs from STAC rather than guessing paths

from pystac_client import Client

client = Client.open("https://earth-search.aws.element84.com/v1")
search = client.search(
    collections=["sentinel-2-l2a"],
    bbox=list(AOI),
    datetime="2025-06-01/2025-10-01",
    query={"eo:cloud_cover": {"lt": 60}},
)
items = sorted(search.items(), key=lambda i: i.properties["eo:cloud_cover"])
assets = {k: v.href for k, v in items[0].assets.items()}

The .tif hrefs are Cloud-Optimised GeoTIFFs, so a windowed read fetches only the tiles it needs — a small fraction of the full band, which for the red band of this scene is 200 MB on disk (Example 4). How to read a COG from a URL without downloading the whole file counts the bytes a window read actually transfers.

Use a generous cloud filter. Over 127 scenes of a 10 km window on Snowdon, scene-level cloud cover correlated only −0.656 with clarity over the window (measured here), so a tight filter discards scenes that were clear over your area.

2. Cut the window before anything else

Reading a full band gives a 10,980 × 10,980 array — 241 MB as the uint16 it is stored as, and 482 MB once cast to float32. Six bands is 2.9 GB before you have done anything.

Cut first. rasterio.windows.from_bounds(...).round_offsets().round_lengths() snaps the window to whole pixels; without the rounding you get fractional offsets and a shape that is off by one against another band.

3. Pick one reference grid and resample everything onto it

Sentinel-2 bands genuinely do not share a grid. For the same AOI:

red    (B04, 10 m): (1017, 1087)   origin (422200.0, 5886250.0)
swir16 (B11, 20 m): (509, 543)     origin (422200.0, 5886260.0)

The origins differ by 10 m, so even doubling the 20 m array is off by a row. Passing out_shape to ds.read() makes GDAL resample on read, using the window's own transform — correct by construction, and no extra memory.

Choose the reference band deliberately:

• 10 m reference keeps the detail of the visible and near-infrared bands and upsamples the shortwave-infrared. Best for mapping.
• 20 m reference downsamples the fine bands and invents nothing. Best for statistics and for any classifier using SWIR.

4. Resample class bands with nearest neighbour

Interpolating class codes is meaningless: halfway between class 4 (vegetation) and class 6 (water) is class 5 (bare soil).

resampling = Resampling.nearest if name == "scl" else Resampling.bilinear

5. Convert to reflectance — after checking the conversion

water = bands["red"][bands["scl"] == 6]
print(f"water red reflectance: {np.median(water):.3f}")

water red reflectance: 0.004

Water should be dark and positive in every band. If it comes out negative you have applied an offset that was already applied — the archive measured here declares offset: -0.1 in its metadata but ships data that already has it. See Radiometric levels explained.

6. Build the mask, then apply it once

UNUSABLE = (0, 1, 3, 8, 9, 10)
usable = ~np.isin(bands["scl"], UNUSABLE)
for name in bands:
    if name != "scl":
        bands[name] = np.where(usable, bands[name], np.nan)
print(f"{usable.mean():.1%} usable")

100.0% usable

A complete loader with metadata alongside

The data used on this page

Sentinel-2 L2A scenes of tile 30UVD from Earth Search, read straight from its Cloud-Optimised GeoTIFFs, for an area around Snowdon, 4.16–4.00° W, 53.03–53.12° N. The setup runs the search in step 1 once — June to September 2025, scene cloud cover under 60% — and hands the quick answer the band hrefs of the least cloudy scene as assets.

from pystac_client import Client

AOI = (-4.16, 53.03, -4.00, 53.12)
items = sorted(Client.open("https://earth-search.aws.element84.com/v1").search(
    collections=["sentinel-2-l2a"], bbox=list(AOI), datetime="2025-06-01/2025-10-01",
    query={"eo:cloud_cover": {"lt": 60}}).items(), key=lambda i: i.properties["eo:cloud_cover"])
assets = {k: v.href for k, v in items[0].assets.items()}
print(len(items), "scenes; least cloudy", items[0].id, f"({items[0].properties['eo:cloud_cover']:.1f}% cloud)")

27 scenes; least cloudy S2C_30UVD_20250712_0_L2A (0.1% cloud)

Example 1 — a complete loader with the metadata carried alongside

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

CLASS_BANDS = {"scl"}
UNUSABLE = (0, 1, 3, 8, 9, 10)


def load_scene(item, aoi, names=("blue", "green", "red", "nir"),
               reference="red", scale=1e-4, mask=True):
    """Load one STAC item as a dict of float32 arrays on one grid."""
    assets = {k: v.href for k, v in item.assets.items()}
    needed = list(dict.fromkeys(list(names) + (["scl"] if mask else [])))

    with rasterio.open(assets[reference]) as ds:
        crs = ds.crs
        bounds = transform_bounds("EPSG:4326", crs, *aoi)
        window = rasterio.windows.from_bounds(
            *bounds, transform=ds.transform).round_offsets().round_lengths()
        shape = (int(window.height), int(window.width))
        transform = ds.window_transform(window)

    bands = {}
    for name in needed:
        with rasterio.open(assets[name]) as ds:
            win = rasterio.windows.from_bounds(
                *transform_bounds("EPSG:4326", ds.crs, *aoi),
                transform=ds.transform).round_offsets().round_lengths()
            arr = ds.read(1, window=win, out_shape=shape,
                          resampling=Resampling.nearest if name in CLASS_BANDS
                          else Resampling.bilinear)
        bands[name] = arr if name in CLASS_BANDS else arr.astype("float32") * scale

    meta = {
        "id": item.id,
        "datetime": item.properties["datetime"],
        "crs": str(crs),
        "transform": list(transform)[:6],
        "shape": shape,
        "reference_band": reference,
        "scale": scale,
        "scene_cloud": item.properties.get("eo:cloud_cover"),
    }

    if mask:
        usable = ~np.isin(bands["scl"], list(UNUSABLE))
        for name in names:
            bands[name] = np.where(usable, bands[name], np.nan)
        meta["mask_drop_classes"] = list(UNUSABLE)
        meta["usable_fraction"] = round(float(usable.mean()), 4)

    print(f"  {item.id}  {shape}  usable {meta.get('usable_fraction', 1):.1%}")
    return bands, meta


scenes = [load_scene(item, AOI, names=("red", "nir")) for item in items[:12]]

  S2C_30UVD_20250712_0_L2A  (1017, 1087)  usable 100.0%
  S2B_30UVD_20250710_0_L2A  (1017, 1087)  usable 47.7%
  S2A_30UVD_20250922_0_L2A  (1017, 1087)  usable 48.8%
  S2A_30UVD_20250925_1_L2A  (1017, 1087)  usable 54.7%
  S2B_30UVD_20250925_0_L2A  (1017, 1087)  usable 50.6%
  S2A_30UVD_20250816_1_L2A  (1017, 1087)  usable 62.1%
  S2C_30UVD_20250702_0_L2A  (1017, 1087)  usable 30.9%
  S2C_30UVD_20250602_0_L2A  (1017, 1087)  usable 37.1%
  S2B_30UVD_20250816_0_L2A  (1017, 1087)  usable 63.3%
  S2A_30UVD_20250604_0_L2A  (1017, 1087)  usable 21.4%
  S2A_30UVD_20250912_0_L2A  (1017, 1087)  usable 4.6%
  S2B_30UVD_20250928_0_L2A  (1017, 1087)  usable 1.9%

The metadata dictionary is not decoration. When a composite is built from twenty of these, one assertion that every shape, transform and scale matches catches the misaligned scene before it becomes a smeared average.

Example 2 — loading a stack as an xarray cube instead

import numpy as np
import xarray as xr


def to_cube(scenes, names=("red", "nir")):
    """Turn a list of (bands, meta) into a labelled time-band-y-x cube."""
    first = scenes[0][1]
    times = np.array([np.datetime64(m["datetime"][:19]) for _, m in scenes])
    order = np.argsort(times)

    data = np.stack([
        np.stack([scenes[i][0][n] for n in names])
        for i in order
    ])                                             # (time, band, y, x)

    height, width = first["shape"]
    a, b, c, d, e, f = first["transform"]
    xs = c + (np.arange(width) + 0.5) * a
    ys = f + (np.arange(height) + 0.5) * e

    cube = xr.DataArray(
        data, dims=("time", "band", "y", "x"),
        coords={"time": times[order], "band": list(names), "y": ys, "x": xs},
        attrs={"crs": first["crs"], "scale": first["scale"]},
    )
    print(cube.sizes, f"{cube.nbytes / 1e6:.0f} MB")
    return cube


cube = to_cube(scenes)

Frozen({'time': 12, 'band': 2, 'y': 1017, 'x': 1087}) 106 MB

Once it is a cube, "median over time, ignoring masked pixels" is cube.median("time", skipna=True) and the alignment is guaranteed by the coordinates rather than by your own bookkeeping.

Example 3 — a fast pre-flight over many scenes

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


def rank_scenes(items, aoi, good=(4, 5, 6, 7, 11)):
    """Read only the 20 m classification band to decide what is worth loading."""
    rows = []
    for item in items:
        try:
            with rasterio.open(item.assets["scl"].href) as ds:
                win = rasterio.windows.from_bounds(
                    *transform_bounds("EPSG:4326", ds.crs, *aoi),
                    transform=ds.transform).round_offsets().round_lengths()
                scl = ds.read(1, window=win)
        except rasterio.errors.RasterioIOError:
            continue
        rows.append({"id": item.id, "date": item.properties["datetime"][:10],
                     "clear": round(float(np.isin(scl, list(good)).mean()), 4)})

    rows.sort(key=lambda r: -r["clear"])
    for r in rows[:5]:
        print(f"  {r['date']}  {r['clear']:6.1%}  {r['id']}")
    return rows


ranked = rank_scenes(items, AOI)
print(f"  {len(ranked)} scenes ranked from their classification band alone")

  2025-07-12   99.4%  S2C_30UVD_20250712_0_L2A
  2025-08-16   62.8%  S2B_30UVD_20250816_0_L2A
  2025-08-16   61.7%  S2A_30UVD_20250816_1_L2A
  2025-09-25   53.5%  S2A_30UVD_20250925_1_L2A
  2025-09-25   49.6%  S2B_30UVD_20250925_0_L2A
  27 scenes ranked from their classification band alone

The classification band is 20 m and single-byte, so one window is a few kilobytes. Ranking the 27 scenes here this way costs a fraction of loading even one full scene, and tells you which two or three are worth loading at all.

Example 4 — the sizes quoted on this page

import requests

item = items[0]
for name in ("red", "swir16", "coastal"):
    href = item.assets[name].href
    size = int(requests.head(href, timeout=60, allow_redirects=True).headers["Content-Length"])
    with rasterio.open(href) as ds:
        win = rasterio.windows.from_bounds(*transform_bounds("EPSG:4326", ds.crs, *AOI),
                                           transform=ds.transform).round_offsets().round_lengths()
        print(f"{name:8} {ds.res[0]:.0f} m: full band {ds.width:,} × {ds.height:,}, {size / 1e6:.0f} MB on disk, "
              f"{ds.width * ds.height * 2 / 1e6:.0f} MB as uint16, {ds.width * ds.height * 4 / 1e6:.0f} MB as float32; "
              f"AOI window {int(win.height)} × {int(win.width)}")

# bytes actually fetched for a window read, counted through rasterio's opener hook
import fsspec
from rasterio.windows import Window

fetched = {"bytes": 0}


class Counted:
    def __init__(self, f):
        self.f = f

    def read(self, n=-1):
        chunk = self.f.read(n)
        fetched["bytes"] += len(chunk)
        return chunk

    def __getattr__(self, name):
        return getattr(self.f, name)

    def __enter__(self):
        return self

    def __exit__(self, *exc):
        self.f.close()


http = fsspec.filesystem("http")
for label, pick in (("10 km AOI window", None), ("512 × 512 window", Window(0, 0, 512, 512))):
    fetched["bytes"] = 0
    with rasterio.open(item.assets["red"].href,
                       opener=lambda p, mode="rb": Counted(http.open(p, mode, block_size=2**16,
                                                                     cache_type="none"))) as ds:
        win = pick or rasterio.windows.from_bounds(*transform_bounds("EPSG:4326", ds.crs, *AOI),
                                                   transform=ds.transform)
        ds.read(1, window=win)
    print(f"{label}: {fetched['bytes'] / 1e6:.2f} MB fetched")

red      10 m: full band 10,980 × 10,980, 200 MB on disk, 241 MB as uint16, 482 MB as float32; AOI window 1017 × 1087
swir16   20 m: full band 5,490 × 5,490, 55 MB on disk, 60 MB as uint16, 121 MB as float32; AOI window 509 × 543
coastal  60 m: full band 1,830 × 1,830, 5 MB on disk, 7 MB as uint16, 13 MB as float32; AOI window 170 × 181
10 km AOI window: 5.48 MB fetched
512 × 512 window: 1.36 MB fetched

Why windowed reads are cheap

Why windowed reads are cheap

The assets are Cloud-Optimised GeoTIFFs: internally tiled, with the tile offsets in the header. A windowed read fetches the header and then only the tiles the window touches, using HTTP range requests.

Measured on the red band of the scene above (Example 4), a 512 × 512 window fetched 1.36 MB and the 10 km AOI window 5.48 MB, out of a 200 MB file. The saving grows with the file: the window costs what its tiles cost, whatever the size of the band.

This is why you should cut the window before doing anything else: the cut is not a filter applied after loading, it is what stops the loading.

Why out_shape beats resampling afterwards

Two ways to get a 20 m band onto a 10 m grid: read it and then resample the array, or pass out_shape and let GDAL do it during the read.

The second is better on three counts. It never materialises the intermediate array. It uses the dataset's own transform, so it cannot drift by a half pixel the way manual array scaling does. And where the file has overviews, GDAL can read a lower resolution level directly rather than reading full resolution and throwing pixels away.

Why the mask goes last but is decided first

Masking replaces values with NaN, which forces float32 and propagates through every subsequent operation. Doing it before resampling would smear NaNs across neighbouring pixels through the bilinear kernel, spreading the mask further than intended.

So the order is: read at native resolution, resample to one grid, scale, then mask. But the decision about which classes are unusable belongs at the top of the script as a named constant, because it is the single most consequential parameter in the whole pipeline.

Why 10 m is not automatically the right grid

Upsampling the shortwave-infrared band to 10 m produces an array with 10 m spacing and 20 m information. Every pixel is a real number, but neighbouring pixels are not independent.

For a map that is fine — you want the crispest visible bands. For statistics it is misleading: a classifier trained on upsampled SWIR sees four correlated pixels where there was one measurement, and any confidence interval computed from pixel counts is too narrow by a factor of two.

When in doubt, work at the coarsest resolution of the bands you are actually using.

Nodata zero, casting and rounded windows

• nodata is 0, which is a legal reflectance. Mask it explicitly rather than relying on the value.
• Cast before arithmetic. These are uint16 arrays; dn - 1000 wraps a DN of 1 to 64,537.
• Round the window. Without .round_offsets().round_lengths(), two bands can produce shapes that differ by one pixel.
• Nearest neighbour for class bands, always.
• Do not assume 2× for 20 m bands. The 10 m and 20 m windows over the same AOI here were 1087 and 543 columns — doubling gives 1086.
• Scene-level cloud cover is a weak filter. Use it loosely, then rank on the classification band.
• Set GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR for remote reads; it removed 4 of 10 HTTP requests in a four-window benchmark.
• A scene can straddle two UTM zones. Two items with the same date and different tile IDs are in different CRSs and cannot be stacked without reprojection.

Internal links

• Spectral bands explained: what satellite imagery actually measures — what each band is
• Radiometric levels explained — checking the scale before you trust it
• Cloud masking explained — choosing the class list
• How to resample satellite bands to a common grid — the grid decision in detail
• How to calculate NDVI from Sentinel-2 in Python — the first thing you do with the loaded array
• How to read a COG from a URL without downloading the whole file — why windowed reads are cheap
• How to turn a STAC search into an xarray data cube — the library version of Example 2
• Satellite bands have different shapes or do not align — when this goes wrong

FAQ

How do I load a Sentinel-2 scene in Python?

Read each band's COG with rasterio, cut to your area with a window, resample everything onto one reference grid with out_shape, scale to reflectance, and mask using the scene classification band.

Do I have to download the whole scene?

No. The assets are Cloud-Optimised GeoTIFFs, so a windowed read fetches only the tiles it needs — 1.36 MB for a 512 × 512 window of a 200 MB band, measured.

Why do the bands have different shapes?

Sentinel-2 has 10 m, 20 m and 60 m bands. For a 10 km window that is 1087, 543 and 181 columns, and the 20 m grid is not simply half the 10 m grid.

Which resampling method should I use?

Bilinear for continuous bands, nearest for class bands such as the scene classification layer. Interpolating class codes produces classes that were never observed.

Should I resample up to 10 m or down to 20 m?

Up for mapping, down for statistics. Upsampling gives 10 m spacing with 20 m information, which makes neighbouring pixels correlated and confidence intervals too narrow.

When should I apply the cloud mask?

After reading and resampling, before any index or statistic. Masking earlier smears NaN through the resampling kernel; masking later means cloud has already been averaged into your answer.

How do I stack several dates?

Load each scene onto the same reference grid, assert that the shapes and transforms match, then stack into an xarray cube with a time coordinate. Alignment by coordinate is safer than alignment by array index.