Download Sentinel-2 L2A imagery from the Copernicus Data Space Ecosystem, mosaic it onto a single grid, and compute spectral indices over it — from a command line or a web interface.
Built for agricultural monitoring at provincial scale: tens of MGRS tiles, 10 m resolution, mosaics that cross UTM zone boundaries.
Python 3.11 or newer.
pip install -e ".[app,dev]"rasterio ships its own GDAL, so no system GDAL install is needed. Verify:
sen2auto configcp .env.example .envThen edit .env:
| Variable | Meaning |
|---|---|
SEN2AUTO_CDSE_USERNAME / SEN2AUTO_CDSE_PASSWORD |
Free account from dataspace.copernicus.eu. Needed to download; searching is public. |
SEN2AUTO_DATA_ROOT |
Where products, mosaics, indices and job state live. |
SEN2AUTO_TARGET_CRS |
Blank to use the AOI's own UTM zone. |
SEN2AUTO_RESOLUTION |
Output pixel size in metres. Default 10. |
SEN2AUTO_BLOCK_SIZE |
Processing block edge in pixels. Default 2048. |
.env is gitignored. Credentials never appear in a job request, a log, or a run
directory — the worker reads them from the environment.
Budget the disk. A Sentinel-2 product is about 1 GB zipped and roughly the same again extracted. A two-band 10 m mosaic of a Pakistani province is around 10 GB; twelve bands is closer to 100 GB.
download keeps two scenes per tile by default, so the compositor has a second
observation to fall through to where the best one is cloudy. That doubles the download.
--per-tile 1 halves it again, at the cost of every cloud becoming a hole rather than a
substitution.
streamlit run app/Home.pyFive pages: Search (pick an area and dates, see what exists), Download, Process (mosaic, then indices — or all three stages as one run, on its Everything at once tab), Analysis (maps, histograms, statistics), and Jobs.
The app never processes anything itself. Each page builds a job request and hands it to the CLI as a detached process, then reads the run directory. Closing the browser does not stop a job; reopening the page picks it back up — the run it shows comes from the run directory, not from the browser session, so a refresh does not lose it.
Settings are read once per server. After editing .env, use Reload configuration
in the sidebar rather than reloading the browser.
# What is available?
sen2auto search --aoi punjab.geojson --start 2026-01-01 --end 2026-01-31 --max-cloud 20
# A bounding box works too: minx,miny,maxx,maxy in WGS84
sen2auto search --aoi 74.0,31.8,74.6,32.4 --start 2026-01-01 --end 2026-01-31
# Will it actually cover the AOI? Ask before downloading a hundred gigabytes.
sen2auto coverage --aoi punjab.geojson --start 2026-01-01 --end 2026-01-31 --max-cloud 20
# Fetch the best scenes per tile
sen2auto download --aoi punjab.geojson --start 2026-01-01 --end 2026-01-31 --max-cloud 20
# Composite onto one grid
sen2auto mosaic --aoi punjab.geojson --name punjab-jan --bands B04,B08
# Compute indices
sen2auto index --name punjab-jan --indices NDVI,NDRE,NDMI
# Or all three at once
sen2auto run --aoi punjab.geojson --name punjab-jan \
--start 2026-01-01 --end 2026-01-31 --bands B04,B08 --indices NDVIJob management:
sen2auto runs # recent runs and their state
sen2auto logs <run-id> # a run's events
sen2auto cancel <run-id> # stop at the next block boundary
sen2auto indices # what indices are available
sen2auto config # resolved settings and free spaceSialkot and Narowal districts, January 2026, NDVI. Four commands, start to finish:
sen2auto search --aoi sialkot_narowal.geojson --start 2026-01-01 --end 2026-01-31 --max-cloud 30
sen2auto download --aoi sialkot_narowal.geojson --start 2026-01-01 --end 2026-01-31 --max-cloud 30
sen2auto mosaic --aoi sialkot_narowal.geojson --name sialkot-jan --bands B04,B08
sen2auto index --name sialkot-jan --indices NDVIsearch prints a table of matching products with tile, date, cloud percentage, orbit
and size, then a count of tiles and distinct dates. index prints mean, standard
deviation, min, max and the valid fraction per index. Actual figures depend entirely
on the scenes available for your area and window.
Outputs land under $SEN2AUTO_DATA_ROOT:
downloads/ raw .zip products
extracted/ unpacked .SAFE directories
mosaics/<name>/ B04.tif, B08.tif, SCL.tif, stack.vrt
indices/<name>/ NDVI.tif
runs/<run-id>/ request.json, status.json, events.jsonl, result.json
Re-running the same --name adds to the set rather than starting over.
| Formula | Bands | Use | |
|---|---|---|---|
| NDVI | (NIR − RED) / (NIR + RED) | B08, B04 | General vegetation vigour |
| NDRE | (NIR − RE1) / (NIR + RE1) | B08, B05 | Late-season canopy; less saturation than NDVI |
| GNDVI | (NIR − GREEN) / (NIR + GREEN) | B08, B03 | Chlorophyll-sensitive |
| SAVI | ((NIR − RED)/(NIR + RED + L)) × (1 + L) | B08, B04 | Sparse cover; L = 0.5 |
| MSAVI2 | self-adjusting soil correction | B08, B04 | Sparse cover without choosing L |
| EVI | 2.5(NIR − RED)/(NIR + 6RED − 7.5BLUE + 1) | B08, B04, B02 | Dense canopy; aerosol-resistant |
| CIre | (NIR / RE1) − 1 | B08, B05 | Chlorophyll content |
| CIgreen | (NIR / GREEN) − 1 | B08, B03 | Chlorophyll content |
| NDWI | (GREEN − NIR) / (GREEN + NIR) | B03, B08 | Open water |
| NDMI | (NIR − SWIR1) / (NIR + SWIR1) | B08, B11 | Canopy water content |
| NDBI | (SWIR1 − NIR) / (SWIR1 + NIR) | B11, B08 | Built-up area |
Adding one is a single entry in sen2auto/indices/definitions.py; masking and nodata
handling are the registry's job, not the index's.
Seven things this pipeline is careful about, because each of them silently corrupts results when it is not — and silence is the common thread. Every one of these failures produces a mosaic that looks like a mosaic.
The BOA offset is applied. Reflectance is (DN + BOA_ADD_OFFSET) / BOA_QUANTIFICATION_VALUE, with the offset read per product from MTD_MSIL2A.xml. It
is −1000 for processing baseline ≥ 04.00 and absent before it. Because the offset is
additive it does not cancel in a normalised difference, so an uncorrected NDVI is
biased everywhere — for a typical crop pixel, by roughly 0.17.
Nodata stays nodata. Index rasters are float32 with NaN nodata. NaN propagates
through arithmetic and cannot collide with a real index value; 0.0 is
indistinguishable from bare soil, and -9999 survives a clip only by accident.
Validity is computed once at read time and carried through every stage rather than
re-derived as array > 0.
Coverage is verified, not assumed. The catalogue's footprint for a partially-imaged
granule is the actual imaged polygon, not the MGRS tile extent, so both "does this set
cover the AOI" and "how much of its tile did this scene observe" are answerable before
anything is downloaded. A granule clipped by the satellite's swath edge is scored down
rather than ranking level with a full one; a tile whose every acquisition is over the
cloud threshold triggers a second search without it, because a 70%-cloud scene
contributes 30% clear ground and a hole contributes none; and whatever still cannot be
covered is reported and written out as coverage_gap.geojson. Counting the tiles you
have says nothing about the area you need.
Cloud is part of the compositing decision, not a later correction. Each scene's SCL
is folded into its validity while compositing, so a cloudy pixel loses to a clear
observation from a lower-ranked scene. Masking only at index time cannot do this: by
then the clear pixel has been discarded, so masking can only turn cloud into a hole.
Validity as DN != 0 is a fill mask, and it reports cloud, cirrus and shadow as good
ground.
Overlaps are resolved by scene quality. Where tiles overlap, each pixel comes from the first valid observation in the best-ranked scene — ranked on cloud cover, distance from the target date, orbit consistency, and how much of its tile the scene actually imaged. Brightest-pixel-wins systematically prefers cloud at every seam.
Scenes are levelled where they meet. One gain and one offset per scene per band are
solved by least squares from the overlaps, measured on clear land only and anchored on
the best-ranked scene. Two correct acquisitions still differ in sun elevation, view
geometry and aerosol load, and first-valid-wins turns that difference into a hard step
exactly on a tile boundary. --rule feathered additionally ramps across the seam
instead of stepping. This is relative normalisation: it makes seams disappear, it does
not make either scene more absolutely correct.
Interior gaps stay gaps. Nothing is filled with a band mean or any other fabricated value. A hole in the data is visible as a hole.
Two further details worth knowing:
- Snow is not masked by default. It is a valid ground observation; masking it
would quietly delete the mountain margin of an AOI rather than removing cloud. Pass
--mask-snowwhen snow is genuinely a nuisance class. - The SCL is mosaicked as a companion raster on the same grid, so tightening cloud
policy costs an index recompute rather than a full re-mosaic. Relaxing it needs a
re-mosaic, since the bands themselves are masked as they are composited — that is
the point, and
--no-mask-cloudsturns it off. - The cloud mask is buffered by 100 m (
--cloud-buffer) before it is applied. SCL under-detects cloud edges and thin cirrus, so an unbuffered mask leaves a bright fringe around every cloud it does find.
sen2auto/
config.py settings from .env / environment
geo.py AOI loading, CRS routing
cdse/ auth, catalogue search, scene selection, coverage, resumable download
safe/ .SAFE parsing, band registry, scene classification
raster/ target grid, warping, GeoTIFF profiles, alignment, VRT stacking
mosaic/ compositing rules, footprint index, seam normalisation, the engine
indices/ index registry and windowed evaluation
runs/ run state, progress events, cancellation (standard library only)
pipeline/ job specs, stages, and the process boundary that runs them
cli.py the interface everything else is a client of
app/ Streamlit pages and components
tests/ 416 tests, including a synthetic .SAFE generator
sen2auto/runs/ deliberately imports nothing outside the standard library. Both the
worker and the UI use it, and the UI has no business loading GDAL to read a progress
file — a broken GDAL should not take down the page that would tell you so.
Mosaicking is destination-driven: the canvas is partitioned into blocks, and for each block the engine asks which scenes overlap it, warps only that block's pixels, composites in memory, and writes the block exactly once. That keeps compression on, makes resume exact at block granularity, and holds peak memory at a few blocks regardless of how large the area is.
Everything warps onto one TargetGrid whose origin is snapped to the resolution.
rasterio.warp.calculate_default_transform is never used: it picks its own
resolution and origin (10.018985 m/px for zone-42 Punjab, in one measured case), which
accumulates about 62 pixels of drift across a province when the result is pasted into
a canvas that assumes exactly 10 m.
Downloads resume with an HTTP Range request rather than starting over. Mosaics
record completed blocks in a checkpoint and skip them on a re-run. Cancelling writes a
sentinel file that the worker checks at each block boundary, so it stops with a
consistent checkpoint rather than mid-write.
Re-running a cancelled job resumes it.
pytest416 tests, no network and no real imagery. Products are synthesised: a real .SAFE
directory tree with namespaced XML and georeferenced rasters, 36×36 px instead of
10980×10980, so the actual parsing and raster code paths run rather than being mocked.
HTTP is mocked with responses.