--- title: "Monitoring reservoir water quality from space with Sentinel-2 and Python (Peak District case study)" description: "Monitor reservoir chlorophyll and turbidity with Sentinel-2 and Python: NDCI, NDTI, cloud and water masks, and a Peak District case study, 2019 to 2026." url: https://www.technicalgis.com/tutorials/reservoir-water-quality-sentinel-2/ author: TechnicalGIS Team publisher: TechnicalGIS (https://www.technicalgis.com/about/) published: 2026-10-15 updated: 2026-10-15 language: en-GB --- # Monitoring reservoir water quality from space with Sentinel-2 and Python (Peak District case study) > Monitor reservoir chlorophyll and turbidity with Sentinel-2 and Python: NDCI, NDTI, cloud and water masks, and a Peak District case study, 2019 to 2026. - Level: Intermediate - Tested with: Python 3.12 - Source code: [Download scripts](https://www.technicalgis.com/code/reservoir-water-quality-sentinel-2/) ## Key takeaways 1. Sentinel-2's red-edge band B5 makes a chlorophyll proxy possible (NDCI), but the result is a relative index, not a concentration in µg/L. 2. Peaty Peak District water is so dark that MNDWI is negative over most of it. A water mask built from the SCL cloud mask, low near-infrared reflectance and a 40 m inward buffer works better. 3. Of 821 scenes from 2019 to 2026, cloud left 122 usable dates. Ladybower had the highest median NDCI (0.049) and Torside the lowest (0.029). In the Upper Derwent reservoirs NDCI rises in April and again in August and September. 4. Against 18 Environment Agency chlorophyll samples taken within three days of a clear image, NDCI showed no relationship (Spearman ρ = −0.03). At 0.5–9 µg/L, the anomalies flag dates to check; they do not replace sampling. ## Before you start - Comfortable with NumPy, pandas and rasterio (the Sentinel-2 Lake Mead tutorial introduces STAC searches and Cloud Optimised GeoTIFFs) - Python 3.10+ with pystac-client, planetary-computer, rasterio, numpy, scipy, pandas, geopandas, osmnx, requests and matplotlib - An internet connection: imagery is streamed from the Microsoft Planetary Computer and needs no account or API key **Short answer:** to monitor reservoir water quality with Sentinel-2 in Python, search the free Sentinel-2 L2A archive with STAC, mask cloud and shadow with the SCL band, keep only clear water pixels well inside the shoreline, and compute indices from the median water spectrum on each date. Use **NDCI** = (B5 − B4) / (B5 + B4) as a chlorophyll-a proxy, **NDTI** or red reflectance as a turbidity proxy, and a band ratio such as B3 / B2 for colour. Then compare each date with the reservoir’s normal for that month. These are **relative** indices: without local calibration and a water-specific atmospheric correction, they flag _when_ and _where_ to look, not how many µg/L of chlorophyll are present. The Peak District reservoirs supply drinking water to Sheffield, Derby, Nottingham, Leicester and Manchester. Like most upland reservoirs in Britain, they are small, narrow, steep-sided and stained brown by peat. That makes them a hard, realistic test for satellite water-quality monitoring. In this tutorial we follow four of them from January 2019 to October 2026: **Howden, Derwent and Ladybower** in the Upper Derwent Valley (Severn Trent), and **Torside** in Longdendale (United Utilities), using every Sentinel-2 image of the area. You will learn how to: 1. get reservoir outlines from **OpenStreetMap** with `osmnx`, and shrink them to avoid mixed shoreline pixels; 2. search and stream **Sentinel-2 L2A** from the Microsoft Planetary Computer, handling the processing-baseline offset; 3. build a cloud, shadow and **water mask** that works for dark, peaty water, and see why MNDWI does not; 4. compute **NDCI, MCI, NDTI and a colour ratio**, with quality control; 5. build a **multi-year time series**, a **monthly climatology** and a ranked list of **anomalous dates**; 6. map NDCI and turbidity for one date and export a **tidy CSV** per reservoir; 7. check the results against public **Environment Agency** samples. All outputs below are real results from [`reservoir_water_quality.py`](https://www.technicalgis.com/code/reservoir-water-quality-sentinel-2/), run on 15 October 2026 with Python 3.12, pystac-client 0.9.0, rasterio 1.5.2 and osmnx 2.1.1. ## What can a satellite see in a reservoir? A satellite sees only the light that leaves the water surface. Several things in the water change its colour: - **Chlorophyll-a** in algae absorbs red light near 665 nm (Sentinel-2 band B4) and, as algae build up, causes a reflectance peak near 705 nm (B5). NDCI and MCI measure that peak. - **Suspended sediment** scatters light and raises reflectance, most visibly in the red (B4). NDTI compares red with green. - **Coloured dissolved organic matter (CDOM)**, the brown stain from peat, absorbs blue light strongly and makes the water very dark. It is a real water-treatment issue in the Peak District, because colour must be removed before the water can be supplied. Sentinel-2 has 10–20 m pixels, a red-edge band at 705 nm and a revisit of 2–5 days, which makes it the best free sensor for small lakes. Sentinel-3 OLCI has better water bands but 300 m pixels, wider than these reservoirs. Landsat 8 and 9 have no red-edge band, so they cannot compute NDCI. ## Step 1: get the reservoir outlines from OpenStreetMap OpenStreetMap maps reservoirs with the tag `water=reservoir`. `osmnx` downloads every matching polygon in a bounding box from the Overpass API. Public Overpass servers are often busy, so the script retries with backoff and a mirror, and caches the result in a GeoPackage: `reservoir_water_quality.py` ```python RESERVOIRS = ["Howden", "Derwent", "Ladybower", "Torside"] # north to south, then Longdendale CRS = "EPSG:32630" # UTM zone 30N, the native grid of Sentinel-2 tile 30UWE INSET_M = 40 # shrink each outline by 40 m (two 20 m pixels) to avoid mixed shoreline pixels osm = ox.features_from_bbox((-1.95, 53.33, -1.62, 53.55), {"water": "reservoir"}) osm = osm[osm.geometry.type.isin(["Polygon", "MultiPolygon"]) & osm["name"].notna()] osm = osm[osm["name"].isin([f"{r} Reservoir" for r in RESERVOIRS])].reset_index() res = gpd.read_file(outlines).to_crs(CRS) res["reservoir"] = res["name"].str.replace(" Reservoir", "") res = res.dissolve("reservoir").loc[RESERVOIRS].reset_index() res["core"] = res.buffer(-INSET_M) ``` The negative buffer is important. A 20 m pixel on the shoreline mixes water with grass, rock or trees, which are hundreds of times brighter in the near-infrared and red edge. Even one such pixel can dominate a water index. Shrinking each outline by 40 m keeps only pixels that are entirely water. `Output — reservoir outlines (OpenStreetMap)` ```text reservoir id area_ha core_ha Howden 955357 62.4 33.6 Derwent 955161 69.5 42.9 Ladybower 1959653 198.1 124.8 Torside 8598540 62.6 42.0 ``` The inset removes 33–46% of each reservoir’s area. That is the price of working with long, narrow lakes. We then rasterise the cores onto one 20 m grid that covers all four reservoirs: `Output — analysis grid` ```text 817 x 738 pixels of 20 m (16.3 x 14.8 km), EPSG:32630 Howden 838 core water pixels Derwent 1073 core water pixels Ladybower 3129 core water pixels Torside 1046 core water pixels ``` ![Sentinel-2 true-colour image of the northern Peak District with white outlines and coloured labels for four reservoirs. Torside, 63 ha, is in the top-left corner in Longdendale. Howden (62 ha), Derwent (70 ha) and Ladybower (198 ha) form a chain down the Upper Derwent Valley on the right. An inset map of Great Britain marks the location in northern England.](https://www.technicalgis.com/_astro/01-locator-map.1c9I_RDw_yglzR.webp) _The four study reservoirs on 12 July 2026, with areas from the OpenStreetMap outlines. The unlabelled dark lakes are other reservoirs that we leave out._ ## Step 2: search the Sentinel-2 archive with STAC The Microsoft Planetary Computer holds the whole Sentinel-2 L2A archive as Cloud Optimised GeoTIFFs. All four reservoirs fall inside one tile, **30UWE**, so we search for that tile only. A few dates appear twice, once from the original processing and once from ESA’s reprocessing, so we keep one scene per date and orbit: `reservoir_water_quality.py` ```python STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1" PERIOD = "2019-01-01/2026-10-14" catalog = pystac_client.Client.open(STAC_URL) search = catalog.search(collections=["sentinel-2-l2a"], bbox=bbox_ll, datetime=PERIOD, query={"eo:cloud_cover": {"lt": 95}, "s2:mgrs_tile": {"eq": "30UWE"}}) items = list(search.items()) # Some dates exist twice (an original and a reprocessed baseline): keep the newest processing. best = {} for it in items: key = (it.datetime.date(), it.properties["sat:relative_orbit"]) if key not in best or it.properties["s2:processing_baseline"] > best[key].properties["s2:processing_baseline"]: best[key] = it items = sorted(best.values(), key=lambda i: i.datetime) ``` `Output — STAC search` ```text 821 Sentinel-2 L2A scenes of tile 30UWE with < 95% cloud, 02 Jan 2019 to 14 Oct 2026 {'02.12': 257, '03.00': 75, '04.00': 51, '05.09': 16, '05.10': 149, '05.11': 171, '05.12': 93, '05.13': 9} ``` The scene-level cloud cover describes the whole 110 km tile, not our reservoirs, so we keep almost everything and decide per reservoir in the next step. `reservoir_water_quality.py` ```python def to_reflectance(dn: np.ndarray, item) -> np.ndarray: """L2A digital numbers to surface reflectance. Baseline 04.00+ (from 25 Jan 2022) adds a -1000 offset.""" offset = -1000 if float(item.properties["s2:processing_baseline"]) >= 4.0 else 0 return np.where(dn > 0, (dn.astype("float32") + offset) / 10000, np.nan) ``` ## Step 3: mask cloud, shadow and anything that is not open water For each scene we read the SCL band first. It is small and tells us whether the other bands are worth downloading. Cloud, cloud shadow, cirrus, snow and saturated pixels are masked. Then the mask grows by 100 m, because the edges of clouds and shadows are rarely classified well. A reservoir counts as clear only if at least half of its core pixels survive. Each band is read with a rasterio window and `out_shape`, which averages the 10 m bands to 20 m as it reads. Only the part of each COG we need is downloaded. Each asset URL is signed at read time with `planetary_computer.sign`, because the tokens expire during a long run. `reservoir_water_quality.py` ```python BAD_SCL = [0, 1, 3, 8, 9, 10, 11] # no data, saturated, cloud shadow, cloud (medium, high), cirrus, snow CLOUD_BUFFER_PX = 5 # grow cloud and shadow by 100 m: their edges are rarely classified well MIN_CLEAR = 0.5 # need at least half of a reservoir's core pixels to be clear NIR_MAX = 0.04 # water is dark in near-infrared; exposed mud, weed and boats are not def process(item) -> list[dict]: """Cloud-mask one scene and return the median water spectrum of each reservoir.""" p = item.properties base = {"date": item.datetime.strftime("%Y-%m-%d"), "item_id": item.id, "orbit": p["sat:relative_orbit"], "baseline": p["s2:processing_baseline"], "sun_zenith": round(p["s2:mean_solar_zenith"], 1)} scl = read_band(item, "SCL") bad = ndimage.binary_dilation(np.isin(scl, BAD_SCL[1:]), iterations=CLOUD_BUFFER_PX) | (scl == 0) clear_frac = {i: float((~bad[labels == i]).mean()) for i in range(1, len(RESERVOIRS) + 1)} rows = [dict(base, reservoir=n, clear_frac=round(clear_frac[i], 3)) for i, n in enumerate(RESERVOIRS, 1)] if max(clear_frac.values()) < MIN_CLEAR: return rows # too cloudy: skip the other bands entirely r = {b: to_reflectance(read_band(item, b), item) for b in BANDS} with np.errstate(divide="ignore", invalid="ignore"): mndwi = (r["B03"] - r["B11"]) / (r["B03"] + r["B11"]) water = ~bad & (r["B08"] < NIR_MAX) & np.isfinite(r["B11"]) # NIR-dark: no exposed bank, weed or boats for i, row in enumerate(rows, 1): core, w = labels == i, water & (labels == i) clear = core & ~bad row.update(n_pixels=int(w.sum()), water_frac=round(float(w.sum() / core.sum()), 3), scl_water_frac=round(float((scl[clear] == 6).mean()), 3) if clear.any() else np.nan, mndwi_pos_frac=round(float((mndwi[clear] > 0).mean()), 3) if clear.any() else np.nan) if clear_frac[i] >= MIN_CLEAR and w.sum() >= 20: row.update({f"rho_{b}": round(float(np.median(r[b][w])), 5) for b in BANDS}) return rows ``` We store the **median spectrum** of each reservoir’s water pixels and compute the indices from it later. A ratio of two medians is far more stable than the median of per-pixel ratios when the reflectances are close to zero. Running in 8 threads, all 821 scenes took about 90 seconds, and the results are cached in a CSV so the analysis can be rerun without downloading anything. ### Why MNDWI fails on peaty water The usual strict water test is MNDWI = (B3 − B11) / (B3 + B11) above 0 or 0.3, as in the [Lake Mead tutorial](https://www.technicalgis.com/tutorials/sentinel-2-stac-lake-mead/). We tried it first, and it rejected almost all the water. Peat-stained water is so dark that its green reflectance is about 1%, the same as the residual SWIR signal left by the atmospheric correction: `Output — why MNDWI fails on peaty water` ```text rho_B03 rho_B11 mndwi scl_water_frac mndwi_pos_frac reservoir Howden 0.0100 0.0121 -0.0942 0.9335 0.0925 Derwent 0.0100 0.0129 -0.0920 0.8250 0.0920 Ladybower 0.0110 0.0121 -0.1394 0.8570 0.0890 Torside 0.0101 0.0142 -0.1662 0.8540 0.0055 ``` On usable dates, median MNDWI is _negative_ over all four reservoirs, and fewer than 10% of clear core pixels have MNDWI above zero. The SCL classifier calls 83–93% of the same pixels water. So we combine three tests: the OSM outline shrunk by 40 m, no cloud or shadow in the SCL, and **low near-infrared** (B8 below 0.04). The NIR test removes exposed mud during drawdown, weed beds, boats and any shoreline pixels that remain. ## Step 4: compute the indices and control quality The indices come from the median water spectrum: `reservoir_water_quality.py` ```python def add_indices(df: pd.DataFrame) -> pd.DataFrame: """Relative indices computed from the median water spectrum of each reservoir and date.""" b = {k: df[f"rho_{k}"] for k in BANDS} df["ndci"] = (b["B05"] - b["B04"]) / (b["B05"] + b["B04"]) # chlorophyll-a proxy df["mci"] = b["B05"] - b["B04"] - (b["B06"] - b["B04"]) * (705 - 665) / (740 - 665) # red-edge peak height df["ndti"] = (b["B04"] - b["B03"]) / (b["B04"] + b["B03"]) # turbidity proxy df["red"] = b["B04"] # turbidity proxy (reflectance) df["green_blue"] = b["B03"] / b["B02"] # colour (CDOM, peat) proxy df["mndwi"] = (b["B03"] - b["B11"]) / (b["B03"] + b["B11"]) # kept only as a diagnostic return df ``` Then every reservoir-date gets a quality label. The first failing test wins: `reservoir_water_quality.py` ```python MAX_ZENITH = 70 # low winter sun: long paths through the atmosphere and weak water signal MIN_WATER = 0.6 # at least 60% of a reservoir's core pixels must pass the cloud and water tests MAX_SWIR = 0.03 # water absorbs SWIR almost completely: a bright median B11 means haze, glint or thin cloud MIN_RHO = 0.003 # B2-B5 darker than this: ratios are dominated by noise (Sen2Cor water error ~0.002-0.005) def qc_reason(row) -> str: if row["clear_frac"] < MIN_CLEAR or pd.isna(row["rho_B04"]): return "cloud or shadow" if row["sun_zenith"] > MAX_ZENITH: return "low sun" if row["water_frac"] < MIN_WATER: return "too few water pixels" if row["rho_B11"] > MAX_SWIR: return "haze or glint" if (row[BANDS_RHO[:-1]] <= 0).any(): return "negative reflectance" if (row[["rho_B02", "rho_B03", "rho_B04", "rho_B05"]] < MIN_RHO).any(): return "signal too dark" return "ok" ``` `Output — quality control` ```text reservoir Howden Derwent Ladybower Torside qc cloud or shadow 663 684 676 666 haze or glint 3 2 2 5 low sun 35 26 30 25 negative reflectance 10 10 5 9 ok 80 80 86 80 signal too dark 13 9 7 11 too few water pixels 17 10 15 25 326 usable reservoir-dates on 122 dates, 24 Feb 2019 to 22 Sep 2026 ``` This is the reality of optical monitoring in the Pennines: **about 82% of reservoir-dates are lost to cloud or shadow**, and only about one scene in seven is usable. Winter is almost empty, because for most of November to January the sun is less than 20° above the horizon at the time of the overpass, and the scenes are usually cloudy anyway. The record is uneven too, from 7 usable dates in 2024 to 23 in 2026. The negative-reflectance and too-dark tests matter for NDCI. When B4 and B5 are both about 0.003, an error of 0.001 in either changes NDCI by about 0.15, which is more than the whole seasonal signal. ## Step 5: build the time series and the monthly climatology ![Four stacked time-series panels, one per reservoir, of NDCI from 2019 to 2026, with summer months shaded yellow. Each reservoir's dots scatter mostly between −0.05 and 0.12 around a dashed median line near 0.03 to 0.05. Circled outliers mark high values, several of them on the same dates in more than one reservoir, such as May 2023 and April 2026.](https://www.technicalgis.com/_astro/05-ndci-timeseries.CIvHLhjL_M2K1S.webp) _NDCI for every usable date. The dashed line is each reservoir's median, yellow bands are June–August, and circles mark dates with a robust z-score above 2 (Step 6). The gaps are cloud._ The individual dates are noisy, so it is easier to see the seasonal pattern in a monthly climatology: the median and interquartile range of each calendar month over all years. `Output — median NDCI by season (n in brackets)` ```text spring summer autumn Howden 0.035 (42) 0.045 (29) 0.071 (7) Derwent 0.044 (44) 0.038 (30) 0.038 (4) Ladybower 0.054 (48) 0.036 (32) 0.114 (3) Torside 0.042 (43) 0.007 (30) 0.042 (6) ``` `Output — whole-record medians` ```text ndci mci ndti red green_blue n reservoir Howden 0.0412 0.0004 -0.0308 0.0094 0.8301 80 Derwent 0.0421 0.0002 -0.0395 0.0098 0.8525 80 Ladybower 0.0490 0.0003 -0.0355 0.0096 0.8828 86 Torside 0.0287 0.0007 -0.0078 0.0098 0.8245 80 ``` ![Three line charts of monthly medians from March to September with shaded interquartile ranges, one line per reservoir. NDCI: the three Upper Derwent reservoirs rise in April, dip in June and July, and rise again in August and September; Torside stays lower and dips to about zero in August. NDTI: values fall from March to June and vary in late summer. B3/B2: a peak in April and another in July for all reservoirs.](https://www.technicalgis.com/_astro/06-monthly-climatology.BPTg5LpG_hUHWk.webp) _Monthly median (line) and interquartile range (shading) of NDCI, NDTI and the B3/B2 colour ratio, 2019–2026, for months with at least three usable dates._ What the indices suggest, with the caveats in mind: - **Ladybower has the highest median NDCI** (0.049 over 86 dates). It is the largest and lowest of the three Upper Derwent reservoirs and receives the water that has passed through Howden and Derwent. - **The Upper Derwent reservoirs show two NDCI peaks**: one in April (Ladybower’s April median is 0.059) and a larger one in August and September (Ladybower 0.064 in August, 0.114 in September, though from only 3 September dates). That pattern is consistent with the classic temperate-lake cycle of a spring diatom bloom and a late-summer bloom, with a clear-water phase in between. - **Torside behaves differently.** It has the lowest NDCI (0.029) and the highest turbidity proxy (NDTI −0.008, against −0.031 to −0.040 for the others), and its summer NDCI is close to zero. On the map date its north-eastern arm had exposed bed (Step 7). Shallow, drawn-down margins can stir up sediment, which would fit a more turbid, less algal signal, but the satellite alone cannot prove the cause. - **The colour ratio B3/B2** peaks in April and July. In theory, more peat colour absorbs more blue and raises B3/B2. In practice, Sen2Cor tends to over-correct the blue band over dark water, so median B3/B2 is below 1 even though the water is brown. Treat it as the weakest of the indices. ## Step 6: rank the dates most worth investigating A water manager does not need a perfect concentration every five days. They need to know **which dates look unusual**, so someone can check the sampling records or go out with a sonde. We compare each date with the same reservoir’s typical value for that calendar month, using a **robust z-score** (median and median absolute deviation). Months with fewer than five observations fall back to the whole record: `reservoir_water_quality.py` ```python good["month"] = good["date"].dt.month grp = good.groupby(["reservoir", "month"])["ndci"] good["clim_med"] = grp.transform("median") good["clim_mad"] = grp.transform(lambda x: (x - x.median()).abs().median()) good["clim_n"] = grp.transform("size") # Months with fewer than five observations fall back to the reservoir's whole-record statistics. whole = good.groupby("reservoir")["ndci"] few = good["clim_n"] < 5 good.loc[few, "clim_med"] = whole.transform("median")[few] good.loc[few, "clim_mad"] = whole.transform(lambda x: (x - x.median()).abs().median())[few] good["ndci_z"] = (good["ndci"] - good["clim_med"]) / (1.4826 * good["clim_mad"]) # If most reservoirs spike on the same day, suspect the atmosphere (haze, aerosol model) before the water. high = good[good["ndci_z"] > 2] good["n_high_same_day"] = good["date"].map(high.groupby("date").size()).fillna(0).astype(int) good["flag"] = np.where(good["n_high_same_day"] >= 3, "check atmosphere (≥3 reservoirs high)", np.where(good["ndci_z"] > 2, "investigate", "")) ``` The last check is the most useful one. Separate reservoirs rarely bloom on exactly the same day, but an error in the atmospheric correction affects the whole scene. So when three or more reservoirs are high on the same date, the atmosphere is the first suspect. `Output — dates most worth investigating (NDCI robust z > 2)` ```text date reservoir ndci ndci_z ndti green_blue sun_zenith flag 2023-05-21 Derwent 0.118 4.873 -0.088 0.947 34.8 check atmosphere (≥3 reservoirs high) 2025-04-10 Howden 0.115 3.922 -0.011 1.093 46.8 investigate 2019-09-17 Ladybower 0.205 3.484 0.000 0.660 51.9 investigate 2025-05-13 Derwent 0.093 3.331 0.025 1.167 35.9 investigate 2025-04-30 Derwent 0.131 3.261 -0.049 0.970 39.9 investigate 2020-08-12 Derwent 0.116 3.181 -0.121 0.866 40.0 investigate 2026-04-08 Howden 0.099 3.038 0.045 1.455 47.2 check atmosphere (≥3 reservoirs high) 2026-04-08 Torside 0.139 2.997 0.087 1.738 47.2 check atmosphere (≥3 reservoirs high) 2026-04-08 Ladybower 0.122 2.960 0.000 1.418 47.2 check atmosphere (≥3 reservoirs high) 2023-05-21 Howden 0.109 2.604 -0.106 0.932 34.8 check atmosphere (≥3 reservoirs high) 2026-03-04 Howden 0.098 2.602 0.150 0.515 61.1 investigate 2025-07-12 Howden 0.143 2.576 0.069 1.373 33.0 investigate ``` Spring 2025 stands out: Howden on 10 April and Derwent on 30 April and 13 May were all well above their usual April and May values, while the other reservoirs were not. That fits a real spring bloom. The two scene-wide spikes, on 21 May 2023 and 8 April 2026, are more likely to be haze or an aerosol error, though a regional bloom after a warm spell is possible. In both cases the next step is the same: look at the true-colour image and check it against sampling records. ## Step 7: map NDCI and turbidity for one clear summer date For the maps, the script chooses the summer date on which all four reservoirs had the most clear water and the least SWIR haze: `Output — map date` ```text 12 July 2026 S2B_MSIL2A_20260712T112109_R037_T30UWE_20260712T133747 sun zenith 33.0° Howden NDCI p10/median/p90 0.019 0.052 0.112 | NDTI median 0.028 Derwent NDCI p10/median/p90 0.023 0.059 0.110 | NDTI median 0.009 Ladybower NDCI p10/median/p90 0.023 0.056 0.100 | NDTI median 0.011 Torside NDCI p10/median/p90 0.017 0.041 0.070 | NDTI median 0.022 ``` ![Two true-colour panels at 10 m. Left: the three Upper Derwent reservoirs as a narrow, nearly black ribbon through green, wooded valleys, each outlined in white with a yellow dashed inner line 40 m inside the shore. Right: Torside Reservoir in Longdendale, also nearly black, with a pale strip of exposed shore along its north-eastern arm.](https://www.technicalgis.com/_astro/02-true-colour.CThHLgNX_ZwJdok.webp) _True colour on 12 July 2026. The dashed yellow line is the 40 m inset core used for the statistics. The water is almost black, typical of peat-stained upland reservoirs; the pale strip at Torside's north-eastern end is exposed bed during drawdown, which the NIR test removes._ Per-pixel ratios of such dark water are noisy, so for mapping only we average the bands over 3 × 3 water pixels (60 m) before computing the indices: `reservoir_water_quality.py` ```python def smooth(a: np.ndarray, mask: np.ndarray) -> np.ndarray: """3 x 3 mean over water pixels only (60 m), to tame pixel noise in very dark water.""" v = ndimage.uniform_filter(np.where(mask, a, 0), 3) w = ndimage.uniform_filter(mask.astype("float32"), 3) with np.errstate(divide="ignore", invalid="ignore"): return np.where(mask, v / w, np.nan) s = {b: smooth(r[b], water) for b in ("B03", "B04", "B05")} with np.errstate(divide="ignore", invalid="ignore"): ndci_map = (s["B05"] - s["B04"]) / (s["B05"] + s["B04"]) ndti_map = (s["B04"] - s["B03"]) / (s["B04"] + s["B03"]) ``` ![NDCI map over a greyscale background, coloured with the Spectral palette from blue (low) through yellow to red (high). The three Upper Derwent reservoirs are a speckled mix of green and yellow with a few red pixels at the edges. In Torside, the open water is mostly blue-green, with a band of dark red along the southern shore and at the north-eastern tip.](https://www.technicalgis.com/_astro/03-ndci-map.DRDuEwi3_Z1Lg1a1.webp) _NDCI on 12 July 2026, smoothed to 60 m. Values are relative. The red band along Torside's southern shore is next to steep, dark, wooded slopes and is most likely an adjacency and shadow effect rather than algae._ ![NDTI map over a greyscale background, coloured from teal (low) through white to brown (high). Howden's main basin and most of Torside are light to mid brown; Derwent and Ladybower are mostly white to pale teal, with brown patches at Ladybower's dam end.](https://www.technicalgis.com/_astro/04-turbidity-map.Ci1inS_-_10vcD.webp) _The NDTI turbidity proxy on the same date. Howden and Torside are relatively turbid compared with Derwent and Ladybower; brighter, browner water is often found where inflowing streams carry sediment and near shallow, drawn-down margins._ The maps show why the inset matters. Even 40 m from the shore, the pixels next to steep wooded slopes have higher NDCI than the open water. The **adjacency effect** (light scattered from bright land into the view of dark water) and terrain shadow both add red-edge signal that is not algae. For a reservoir only 300–400 m wide, that affects a large share of the surface. ## Step 8: export a tidy CSV per reservoir Every scene is kept in the export, with a `qc` column, so you can see why a date was rejected and change the thresholds later: `reservoir_water_quality.py` ```python keep = ["date", "reservoir", "item_id", "orbit", "baseline", "sun_zenith", "clear_frac", "water_frac", "n_pixels", *BANDS_RHO, "ndci", "mci", "ndti", "red", "green_blue", "mndwi", "qc"] tidy = obs[keep].merge(good[["item_id", "reservoir", "ndci_z", "flag"]], on=["item_id", "reservoir"], how="left") for name in RESERVOIRS: f = OUT / f"{name.lower()}_sentinel2_water_quality.csv" t = tidy[tidy["reservoir"] == name] t.round(5).to_csv(f, index=False, date_format="%Y-%m-%d") ``` `Output — tidy CSV exports` ```text howden_sentinel2_water_quality.csv 821 rows, 80 usable derwent_sentinel2_water_quality.csv 821 rows, 80 usable ladybower_sentinel2_water_quality.csv 821 rows, 86 usable torside_sentinel2_water_quality.csv 821 rows, 80 usable ``` Each row is one scene and one reservoir, with the clear and water fractions, the seven median reflectances, the indices, the quality label and, for usable dates, the robust z-score and flag. The files load straight into pandas, R, Excel or a dashboard. ## Step 9: check against Environment Agency samples Most routine monitoring of these reservoirs is done by the water companies: Severn Trent for the Upper Derwent and United Utilities for Longdendale. Their operational data is not published. However, the **Environment Agency Water Quality Archive** is public and has an open API. Searching it found sampling points at all four reservoirs, with chlorophyll-a and colour: `reservoir_water_quality.py` ```python EA = "https://environment.data.gov.uk/water-quality" POINTS = {"MD-49714900": "Howden", "MD-49714220": "Derwent", "MD-LSN0395": "Ladybower", "NW-88021743": "Torside"} DETERMINANDS = {"7887": "chlorophyll_ugl", "0072": "colour_hazen"} resp = requests.get(f"{EA}/sampling-point/{point}/observation", headers={"Accept": "text/csv"}, params={"dateFrom": "2019-01-01", "dateTo": "2026-10-14", "determinand": code, "limit": 250, "skip": skip}, timeout=60) ``` `Output — Environment Agency samples (2019 onwards)` ```text reservoir var n min max Derwent chlorophyll_ugl 47 2019-01-16 2026-08-13 Derwent colour_hazen 19 2019-01-16 2022-07-19 Howden chlorophyll_ugl 11 2019-01-16 2020-03-19 Howden colour_hazen 11 2019-01-16 2020-03-19 Ladybower chlorophyll_ugl 5 2026-04-20 2026-08-05 Ladybower colour_hazen 5 2026-04-20 2026-08-05 Torside chlorophyll_ugl 15 2019-02-11 2020-03-18 Torside colour_hazen 15 2019-02-11 2020-03-18 ``` We pair each sample with the nearest usable Sentinel-2 date within ±3 days. Values below the detection limit, such as `<0.5`, are kept at the limit and marked: `Output — satellite vs sample matches (within ±3 days)` ```text reservoir sample_date s2_date value censored ndci mci Howden 2019-02-21 2019-02-24 0.50 True 0.0279 0.0023 Howden 2019-03-25 2019-03-28 0.60 False 0.0340 0.0012 Howden 2020-03-19 2020-03-22 0.50 True 0.0680 0.0006 Derwent 2019-02-21 2019-02-24 0.50 True 0.0333 0.0025 Derwent 2019-03-25 2019-03-28 0.50 True 0.0299 0.0011 Derwent 2021-06-18 2021-06-15 3.00 False 0.0451 -0.0004 Derwent 2022-04-20 2022-04-21 0.76 False 0.0352 0.0002 Derwent 2022-06-23 2022-06-20 8.70 False -0.0145 0.0013 Derwent 2022-08-10 2022-08-12 1.80 False 0.0548 0.0017 Derwent 2024-03-21 2024-03-24 0.59 False 0.0000 0.0004 Derwent 2025-05-15 2025-05-13 1.70 False 0.0925 0.0001 Derwent 2025-06-17 2025-06-17 1.30 False 0.0571 0.0004 Derwent 2026-04-21 2026-04-23 1.40 False 0.0704 -0.0005 Ladybower 2026-04-20 2026-04-23 1.80 False 0.0714 -0.0008 Ladybower 2026-07-10 2026-07-09 4.50 False 0.0159 0.0010 Ladybower 2026-08-05 2026-08-03 9.00 False -0.0435 0.0013 Torside 2019-04-24 2019-04-22 2.70 False 0.0649 -0.0000 Torside 2019-08-28 2019-08-25 5.00 False 0.0687 -0.0001 chlorophyll vs NDCI n = 18 Spearman rho = -0.03 p = 0.90 chlorophyll vs MCI n = 18 Spearman rho = -0.27 p = 0.28 colour vs B3/B2 n = 12 Spearman rho = -0.04 p = 0.91 ``` ![Two scatter plots. Left: Environment Agency chlorophyll-a from 0.5 to 9 µg/L against Sentinel-2 NDCI from −0.04 to 0.09, with no visible trend; the two highest chlorophyll values have the lowest NDCI. Right: filtered colour from about 40 to 140 Hazen against B3/B2 from 0.65 to 1.15, also with no trend.](https://www.technicalgis.com/_astro/07-in-situ-check.BtauiH80_26gVyw.webp) _Sentinel-2 indices against Environment Agency samples taken within three days. Neither NDCI nor the colour ratio tracks the measured values (Spearman ρ = −0.03 and −0.04). Open circles are samples below the detection limit._ **The honest result: at these reservoirs, Sentinel-2 NDCI from Sen2Cor reflectance does not track measured chlorophyll.** The two highest samples, 8.7 µg/L at Derwent on 23 June 2022 and 9.0 µg/L at Ladybower on 5 August 2026, came with some of the _lowest_ NDCI values. There are several reasons: - **The chlorophyll is low.** Every sample is between 0.5 and 9 µg/L. NDCI was developed for productive, turbid water, and its red-edge peak only becomes clear above roughly 10–20 µg/L. Below that, peat colour and residual atmospheric error control the signal. - **The water is very dark.** With B4 and B5 near 0.005–0.01, an atmospheric-correction error of 0.002 changes NDCI by more than the whole seasonal range. - **Small, uneven samples.** There are only 18 pairs, the samples are taken at the dam or a jetty rather than the open water, and ±3 days is long enough for a bloom to move or die back. That does not make the workflow useless. It shows what the indices can and cannot do. With this level of correction, use them to find **unusual dates and spatial patterns** that deserve a closer look, and keep the time series as a consistent record you can calibrate later. Do not report them as concentrations. ## Pitfalls and limitations - **Relative, not calibrated.** NDCI, MCI, NDTI and B3/B2 are proxies. Turning them into µg/L or NTU needs a local regression against matched samples, and the result above shows that would not work yet with Sen2Cor reflectance. - **Atmospheric correction.** Sen2Cor is designed for land. Over dark water it often over-corrects the blue band and leaves errors of a few thousandths in reflectance, which is a large share of the water signal. The standard next steps are **ACOLITE** (dark spectrum fitting) or **C2RCC** in ESA SNAP, both made for inland and coastal water, run on Level-1C data. - **Cloud and shadow.** About 82% of reservoir-dates were lost. Shadows from small cumulus clouds that SCL misses are a common cause of false dark or bright water; the 100 m cloud buffer helps but is not perfect. - **Adjacency and terrain shadow.** In narrow valleys, scattered light from bright land and shadows from steep slopes reach well beyond 40 m, as the NDCI map at Torside shows. A 60–100 m inset or an adjacency correction (ACOLITE has one) is safer for science-grade work. - **Sunglint.** In summer the sun is high and wind-roughened water can reflect sunlight straight towards the sensor. The SWIR test (median B11 above 0.03) catches strong glint, but weak glint still lifts all bands. - **Shallow water and drawdown.** In shallow margins and during low water, the bed contributes to the signal. The NIR test removes exposed mud, but not water a few tens of centimetres deep over a bright bed. - **Peat colour (CDOM).** Brown water absorbs blue and green strongly and lowers the signal in every band, which leaves less room for the chlorophyll signal. The colour ratio is unreliable while the blue band is poorly corrected. - **Low winter sun.** Above a solar zenith angle of 70° the atmospheric path is long and the signal weak, so winter has almost no data. ## Frequently asked questions ### Can Sentinel-2 measure chlorophyll? Not directly. Sentinel-2 measures reflected light, and its red-edge band B5 (705 nm) responds to chlorophyll, so indices such as NDCI and MCI give a chlorophyll _proxy_. In productive, turbid lakes with tens of µg/L, published algorithms calibrated with field samples can estimate concentrations reasonably well. In clear or peat-stained water with low chlorophyll, such as these Peak District reservoirs, the signal is too weak for reliable estimates without water-specific atmospheric correction and local calibration. ### How do you monitor algal blooms from satellite? Use a sensor with a red-edge band (Sentinel-2 MSI or Sentinel-3 OLCI), mask cloud and shoreline pixels, compute a chlorophyll index such as NDCI, MCI or, with OLCI, the cyanobacteria index, and compare each date with that lake’s own seasonal normal. Treat large anomalies as alerts to investigate with sampling, as in Step 6. For larger lakes, the Copernicus Land Monitoring Service publishes ready-made lake water-quality products such as turbidity and trophic state. ### What is NDCI? The Normalised Difference Chlorophyll Index, (B5 − B4) / (B5 + B4) for Sentinel-2, was published by Mishra and Mishra in 2012. Chlorophyll absorbs red light (B4, 665 nm), while algal cells reflect strongly just beyond the absorption, near 705 nm (B5). More algae push NDCI up. Values are relative and depend on the atmospheric correction. ### Which satellite is best for monitoring small reservoirs? Sentinel-2 is the best free option for reservoirs smaller than a few square kilometres, with 10–20 m pixels, a red-edge band and a 2–5 day revisit. Landsat 8 and 9 have 30 m pixels and no red-edge band, so they suit turbidity and temperature rather than chlorophyll. Sentinel-3 OLCI has excellent water bands but 300 m pixels, which is too coarse for lakes narrower than about 1 km. Commercial sensors such as PlanetScope have more frequent images but weaker spectral and radiometric quality. ### Why is my reservoir’s MNDWI negative? Very dark water, such as peat-stained upland reservoirs, reflects only about 1% of green light. After atmospheric correction, the residual signal in the SWIR band can be just as large, so MNDWI is close to zero or negative. Use a known outline shrunk inwards, the SCL cloud mask and a near-infrared threshold instead, as in Step 3. ### Is there public water quality data for UK reservoirs? Partly. The Environment Agency Water Quality Archive is open and includes some reservoir samples, as used in Step 9. Most operational monitoring of drinking-water reservoirs is done by the water companies and is not published. The UK Centre for Ecology & Hydrology runs long-term monitoring at some lakes, mainly in Cumbria and Scotland. ## Next steps - Process Level-1C scenes with **ACOLITE** or **C2RCC** and rerun the comparison with the Environment Agency samples; a water-specific correction is the single biggest improvement. - Add **Sentinel-3 OLCI** for a large lake nearby, or switch to a larger lake such as Rutland Water, where 300 m pixels are usable. - Add **Landsat 8 and 9** to fill gaps in the turbidity record, after harmonising the bands. - Map the reservoirs’ water area through the year, as in [Measuring Lake Mead’s shrinking shoreline with Sentinel-2 and STAC](https://www.technicalgis.com/tutorials/sentinel-2-stac-lake-mead/) or our earlier [Lake Mead drought post](https://www.technicalgis.com/2022/08/11/monitoring-lake-mead-drought/). - When cloud blocks every optical image, track water extent with radar instead: [Sentinel-1 SAR flood mapping in Python](https://www.technicalgis.com/tutorials/sentinel-1-sar-flood-mapping/). - Get the inflowing streams and catchments with [OSMnx](https://www.technicalgis.com/tutorials/street-networks-osmnx/) or OS Open Rivers, and relate the colour signal to peat cover upstream. ## Data and licences - **Sentinel-2 L2A**: contains modified Copernicus Sentinel data 2019–2026, processed by ESA. Accessed as Cloud Optimised GeoTIFFs from the [Microsoft Planetary Computer](https://planetarycomputer.microsoft.com/dataset/sentinel-2-l2a) on 15 October 2026. - **Reservoir outlines**: © OpenStreetMap contributors, available under the Open Database Licence (ODbL). - **Water quality samples**: Environment Agency Water Quality Archive, © Environment Agency, licensed under the Open Government Licence v3.0. Accessed on 15 October 2026. - **Country outlines**: Natural Earth, public domain. - References: Mishra, S. and Mishra, D. R. (2012), “Normalized difference chlorophyll index”, _Remote Sensing of Environment_ 117; Gower, J. et al. (2005), “Detection of intense plankton blooms using the 709 nm band of the MERIS imaging spectrometer”, _International Journal of Remote Sensing_ 26(9); Lacaux, J. P. et al. (2007), _Remote Sensing of Environment_ 106(1). - Code: [`reservoir_water_quality.py`](https://www.technicalgis.com/code/reservoir-water-quality-sentinel-2/), tested with Python 3.12.0, pystac-client 0.9.0, planetary-computer 1.0.0, rasterio 1.5.2, osmnx 2.1.1 and pandas 3.0.6. ## Key terms - **NDCI**: Normalised Difference Chlorophyll Index, (B5 − B4) / (B5 + B4). It compares the red-edge reflectance peak caused by algae with the red absorption by chlorophyll (Mishra and Mishra, 2012). - **NDTI**: Normalised Difference Turbidity Index, (B4 − B3) / (B4 + B3). Suspended sediment raises red reflectance relative to green (Lacaux et al., 2007). - **CDOM**: Coloured dissolved organic matter: dissolved humic substances, mostly from peat in the uplands, that stain water brown and absorb blue light strongly. - **SCL**: The Sentinel-2 L2A scene classification layer: a 20 m band that labels each pixel as cloud, shadow, water, vegetation, snow and so on. - **Adjacency effect**: Light scattered in the atmosphere from bright land next to a lake brightens the dark water pixels near the shore, which biases water indices. - **Robust z-score**: How far a value is from the median, in units of 1.4826 × the median absolute deviation. Unlike the usual z-score, it is not inflated by the outliers it is trying to find. ## Data for this tutorial - [Microsoft Planetary Computer](https://www.technicalgis.com/data-sources/#src-microsoft-planetary-computer) · Sentinel-2 L2A as Cloud Optimised GeoTIFFs via STAC - [Copernicus Data Space Ecosystem](https://www.technicalgis.com/data-sources/#src-copernicus-data-space-ecosystem) · Sentinel-2 L1C for ACOLITE or C2RCC processing - [OpenStreetMap](https://www.technicalgis.com/data-sources/#src-openstreetmap) · Reservoir outlines (water=reservoir) - [Natural Earth](https://www.technicalgis.com/data-sources/#src-natural-earth) · Country outlines for the locator inset Licences and formats for each source are listed in the [GIS data sources catalogue](https://www.technicalgis.com/data-sources/). --- Source: https://www.technicalgis.com/tutorials/reservoir-water-quality-sentinel-2/