Draw an area of interest on a map, pick a date range, and get a cloud-masked Sentinel-2 spectral index time series, an animated GIF and a downloadable CSV and GeoTIFF — from open data, with no account and no local imagery.
Answering "how did this field green up last summer?" is a half-day job with the usual tooling.
You search a STAC catalog, work out whether the collection calls the red band red or B04,
discover that scenes after January 2022 carry a -1000 radiometric offset, write the SCL masking
by hand, remember that your AOI straddles two Sentinel-2 tiles so half your dates appear twice
with partial coverage, and then find that the mean you plotted for 14 June came from 3% of the
AOI because the rest was cloud.
None of that is interesting work, and all of it is easy to get subtly wrong in a way that still produces a plausible-looking chart. This packages the whole path — search, stack, mask, reduce, render — behind a map and a date picker, and it shows the valid-pixel fraction next to every point so you can tell a real observation from a number computed off a handful of pixels.
Python 3.12 or newer. Clone the repository and install the dependencies; there is no packaging
step and nothing to pip install from an index.
git clone https://github.com/python-remote-sensing/ndvi-timeseries-web.git
cd ndvi-timeseries-web
python3.12 -m venv .venv
source .venv/bin/activate
pip install -r requirements.txt./run.sh # http://localhost:8501
./run.sh --server.port 8600 # extra arguments go to `streamlit run`Equivalently, streamlit run ndvi_app/app.py from the repository root.
The layout:
┌─ Sentinel-2 index time series ────────────────────────────────────────────────┐
│ │
│ ┌─ Settings ────────┐ ┌─ 1. Area of interest ──────┐ ┌─ 2. Confirm and run ┐│
│ │ Index │ │ ▛▚ ┌───────────────────┐ │ │ ││
│ │ [NDVI ▾] │ │ ▙▟ │ ▄▄▄▄▄▄▄▄ │ │ │ ...or upload a ││
│ │ │ │ ⬛ │ █ drawn █ │ │ │ GeoJSON polygon ││
│ │ Date range │ │ ✎ │ █ AOI █ │ │ │ [ Browse files ] ││
│ │ [2024-04-01] │ │ 🗑 │ ▀▀▀▀▀▀▀▀ │ │ │ ││
│ │ [2024-09-30] │ │ │ OpenStreetMap│ │ │ AOI area ││
│ │ │ │ └───────────────────┘ │ │ 3.5 km² ││
│ │ Max cloud cover │ │ polygon / rectangle tools │ │ Auto resolution: 10 m││
│ │ 0 ──●───── 100 │ │ on the left, layer switch │ │ ││
│ │ 40 % │ │ for satellite basemap │ │ [ Run ] ││
│ │ │ └────────────────────────────┘ └─────────────────────┘│
│ │ ▸ Advanced │ │
│ └───────────────────┘ Scenes 11 │ Dates with data 6/6 │ AOI 3.5 km² │ 10 m │
│ │
│ ┌ Time series ┬ Animation ┬ Scene table ┐ │
│ │ 0.56 ┤ ● ● │
│ │ 0.48 ┤ ● │
│ │ 0.40 ┤ ● ● │
│ │ 0.34 ┤ ● ← point colour = valid │
│ │ └────┬─────┬─────┬─────┬─────┬─────┬── pixel fraction (0–1) │
│ │ 06-01 06-08 06-15 06-22 07-01 07-08 │
│ │ valid ▉▉▉▉ ▃▃▃▃ ▉▉▉▉ ▉▉▉▉ ▉▉▉▉ ▃▃▃▃ ← short bars are partial coverage │
│ └──────────────────────────────────────────────────────────────────────────┘ │
│ │
│ [ Download CSV ] [ Download GIF ] [ Download GeoTIFF ] │
└───────────────────────────────────────────────────────────────────────────────┘
The Animation tab holds the GIF: one frame per scene, a fixed colour ramp so brightness changes mean real index changes, a date label and a colorbar. The Scene table tab lists every STAC item that was used with its id, acquisition time and scene cloud cover.
The same pipeline runs headless, for cron jobs and Prefect flows:
python -m ndvi_app \
--aoi field.geojson \
--start 2024-06-01 --end 2024-07-15 \
--index ndvi --max-cloud 25 \
--out-dir results[ 2%] Validating the area of interest
[ 8%] Searching Earth Search (Element 84 / AWS) for Sentinel-2 L2A scenes
[ 20%] Found 11 scene(s); building the datacube
[ 35%] Applying SCL cloud masking and computing the index
[ 50%] Reading pixels from the cloud-optimized GeoTIFFs
[ 85%] Reducing to an AOI mean per date
[ 95%] Building the composite
[100%] Done
NDVI (vegetation) over 3.5 km2 at 10 m
11 scene(s) from Earth Search (Element 84 / AWS), 6/6 dates with valid pixels
date mean std valid%
-----------------------------------------
2024-06-01 0.4175 0.4179 100.0%
2024-06-06 0.5619 0.3214 34.5%
2024-06-19 0.3866 0.2220 100.0%
2024-07-04 0.3959 0.3201 100.0%
2024-07-11 0.3393 0.3952 100.0%
2024-07-14 0.5179 0.2992 34.5%
wrote csv: results/ndvi.csv
wrote gif: results/ndvi.gif
wrote tif: results/ndvi_composite.tif
wrote chart: results/ndvi_chart.png
The two 34.5% dates are real: that AOI sits on a Sentinel-2 tile boundary and on those passes only one tile covered it. The column is there precisely so you can see that rather than wonder why the curve jumps.
--json gives a machine-readable summary instead:
python -m ndvi_app --aoi field.geojson --start 2024-06-01 --end 2024-07-15 --json --quiet{
"index": "ndvi",
"catalog": "earth-search",
"aoi_area_km2": 3.489,
"resolution_m": 10.0,
"grid": {"height": 224, "width": 160},
"scenes_found": 11,
"dates": 6,
"dates_with_valid_pixels": 6,
"mean_range": [0.339251, 0.561922],
"warnings": [],
"outputs": {},
"series": [
{"date": "2024-06-01", "mean": 0.417468, "valid_fraction": 1.0},
{"date": "2024-06-06", "mean": 0.561922, "valid_fraction": 0.3454}
]
}The CSV carries the full per-date statistics:
date,ndvi_mean,ndvi_median,ndvi_std,ndvi_min,ndvi_max,valid_pixels,aoi_pixels,valid_fraction
2024-06-01,0.417468,0.479959,0.417861,-0.675676,0.999544,34982,34982,1.000000
2024-06-06,0.561922,0.617385,0.321365,-0.973799,0.999565,12083,34982,0.345406
2024-06-19,0.386620,0.381642,0.222035,-0.091300,0.826546,34982,34982,1.000000The AOI bounding box, date range and cloud ceiling become a single STAC query against
Earth Search or the Planetary Computer. The cloud
filter is a query term on eo:cloud_cover, so cloudy scenes are excluded by the API and are
never fetched at all.
The matching items go to stackstac.stack, which returns a lazy dask-backed
(time, band, y, x) array. No pixel is read until the reduction runs, and only the AOI window
of each Cloud-Optimized GeoTIFF is ever requested — HTTP range reads, not whole scenes.
Earth Search publishes assets as red, nir, swir16. The Planetary Computer publishes the
same bands as B04, B08, B11. Index formulas here are written against canonical names, and
stac_source.CatalogSpec.asset_map does the translation, so adding a catalog is a dict entry
rather than a rewrite of the band math.
ESA processing baseline 04.00 (operational from 25 January 2022) shifted Sentinel-2 L2A digital numbers by -1000 so that negative retrievals survive unsigned encoding. Ignore it and every reflectance is 0.1 too high.
For NDVI that barely matters — a normalized difference is close to scale-invariant. For EVI it
matters a lot, because the +1 soil-adjustment term in the denominator is not. So the per-scene
offset is read from s2:processing_baseline (or from earthsearch:boa_offset_applied when the
catalog says it already handled it) and applied before the 1e-4 reflectance scaling, per
scene, since a multi-year request spans the baseline change.
The Scene Classification Layer is stacked on the same grid, at the same resolution, with the same resampling as the reflectance bands — masking against a separately-loaded SCL is a classic source of one-pixel misregistration. Discarded by default:
| SCL | Class | Masked by default |
|---|---|---|
| 0 | no_data | yes |
| 1 | saturated_or_defective | yes |
| 2 | dark_area_pixels | no |
| 3 | cloud_shadow | yes |
| 4 | vegetation | no |
| 5 | not_vegetated | no |
| 6 | water | no |
| 7 | unclassified | no |
| 8 | cloud_medium_probability | yes |
| 9 | cloud_high_probability | yes |
| 10 | thin_cirrus | yes |
| 11 | snow | yes |
Class 7 (unclassified) is kept on purpose. Sen2Cor assigns it freely over bright bare soil,
and dropping it biases arid-region time series towards the vegetated pixels that happen to
classify cleanly. stac_source.STRICT_MASK_CLASSES is available if you want it gone.
Masked pixels become NaN, and the AOI polygon is rasterized onto the output grid so that only pixels genuinely inside it count. That matters for the denominator: with a bounding-box denominator a diagonal AOI would never report more than about 50% valid even on a perfectly clear day.
Each date then yields mean, median, std, min, max, the valid pixel count and the valid fraction. A fully-masked date is kept, with NaN statistics and a valid fraction of 0 — dropping it would hide the fact that the satellite did observe that day, which is exactly what someone checking for a data gap needs to see.
Scenes acquired on the same calendar day (adjacent orbits over one AOI) are merged into a single row, weighted by each scene's valid pixel count, so a tile-boundary AOI produces one honest observation per pass rather than a sawtooth.
Two limits are enforced before any network call:
- Area. Above 2500 km² the request is refused with the actual area and the limit in the message. At 10 m that is already 25 million pixels per band per date.
- Pixel budget. Below the area limit but still large, the output grid is coarsened to a multiple of 10 m so one date stays under about 4 million pixels, and you are told it happened.
The output CRS is pinned to the AOI's UTM zone rather than inferred from item metadata. Earth
Search v1 publishes proj:code where older stackstac releases expect proj:epsg, and letting
the stack guess fails outright with "cannot pick a common CRS".
All data logic lives in pure, importable modules; ndvi_app/app.py only reads widgets and draws
results. That is what makes the suite hermetic — the tests exercise the real code paths without
a Streamlit runtime and without the network.
| Module | Responsibility |
|---|---|
ndvi_app/aoi.py |
GeoJSON parsing, geometry validation, geodesic area, guards, UTM zone |
ndvi_app/indices.py |
Index registry, band math, reflectance scaling |
ndvi_app/stac_source.py |
Catalog registry, query construction, cube loading, SCL masking |
ndvi_app/timeseries.py |
AOI rasterization, per-date reduction, compositing, CSV |
ndvi_app/render.py |
GIF frames, chart PNG, in-memory GeoTIFF |
ndvi_app/pipeline.py |
Orchestration shared by the app and the CLI |
ndvi_app/app.py |
Streamlit view layer only |
ndvi_app/cli.py |
argparse front end |
| Flag | Default | Meaning |
|---|---|---|
--aoi PATH |
required | GeoJSON file with one polygon in EPSG:4326 |
--start YYYY-MM-DD |
required | First acquisition date, inclusive |
--end YYYY-MM-DD |
required | Last acquisition date, inclusive |
--index {ndvi,ndwi,ndmi,evi} |
ndvi |
Spectral index to compute |
--catalog {earth-search,planetary-computer} |
earth-search |
STAC endpoint |
--max-cloud FLOAT |
40.0 |
Scene cloud cover ceiling, percent |
--mask-classes LIST |
0,1,3,8,9,10,11 |
SCL classes to discard, numbers or names |
--resolution FLOAT |
auto | Output pixel size in metres |
--max-scenes INT |
60 |
Ceiling on scenes stacked |
--min-valid-fraction FLOAT |
0.0 |
Drop dates below this valid fraction |
--composite {median,mean,max,min} |
median |
Temporal reducer for the GeoTIFF |
--no-merge-same-day |
off | One row per scene instead of per calendar day |
--out-dir DIR |
— | Write all four products into a directory |
--csv/--gif/--tif/--chart PATH |
— | Explicit path for one product |
--frame-ms INT |
700 |
Milliseconds per GIF frame |
--json |
off | Machine-readable summary on stdout |
--quiet |
off | Suppress progress on stderr |
--list-indices |
— | List the available indices and exit |
Exit codes: 0 success, 1 no usable data or a runtime failure, 2 bad arguments.
| Name | Formula | Bands |
|---|---|---|
ndvi |
(NIR − Red) / (NIR + Red) |
nir, red |
ndwi |
(Green − NIR) / (Green + NIR) |
green, nir |
ndmi |
(NIR − SWIR₁.₆) / (NIR + SWIR₁.₆) |
nir, swir16 |
evi |
2.5 (NIR − Red) / (NIR + 6 Red − 7.5 Blue + 1) |
nir, red, blue |
from datetime import date
from ndvi_app.pipeline import run_timeseries
result = run_timeseries(
{"type": "Polygon", "coordinates": [[[9.70, 45.05], [9.72, 45.05],
[9.72, 45.07], [9.70, 45.07],
[9.70, 45.05]]]},
date(2024, 6, 1), date(2024, 7, 15),
index="ndvi", max_cloud=25,
progress=lambda fraction, message: print(f"{fraction:.0%} {message}"),
)
print(result.frame[["date", "mean", "valid_fraction"]])
print(result.usable_dates, "dates with data")result also carries cube (the masked (time, y, x) index array), composite, transform,
crs, aoi_mask, scenes and warnings. The pure helpers are useful on their own:
from ndvi_app.indices import compute_index
from ndvi_app.stac_source import scl_valid_mask
from ndvi_app.timeseries import reduce_timeseries, to_csv_bytes
from ndvi_app.render import render_gif, to_geotiff_bytes- Sentinel-2 L2A only. No Landsat, no Sentinel-1, no harmonized cross-sensor series.
- Optical, so cloud-limited. In a wet season you may get two usable dates a month. The tool reports that honestly rather than interpolating over it.
- Sen2Cor's SCL is imperfect. Thin cirrus and cloud edges leak through, and the classifier confuses bright soil, buildings and snow. Treat a single anomalous date with suspicion; that is what the min/max/std columns are for.
- No gap filling, no smoothing, no phenology fitting. The output is observations, not a fitted curve. Savitzky-Golay or a harmonic fit belongs in your analysis, not in the fetch.
- The area guard is deliberately conservative at 2500 km². Regional or national work wants a Dask cluster and a batch pipeline, not a Streamlit session.
- AOIs crossing the antimeridian or a UTM zone edge are handled by pinning to the centroid's zone; a very wide AOI will be reprojected some distance from its natural zone.
- The GIF has one frame per scene, while the chart has one row per calendar day. On a tile-boundary AOI the animation therefore shows more frames than the chart shows points.
- In-memory reduction. The cube is computed into RAM before reduction, which is why the pixel budget exists. Streaming reductions would lift that limit at the cost of complexity.
Background on the mechanisms this tool automates, from the site it is maintained alongside:
- Masking clouds with the Sentinel-2 SCL band — the class-by-class reasoning behind the default mask set used here.
- Using pystac-client to filter Sentinel-2 imagery by date — how the datetime and
queryterms in the STAC search are built. - stackstac vs odc-stac for STAC to array — the trade-offs behind stacking items into a lazy datacube.
- Calculating NDVI directly from xarray DataArrays — band math that keeps coordinates and dask chunking intact.
- Creating monthly NDVI composites with xarray resample — the next step once you want a regular series rather than irregular acquisitions.
pip install -r requirements-dev.txt
python -m pytest # 216 tests, no network
ruff check .
ruff format --check .The suite is hermetic by construction: STAC responses come from fake clients, datacubes are synthetic xarray arrays on a real UTM grid, and GIF and GeoTIFF outputs are read back with Pillow and rasterio to assert on frame counts, dimensions, CRS, transform and pixel values. CI runs it on Python 3.12, 3.13 and 3.14.
Issues and pull requests are welcome. Please keep data logic out of app.py, add tests that
assert on values rather than on "it ran", and run ruff format . before committing.
MIT — see LICENSE.
Built and maintained alongside Python Remote Sensing & Raster Processing Pipelines.