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:
- Which bands, and which are at 10 m, 20 m and 60 m?
- What window — the whole 110 km tile, or an area of interest?
- Which grid do the coarse bands get resampled onto, and how?
- What conversion turns the integers into reflectance?
- 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, profilebands, 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)}
Step-by-step solution
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 10 km window costs a few hundred kilobytes rather than the 209 MB of the full band. See How to read a COG from a URL without downloading the whole file.
Use a generous cloud filter. Scene-level cloud cover correlates only −0.654 with clarity over a small window, 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 per band as float32. Six bands is 1.4 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.bilinear5. Convert to reflectance — after checking the conversion
water = bands["red"][bands["scl"] == 6]
print(f"water red reflectance: {np.median(water):.3f}") # 0.008Water 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")48.9% usable
Code examples
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 S2A_30UVD_20250922_0_L2A (1017, 1087) usable 48.9%
S2C_30UVD_20250712_0_L2A (1017, 1087) usable 99.4%
S2B_30UVD_20250816_0_L2A (1017, 1087) usable 62.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 cubeFrozen({'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 rowsThe classification band is 20 m and single-byte, so one window is a few kilobytes. Ranking 120 scenes this way costs a fraction of loading even one full scene, and tells you which two or three are worth loading at all.
Explanation
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 one of these bands, a 512 × 512 window cost 1.59 MB in 5 requests against a 33 MB file. The saving grows with the file — the same window on the full 209 MB band costs the same 1.59 MB.
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.
Edge cases or notes
nodatais 0, which is a legal reflectance. Mask it explicitly rather than relying on the value.- Cast before arithmetic. These are
uint16arrays;dn - 1000wraps 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_DIRfor 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 — around 1.6 MB for a 512 × 512 window regardless of whether the file is 33 MB or 209 MB.
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.