How to Search and Download Satellite Imagery with STAC in Python
Problem statement
You need a summer NDVI map for one city. The old route was a portal account, a search form, a queued order and a 1.2 GB ZIP for a scene 200 times larger than your study area.
With STAC the whole thing is a script — but only if you get four things right, and each of them fails silently:
- The wrong product level. L1C looks identical to L2A and produces a time series that measures the atmosphere.
- No cloud masking. A scene at 0.06% cloud can have its only cloud over your city.
- The whole scene downloaded. Reading one window instead is a hundred times faster and needs no disk.
- The wrong scale factor. Reflectance is stored as integers, and the conversion is not always what the metadata claims.
That last one is the nastiest. Get it wrong and NDVI comes out above 1.0 — physically impossible, arithmetically clean, and easy to miss if you never check the range.
Quick answer
Search, mask, window-read, and check the result is physically possible:
import os
os.environ["GDAL_DISABLE_READDIR_ON_OPEN"] = "EMPTY_DIR"
os.environ["AWS_NO_SIGN_REQUEST"] = "YES"
import numpy as np
from pystac_client import Client
catalog = Client.open("https://earth-search.aws.element84.com/v1")
item = next(catalog.search(
collections=["sentinel-2-l2a"],
bbox=[-2.30, 53.44, -2.18, 53.52],
datetime="2026-06-01/2026-08-01",
query={"eo:cloud_cover": {"lt": 20}},
).items())
print(item.id, item.datetime.date(), round(item.properties["eo:cloud_cover"], 2))S2B_30UWE_20260712_0_L2A 2026-07-12 0.06
Then read only your window, and assert the NDVI is in range:
assert -1.0 <= np.nanmin(ndvi) and np.nanmax(ndvi) <= 1.0, "scale factor is wrong"
Step-by-step solution
1. Search with the constraints that eliminate the most scenes
search = catalog.search(
collections=["sentinel-2-l2a"],
bbox=[-2.30, 53.44, -2.18, 53.52],
datetime="2026-06-01/2026-08-01",
query={"eo:cloud_cover": {"lt": 20}},
)
print(search.matched(), "scenes")5 scenes
Twenty-three scenes intersect that box in that window; five are under 20% cloud. Filtering on the server costs nothing and saves you from scoring twenty-three footprints.
2. Configure GDAL before opening anything
Two environment variables change remote reads from "slow" to "fast":
os.environ["GDAL_DISABLE_READDIR_ON_OPEN"] = "EMPTY_DIR"
os.environ["AWS_NO_SIGN_REQUEST"] = "YES"The first stops GDAL listing the entire S3 prefix looking for sidecar files — on a bucket with millions of objects that can take longer than the read. The second says the data is public, so no credentials are needed.
Set them before importing rasterio, or GDAL may have already cached its configuration.
3. Read a window, not a file
import pyproj
import rasterio
from rasterio.windows import from_bounds
STUDY = (-2.30, 53.44, -2.18, 53.52)
with rasterio.open(item.assets["red"].href) as src:
print(f"{src.width} x {src.height}, {src.crs}, {src.dtypes[0]}")
to_scene = pyproj.Transformer.from_crs("EPSG:4326", src.crs, always_xy=True)
left, bottom = to_scene.transform(STUDY[0], STUDY[1])
right, top = to_scene.transform(STUDY[2], STUDY[3])
window = from_bounds(left, bottom, right, top, src.transform)
red = src.read(1, window=window)
print(red.shape)10980 x 10980, EPSG:32630, uint16
(898, 787)
The bbox you searched with is WGS84; the scene is UTM. Reprojecting the corners is mandatory — passing degrees to from_bounds on a UTM transform produces a window somewhere near the origin, and rasterio will happily return an empty array for it.
4. Mask with the scene classification band
Sentinel-2 L2A ships an scl asset — a per-pixel classification at 20 m:
CLASSES = {
0: "no data", 1: "saturated", 2: "dark area", 3: "cloud shadow",
4: "vegetation", 5: "not vegetated", 6: "water", 7: "unclassified",
8: "cloud medium", 9: "cloud high", 10: "thin cirrus", 11: "snow",
}
KEEP = [4, 5, 6, 7] # everything that is a real surface observationClasses 3, 8, 9 and 10 are cloud and shadow. Keeping them is what puts a bright white blob in the middle of your NDVI map with a value of about zero.
The SCL band is 20 m and the visible bands are 10 m, so it must be resampled to match — with nearest neighbour, because the values are class codes and averaging two class codes produces a third, meaningless one. This is the same trap as choosing a resampling method.
5. Convert to reflectance, then check it is physically possible
Reflectance is stored as uint16. The conversion is nominally DN / 10000, and since processing baseline 04.00 Sentinel-2 also carries an offset of −1000. But the metadata is not always right about which has already been applied:
print(item.assets["red"].extra_fields.get("raster:bands"))[{'nodata': 0, 'data_type': 'uint16', 'spatial_resolution': 10,
'scale': 0.0001, 'offset': -0.1}]
Apply that offset to this product and NDVI comes out with a median above 1.0 — impossible. Check against physics instead:
veg = scl == 4
print("red DN on vegetation:", np.percentile(red[veg], [5, 50, 95]).round(0))red DN on vegetation: [ 239. 601. 1022.]
A median of 601 is 6% reflectance at DN / 10000 — exactly right for vegetation in the red. Subtract another 1000 and it becomes negative, which no surface does. The offset is already applied in this product, and the declared metadata is stale.
Code examples
Example 1 — a complete windowed NDVI, masked and verified
import os
os.environ["GDAL_DISABLE_READDIR_ON_OPEN"] = "EMPTY_DIR"
os.environ["AWS_NO_SIGN_REQUEST"] = "YES"
import numpy as np
import pyproj
import rasterio
from rasterio.enums import Resampling
from rasterio.windows import from_bounds
from pystac_client import Client
STUDY = (-2.30, 53.44, -2.18, 53.52) # WGS84 degrees
SURFACE_CLASSES = [4, 5, 6, 7] # veg, bare, water, unclassified
def read_window(item, asset, bounds, *, out_shape=None):
"""Read only the pixels covering `bounds` (WGS84) from a remote COG."""
with rasterio.open(item.assets[asset].href) as src:
to_scene = pyproj.Transformer.from_crs("EPSG:4326", src.crs, always_xy=True)
left, bottom = to_scene.transform(bounds[0], bounds[1])
right, top = to_scene.transform(bounds[2], bounds[3])
window = from_bounds(left, bottom, right, top, src.transform)
kwargs = {"window": window}
if out_shape is not None:
kwargs.update(out_shape=out_shape, resampling=Resampling.nearest)
data = src.read(1, **kwargs)
return data, src.crs, src.window_transform(window)
catalog = Client.open("https://earth-search.aws.element84.com/v1")
item = next(catalog.search(
collections=["sentinel-2-l2a"], bbox=list(STUDY),
datetime="2026-06-01/2026-08-01",
query={"eo:cloud_cover": {"lt": 20}},
).items())
red, crs, transform = read_window(item, "red", STUDY)
nir, _, _ = read_window(item, "nir", STUDY)
scl, _, _ = read_window(item, "scl", STUDY, out_shape=red.shape) # 20 m -> 10 m, nearest
usable = np.isin(scl, SURFACE_CLASSES) & (red > 0) & (nir > 0)
r = red.astype("float32") / 10000.0
n = nir.astype("float32") / 10000.0
ndvi = np.where(usable, (n - r) / (n + r), np.nan)
lo, hi = np.nanmin(ndvi), np.nanmax(ndvi)
assert -1.0 <= lo and hi <= 1.0, f"NDVI out of range [{lo:.3f}, {hi:.3f}] — check the scale factor"
print(f"{item.id} {red.shape} {usable.mean():.1%} usable")
print(f"NDVI p5 {np.nanpercentile(ndvi, 5):.3f} median {np.nanmedian(ndvi):.3f} "
f"p95 {np.nanpercentile(ndvi, 95):.3f}")S2B_30UWE_20260712_0_L2A (898, 787) 99.7% usable
NDVI p5 0.058 median 0.322 p95 0.824
The assert is the most valuable line in the block. It costs nothing and it is the only thing standing between a wrong scale factor and a published map.
Example 2 — checking the result against the classification
for code, label in [(4, "vegetation"), (5, "not vegetated"), (6, "water")]:
mask = (scl == code) & usable
if mask.sum():
print(f"SCL {code} {label:15} n={mask.sum():7,} "
f"median NDVI {np.nanmedian(np.where(mask, ndvi, np.nan)):6.3f}")SCL 4 vegetation n=238,122 median NDVI 0.651
SCL 5 not vegetated n=462,968 median NDVI 0.211
SCL 6 water n= 3,468 median NDVI -0.001
This is the check that makes the whole pipeline trustworthy, and it takes three lines. Vegetation should sit well above 0.5, bare and built surfaces around 0.1–0.3, water at or below zero. Water at −0.001 is very nearly a laboratory result.
If water came back at 0.4, the bands are swapped. If vegetation came back at 0.2, the scale factor is wrong. Neither would raise.
Example 3 — saving the window as a georeferenced GeoTIFF
profile = {
"driver": "GTiff",
"height": ndvi.shape[0],
"width": ndvi.shape[1],
"count": 1,
"dtype": "float32",
"crs": crs,
"transform": transform, # the WINDOW transform, not the scene's
"nodata": np.nan,
"compress": "deflate",
"tiled": True,
"blockxsize": 256,
"blockysize": 256,
}
with rasterio.open(f"ndvi_{item.datetime.date()}.tif", "w", **profile) as dst:
dst.write(ndvi.astype("float32"), 1)
dst.update_tags(
source=item.id,
acquired=item.datetime.isoformat(),
cloud_cover=item.properties["eo:cloud_cover"],
scale_note="DN/10000; offset already applied in source COG",
)
with rasterio.open(f"ndvi_{item.datetime.date()}.tif") as check:
print(check.crs, check.res, check.tags()["source"])EPSG:32630 (10.0, 10.0) S2B_30UWE_20260712_0_L2A
src.window_transform(window) is the line people forget. Writing the scene's transform with the window's array puts your 898×787 patch at the scene's origin — 60 km away — with no error at all. See raster and vector do not line up for what that looks like downstream.
The update_tags call embeds provenance in the file itself, so the derived raster still knows which scene and which scale convention produced it.
Explanation
Why the window read is so much faster
The asset is a Cloud-Optimized GeoTIFF: internally tiled, with overviews, header first. GDAL reads the header with one small range request, computes which 1024×1024 tiles intersect your window, and fetches only those.
For a 10,980² scene and an 898×787 window that is under 1% of the file — about six seconds over a normal connection, against several minutes to download the whole band. Nothing is cached to disk unless you ask for it.
The corollary: do not download and then clip. The clip is the read.
Why the SCL band must be resampled with nearest neighbour
SCL values are category codes: 4 is vegetation, 8 is cloud. Resampling from 20 m to 10 m with bilinear interpolation averages neighbouring codes, so a vegetation pixel next to a cloud pixel becomes 6 — water.
Nearest neighbour is the only correct choice for categorical rasters, and it is not the default in every API. This is the same argument as choosing a resampling method for a raster, with sharper consequences because the wrong answer is a plausible-looking class rather than a blurred number.
Why cloud cover in the metadata is not enough
eo:cloud_cover is a single percentage for a 110 km × 110 km granule. Your study area might be 10 km across — 0.8% of it. A scene at 5% cloud can be entirely clear over your city, or entirely obscured, and the number does not distinguish them.
Compute your own, over your own window:
cloudy = np.isin(scl, [3, 8, 9, 10]).mean()
print(f"scene {item.properties['eo:cloud_cover']:.1f}% cloud, study area {cloudy:.1%}")scene 0.1% cloud, study area 0.0%
Use the scene-level figure to shortlist, and the per-window figure to decide.
Why to assert the physical range
Every scale-factor mistake produces arithmetic that runs cleanly and output that is impossible. NDVI is bounded to [−1, 1] by construction whenever both reflectances are non-negative — so a value of 1.04 is proof that at least one input went negative.
That makes the assertion a complete test for a whole class of error: wrong offset, wrong scale, swapped bands, integer overflow. One line, no false positives, and it fails at the point of the mistake rather than three steps downstream in a map nobody can explain.
Edge cases or notes
- Set the GDAL environment variables before importing rasterio. Afterwards, the configuration may already be cached.
- Asset names differ between catalogues. Earth Search uses
red/nir; the Planetary Computer usesB04/B08. Inspectitem.assetsrather than hard-coding. - Planetary Computer assets need signing:
item = planetary_computer.sign(item)before the hrefs resolve. Earth Search does not. - Bands have different native resolutions — 10 m, 20 m and 60 m. Always resample to a common grid explicitly, and choose the method by whether the data is categorical.
nodatais 0 for Sentinel-2 reflectance, and 0 is also a legitimate dark pixel in principle. Maskingred > 0is standard practice and loses almost nothing.- A scene footprint can contain nodata over your area. Orbit-edge granules are partly empty; check the fraction of zeros, not just the footprint.
s2:nodata_pixel_percentageand thes2:*_percentageproperties are in the item metadata and let you screen for this before reading anything.- Reruns are not free. Each windowed read is a fresh set of HTTP requests. Cache the derived arrays, not the scenes.
Internal links
- STAC explained: how satellite imagery catalogues work — the model behind the API
- The raster data model explained — dtype, NoData and transform
- Raster resampling explained — why SCL must use nearest neighbour
- Raster and vector do not line up in Python — what a wrong window transform produces
- How to clip a raster to a polygon in Python — refining the window to a real boundary
- How to merge and mosaic rasters with rasterio — when one scene does not cover the area
- How to calculate zonal statistics in Python — summarising the NDVI you just made
- GIS data sources explained — recording which scene produced which output
FAQ
Do I have to download the whole scene?
No. Assets are Cloud-Optimized GeoTIFFs, so a windowed read fetches only the internal tiles covering your area — typically under 1% of the file.
Why is my NDVI above 1?
A scale-factor error. Reflectance went negative somewhere, usually because an offset was applied that the data had already had applied. Check vegetation red reflectance: it should be a few percent, never negative.
Which bands do I need for NDVI?
Red and near-infrared — red and nir on Earth Search, B04 and B08 elsewhere. Both are 10 m for Sentinel-2, so no resampling is needed between them.
How do I mask clouds?
Read the scl asset and keep classes 4, 5, 6 and 7. Resample it to your working resolution with nearest neighbour, never bilinear.
Why does the search bbox not match the scene CRS?
STAC always uses WGS84 for bbox. The imagery is stored in UTM. Transform the corners before building a read window.
Can I search for something other than optical imagery?
Yes — radar, elevation and land cover collections all live in the same catalogues and search identically. cop-dem-glo-30 is a global DEM on Earth Search.
Should I trust the scale and offset in the metadata?
Verify it. Metadata can be stale, as it is for this product. Check that reflectance values are physically plausible and that any index lands in its defined range.