Fast, Rust-backed best-pixel cloud-free composites straight from STAC/COG sources — give it a bbox + date range, get cloud-free multiband reflectance as numpy arrays (no downloads, no GeoTIFF writes). Supports four STAC sources out of the box:
| Endpoint | Collection(s) | Bands | Quality mask | Output CRS |
|---|---|---|---|---|
pc (default) |
Sentinel-2 L2A | 12 (full S2 set) | SCL | UTM (from AOI) |
earth-search |
Sentinel-2 L2A | 12 (full S2 set) | SCL | UTM (from AOI) |
hls |
HLS L30 + S30 (combined harmonized pool) | 7 common bands | Fmask | UTM (from AOI) |
mcd43a4 |
MCD43A4 v6.1 (MODIS NBAR, daily) | 6 (no rededge) | Mandatory_Quality_Band1 | MODIS Sinusoidal (native) |
- Asynchronously fetches scenes from Planetary Computer (default), the Element84 earth-search S2 L2A mirror, or PC's Harmonized Landsat- Sentinel-2 (HLS v2.0) via STAC search + range-request COG reads.
- Selects the top-k clearest observations per MGRS-tile-aware chunk using a coarse quality scout pass (SCL classes for S2 L2A; Fmask bit-flags for HLS).
- Decodes DEFLATE-compressed COG tiles via libdeflate, undoes TIFF predictor=2, and reprojects to a user-specified UTM grid with an AVX2/FMA bilinear resampler.
- Composes a best-pixel monthly composite ranked by per-pixel quality.
For a 100 × 100 km AOI at 60 m resolution, top-k=3 from Planetary Computer, on a single 16-core Zen 4 node:
| Workload | Time |
|---|---|
| 1 year × 1 month | ~2 s |
| 5 years × 1 month, sequential | ~11 s |
| 5 years × 1 month, 5-thread parallel | ~6 s |
The 5-way parallel floor is set by network throughput to PC, not CPU.
pip install bestpixelA single abi3 wheel runs on CPython 3.9–3.14. Or build locally with maturin:
pip install maturin
cd surface_priors_rs # the crate directory
maturin develop --release --features pythonimport bestpixel as bp
out = bp.build_composite(
bbox=(30.5, 30.5, 31.6, 31.5), # west, south, east, north (WGS84)
datetime="2024-07-01/2024-07-31", # STAC datetime range
resolution=60.0, # metres
endpoint="pc", # "pc" | "earth-search" | "hls" | "mcd43a4" | "auto"
bands=["coastal", "blue", "green", "red", "nir", "swir16", "swir22"],
)
red = out["bands"]["red"] # uint16 ndarray, (H, W)
quality = out["quality"] # uint16, 0=clear, 1=marginal, 2=dark, 65535=nodata
print(out["grid"]) # bounds, epsg, transform — for georeferencingValues are int16 DN; reflectance = DN / 10000 (the Sentinel-2 N0400
+1000 BOA offset is harmonized internally so all years are comparable).
By default each chunk takes a fixed top_k clearest observations. For
adaptive depth, set coverage_target in (0, 1): the scout's coarse
per-chunk cloud masks drive how many scenes each chunk stacks — cloudy
chunks pull more, clear chunks stop early — with a min_k redundancy
floor and max_k cap. Add windowed_fetch=True to read each scene only
over the chunks that need it (far fewer bytes when depth is concentrated,
e.g. an under-observed swath-edge corner):
out = bp.build_composite(
bbox=(30.5, 30.5, 31.6, 31.5),
datetime="2024-07-01/2024-07-31",
coverage_target=0.95, min_k=2, max_k=6, # adaptive per-chunk depth
windowed_fetch=True, # read only the windows that need depth
)Available band names (stable across endpoints):
coastal, blue, green, red, rededge1, rededge2, rededge3, nir, nir08, nir09, swir16, swir22. SCL / Fmask is consumed internally to derive
the quality raster — kept as a discrete class label (nearest
resampling) all the way through, so quality buckets stay categorical.
Pass bands=None (or omit) to fetch all bands the endpoint supports
(12 for S2 L2A, 7 for HLS).
build_monthly_composites(bbox, years=[...], months=[...], ...)
produces one composite per (year, month) combination from a single
STAC search + scout pass. Cleaner than looping build_composite in
Python, and faster because the search + scout work is amortised.
out = bp.build_monthly_composites(
bbox=(30.5, 30.5, 31.6, 31.5),
years=[2018, 2019, 2020],
months=[6, 7, 8], # June / July / August
resolution=60.0,
top_k=3,
endpoint="pc",
)
# 9 composites returned as a list, each carrying its own year/month:
for r in out:
print(r["year"], r["month"], r["bands"]["red"].shape)Per-period tasks are spawned concurrently on the shared runtime, so
the 9-period example above runs in ~22 s wall — about the same as
sequential build_composite calls but with one STAC search instead
of nine.
endpoint="hls" pulls from PC's hls2-l30 + hls2-s30 collections
in a single combined pool and composites them together. HLS already
applies the Roy et al. c-factor NBAR-style normalisation, so Landsat-8/9
OLI and Sentinel-2 MSI observations are bit-comparable.
Only the 7 harmonized common bands are exposed: coastal, blue, green, red, nir, swir16, swir22. The "nir" band uses Landsat's B05 / Sentinel-2's
B8A (narrow NIR — the harmonized choice from Roy 2021), not S2's B08
broad NIR.
out = bp.build_composite(
bbox=(30.5, 30.5, 31.6, 31.5),
datetime="2024-07-01/2024-07-31",
resolution=60.0,
endpoint="hls",
)
# 5-year × 1-month over a 100 km AOI, parallel-5: ~7 sInternally HLS scenes resolve their per-collection asset names
on-the-fly: for the same band slot ("red"), L30 reads B04 and S30
reads B04 (they happen to align), while "nir" reads B05 on L30 and
B8A on S30. Quality scoring uses the bit-packed Fmask: cloud,
cirrus, cloud shadow, snow, and high aerosol bits each weight into a
lower-is-better score that drives best-pixel selection.
endpoint="mcd43a4" pulls the daily 500 m MODIS Nadir BRDF-Adjusted
Reflectance product from PC's modis-43A4-061 collection.
Output stays in MODIS Sinusoidal end-to-end — no on-the-fly
reprojection to UTM. The returned grid dict reports epsg=0 and
carries the proj4 string instead:
out = bp.build_composite(
bbox=(30.5, 30.5, 31.6, 31.5),
datetime="2024-07-01/2024-07-31",
resolution=500.0, # MODIS native
endpoint="mcd43a4",
)
print(out["grid"]["proj4"])
# '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 ...'If you need MCD43A4 stacked against an S2 or HLS UTM grid, gdalwarp
/ rio warp the returned arrays as a post-processing step.
Six harmonized bands are exposed: blue, green, red, nir, swir16, swir22 — mapped to MODIS bands 3/4/1/2/6/7 respectively. MODIS Band 5
(1240 nm) has no S2 / HLS counterpart and is omitted. Quality scoring
uses BRDF_Albedo_Band_Mandatory_Quality_Band1: 0 = full BRDF inversion (best), 1 = magnitude inversion (acceptable), 255 = fill.
To stack MCD43A4 against an S2 or HLS UTM grid, pass
output_crs="utm". The pipeline keeps the source COG reads in
sinusoidal coordinates but bilinearly reprojects to UTM during the
resample step (per-pixel PROJ transform; no extra disk passes).
out = bp.build_composite(
bbox=(30.5, 30.5, 31.6, 31.5),
datetime="2024-07-01/2024-07-31",
resolution=500.0,
endpoint="mcd43a4",
output_crs="utm",
)
print(out["grid"]["epsg"]) # 32636 (UTM 36N)Same flag works on build_monthly_composites. SCL / Fmask / mandatory-
quality stay nearest during reprojection, so categorical labels survive
the CRS change.
build_l1c_composite(...) builds a surface-reflectance composite from
raw L1C top-of-atmosphere scenes that you correct yourself with 6S,
instead of relying on the L2A product's Sen2Cor correction. The recipe
(validated against AERONET): use MAIAC AOD to pick the low-aerosol days,
6S-correct those TOA scenes (MAIAC AOD + per-scene water vapour + geometry),
and best-pixel composite — preferring the lowest-AOD clear pixel.
L1C is JP2-on-requester-pays-s3 (not an HTTP COG), so this path fetches via
Google Earth Engine — pip install 'bestpixel[gee]' and authenticate
earthengine-api (service account via GEE_SERVICE_ACCOUNT /
GEE_SERVICE_ACCOUNT_KEY, or earthengine authenticate). The fetch is
Python (ee.data.getPixels, scout-first windowed); the per-pixel 6S
correction runs in Rust (bestpixel.correct_toa).
It needs a per-scene atmosphere sidecar (MAIAC AOD + water vapour +
geometry → 6S coefficients over a TCWV LUT), built by the scripts under
scripts/ (GEE + SIAC native 6S); see scripts/run_l1c_composite.sh.
import bestpixel as bp
out = bp.build_l1c_composite(
bbox=(31.0, 29.9, 31.5, 30.3), # WGS84 W,S,E,N
datetime=("2022-07-01", "2022-08-01"),
sidecar="cairo_sidecar.json", # 6S coefficients per scene
resolution=60,
low_aod_frac=0.6, # keep the cleanest 60% of days
rank="aod", # lowest-AOD clear pixel wins (cs tie-break)
out="cairo_l1c.tif", # optional scaled int16 GeoTIFF
)
blue = out["bands"]["blue"] # (H, W) surface reflectance
print(out["scenes"], out["grid"]["epsg"])Returns {"bands": {name: (H,W) float32}, "grid": {...}, "count": (H,W), "scenes": [...]}. Bands are the seven the sidecar carries coefficients for
(coastal, blue, green, red, nir08, swir16, swir22). Ranking by AOD (rather
than Cloud Score+, which saturates near 1.0 on clear pixels) keeps the result
deterministic and picks the most reliably-corrected day.
The crate also produces a native binary that writes tiled DEFLATE-compressed GeoTIFFs:
cargo build --release
./target/release/spx-build \
--bbox 30.5 30.5 31.6 31.5 \
--datetime 2024-07-01/2024-07-31 \
--resolution 60 \
--top-k 3 \
--endpoint pc \
--out-dir /tmp/spx-outThe CLI supports a persistent disk cache (--disk-cache <dir>) so
repeated runs over the same AOI skip STAC + scout work.
src/
cog.rs HTTP COG reader (TIFF/IFD, range tiles, DEFLATE+predictor=2)
stac.rs Async STAC search client
endpoint.rs PC / Element84 endpoint config + SAS-token signing
grid.rs UTM grid math, bilinear AVX2 resampler, nearest u8 resampler (SCL)
projx.rs PROJ-backed coordinate transforms
tile_classification.rs Geometry-based exclusive-coverage MGRS chunk classifier
pipeline.rs scout / select / fetch / compose
writer.rs Tiled DEFLATE GeoTIFF output
disk_cache.rs JSON cache for search / scout / partition
py.rs PyO3 module exposing build_composite()
bin/spx_build.rs CLI entrypoint