How to Detect New Forest Clearing with Sentinel-2 in Python
On this page
Problem statement
Annual forest-loss products arrive months after the year they describe. Monitoring needs something faster: compare a recent satellite image with one from before and flag forest that has turned into bare soil or pasture. Sentinel-2 revisits every five days at 10 m, is free, and is served as cloud-optimised GeoTIFFs that can be read by window. The basic detector — a drop in NDVI between two cloud-free dry-season images, inside a forest mask — takes a page of code.
Its difficulties are the ones every change detector has: clouds and their shadows, seasonal changes that look like clearing, and a threshold that trades missed clearings against false alarms. Checking the result against an independent loss map tells you where on that trade-off you are.
Quick answer
Pick two clear dry-season scenes a year apart, compute NDVI for both inside the forest mask, and flag large drops:
import numpy as np, rasterio, pystac_client
from rasterio.windows import from_bounds
from rasterio.warp import transform_bounds
W, S, E, N = -62.95, -10.75, -62.80, -10.60 # ~16 × 16 km inside the Hansen window
cat = pystac_client.Client.open("https://earth-search.aws.element84.com/v1")
def clearest(date_range):
items = cat.search(collections=["sentinel-2-l2a"], bbox=[W, S, E, N], datetime=date_range,
query={"eo:cloud_cover": {"lt": 5}}).item_collection()
return min(items, key=lambda i: i.properties["eo:cloud_cover"])
def ndvi(item):
out = {}
for band in ("red", "nir", "scl"):
with rasterio.open(item.assets[band].href) as src:
win = from_bounds(*transform_bounds("EPSG:4326", src.crs, W, S, E, N), src.transform)
scale = 2 if band == "scl" else 1 # SCL is 20 m
out[band] = src.read(1, window=win, out_shape=(int(win.height * scale),
int(win.width * scale))).astype("float32")
if band == "red":
tf, crs = src.window_transform(win), src.crs
bad = np.isin(out["scl"], [0, 1, 3, 8, 9, 10]) # nodata, saturated, shadow, cloud
v = (out["nir"] - out["red"]) / (out["nir"] + out["red"] + 1e-6)
return np.where(bad, np.nan, v), tf, crs
before = clearest("2023-07-01/2023-08-31")
after = clearest("2024-07-01/2024-08-31")
n0, tf, crs = ndvi(before)
n1, _, _ = ndvi(after)
print(f"before: {before.id} ({before.properties['eo:cloud_cover']:.1f}% cloud)")
print(f"after: {after.id} ({after.properties['eo:cloud_cover']:.1f}% cloud)")
forest_like = n0 > 0.75 # dense green in the dry season
drop = n0 - n1
cleared = forest_like & (drop > 0.35)
ha = 0.01 # a 10 m pixel is 0.01 ha
print(f"forest-like pixels before: {forest_like.sum() * ha:,.0f} ha")
print(f"flagged as cleared (NDVI drop > 0.35): {cleared.sum() * ha:,.0f} ha")before: S2A_20LMP_20230704_0_L2A (0.0% cloud)
after: S2A_20LMP_20240718_0_L2A (0.0% cloud)
forest-like pixels before: 10,852 ha
flagged as cleared (NDVI drop > 0.35): 822 ha
Two images, one subtraction, a threshold. The flagged area is a candidate list — the next step is checking it against something independent.
Dry season, masks, a forest baseline, a threshold
1. Use the dry season
In southern Amazonia the dry season (June–September) gives the clearest images and is when clearing and burning happen. Comparing the same season in consecutive years removes most phenology.
2. Mask clouds and shadows
Sentinel-2 L2A's scene classification layer (SCL) marks clouds, cirrus and cloud shadows. Shadows look like dark bare soil in NDVI terms; masking them prevents false clearings.
3. Build a forest mask from the earlier image
Dense, green vegetation in the dry season — NDVI above about 0.75 — is a practical forest proxy for the baseline. A proper forest map (Hansen, national products) is better where available.
4. Threshold the drop
Clearing takes NDVI from 0.8+ to bare-soil or pasture values (0.2–0.5). A drop of 0.3–0.4 separates clearing from normal variation; the threshold sets the balance between missed and false detections.
5. Clean the result
Remove isolated pixels (a minimum mapping unit of a few pixels) and flag detections next to masked clouds as less certain.
6. Validate
Compare with an independent loss map for the same period, or with a sample checked visually on imagery.
New clearing against Hansen's loss
The data used on this page
The Rondônia window used throughout the forest guides (63.0°–62.5°W, 10.8°–10.3°S): Hansen Global Forest Change v1.12 — tree cover in 2000, the year of loss, and the data mask — read from Google's storage into a three-band rondonia_gfc.tif, as in the Hansen reading guide.
import os
import rasterio
from rasterio.windows import from_bounds
W, S, E, N = -63.0, -10.8, -62.5, -10.3
V = "GFC-2024-v1.12"
if not os.path.exists("rondonia_gfc.tif"):
bands = []
for name in ("treecover2000", "lossyear", "datamask"):
url = f"https://storage.googleapis.com/earthenginepartners-hansen/{V}/Hansen_{V}_{name}_10S_070W.tif"
with rasterio.open(url) as src:
window = from_bounds(W, S, E, N, src.transform)
bands.append(src.read(1, window=window))
profile = src.profile | {"count": 3, "width": bands[0].shape[1], "height": bands[0].shape[0],
"transform": src.window_transform(window), "compress": "deflate",
"driver": "GTiff", "tiled": False}
profile.pop("blockxsize", None)
profile.pop("blockysize", None)
with rasterio.open("rondonia_gfc.tif", "w", **profile) as dst:
for i, band in enumerate(bands, start=1):
dst.write(band, i)
with rasterio.open("rondonia_gfc.tif") as src:
print(src.count, src.shape, src.crs, src.res)3 (2000, 2000) EPSG:4326 (0.00025, 0.00025)
Example 1 — compare with Hansen's 2024 loss
from rasterio.warp import reproject, Resampling
with rasterio.open("rondonia_gfc.tif") as src:
ly = np.zeros(n0.shape, dtype="uint8")
reproject(rasterio.band(src, 2), ly, dst_transform=tf, dst_crs=crs,
resampling=Resampling.nearest)
hansen_2024 = ly == 24
valid = forest_like & ~np.isnan(drop)
tp = (cleared & hansen_2024 & valid).sum()
print(f"Hansen 2024 loss inside the forest mask: {(hansen_2024 & valid).sum() * ha:,.0f} ha")
print(f"detected by NDVI drop: {tp * ha:,.0f} ha "
f"({tp / max((hansen_2024 & valid).sum(), 1):.0%} of Hansen)")
print(f"NDVI detections Hansen has no 2024 loss for: "
f"{(cleared & ~hansen_2024 & valid).sum() * ha:,.0f} ha")Hansen 2024 loss inside the forest mask: 170 ha
detected by NDVI drop: 11 ha (6% of Hansen)
NDVI detections Hansen has no 2024 loss for: 812 ha
Agreement is poor in both directions, and the dates explain much of it. The two images are from early July 2023 and mid-July 2024, while Hansen's 2024 loss runs to December. Clearing and burning peak in August and September, after the second image, so most of Hansen's 2024 loss had not yet happened when it was taken. Clearing between July and December 2023 is in the NDVI difference but belongs to Hansen's 2023. And 2024 was a severe drought year in the Amazon, when browning canopy and burnt vegetation lower NDVI without a stand-replacing loss. A two-date detector measures change between two days, not loss in a calendar year: compare it with a reference over the same interval.
Example 2 — the threshold trade-off
for t in (0.2, 0.3, 0.35, 0.45, 0.55):
det = forest_like & (drop > t) & valid
hit = (det & hansen_2024).sum() / max((hansen_2024 & valid).sum(), 1)
extra = (det & ~hansen_2024).sum() * ha
print(f"drop > {t:.2f}: {det.sum() * ha:>7,.0f} ha flagged, {hit:5.0%} of Hansen 2024 found, "
f"{extra:>7,.0f} ha not in Hansen 2024")drop > 0.20: 2,133 ha flagged, 26% of Hansen 2024 found, 2,090 ha not in Hansen 2024
drop > 0.30: 1,212 ha flagged, 12% of Hansen 2024 found, 1,191 ha not in Hansen 2024
drop > 0.35: 822 ha flagged, 6% of Hansen 2024 found, 812 ha not in Hansen 2024
drop > 0.45: 196 ha flagged, 1% of Hansen 2024 found, 194 ha not in Hansen 2024
drop > 0.55: 10 ha flagged, 0% of Hansen 2024 found, 10 ha not in Hansen 2024
Lower thresholds find more of the reference loss and flag much more that the reference does not have; at no threshold does the detector find even a third of Hansen's 2024 loss here, because the dates do not line up. There is no correct threshold, only a choice about which error is cheaper — and no threshold fixes a mismatch in the period.
Example 3 — a minimum mapping unit
from scipy import ndimage
labels, n = ndimage.label(cleared)
sizes = np.bincount(labels.ravel())[1:]
keep = np.isin(labels, np.where(sizes >= 10)[0] + 1) # at least 0.1 ha
print(f"{n} detected patches; {(sizes >= 10).sum()} of at least 0.1 ha, "
f"covering {keep.sum() * ha:,.0f} ha of the {cleared.sum() * ha:,.0f} ha flagged")2003 detected patches; 741 of at least 0.1 ha, covering 762 ha of the 822 ha flagged
Most noise is in isolated pixels and tiny patches; real clearings are contiguous. A minimum patch size removes much of the noise for a small loss of detections.
Why NDVI drops, and two dates are fragile
Why NDVI drops at clearing
Healthy canopy reflects strongly in the near infrared and absorbs red light; bare soil and burned or dry pasture do neither. Clearing therefore moves NDVI from above 0.8 to 0.2–0.5 in a single step, a much larger change than seasonal variation in evergreen forest.
Why two dates are fragile
Any difference between two single images includes clouds missed by the mask, haze, sun angle and phenology. Operational systems such as Brazil's DETER use many images and time-series methods for exactly this reason. A two-date detector is a transparent starting point, not an alert system.
Why the reference disagrees too
Hansen's loss is also a model output, on a calendar year and at 30 m. Disagreement between two maps measures the uncertainty of both. A sample of visually checked points is the only way to estimate the accuracy of either.
Cloud-free scenes, logging and fire scars
- Cloud-free scenes are rare in the wet season; the method is a dry-season tool in the tropics.
- Selective logging removes individual trees and barely changes NDVI; it needs other methods.
- Fire scars drop NDVI like clearing and may recover; a later image separates them.
- Radar (Sentinel-1) sees through clouds and is used for year-round alerts.
- Harmonisation: processing baseline changes in Sentinel-2 L2A (2022) shifted reflectance offsets; use harmonised collections for long series.
Internal links
- How to calculate NDVI in Python — the index itself
- How to detect change between satellite images — change detection in general
- SAR change detection in Python — the cloud-proof alternative
- How to read Hansen Global Forest Change data — the reference layer
- How to mask clouds in Sentinel-2 — masking in depth
- Global Forest Change explained — what the reference measures
FAQ
How do I detect deforestation with Sentinel-2?
Compare NDVI between cloud-masked dry-season images a year apart, within a forest mask, and flag large drops.
What NDVI drop indicates clearing?
Around 0.3–0.4 from a forest baseline above 0.75; tune the threshold against reference data.
How do I mask clouds in Sentinel-2 L2A?
Use the scene classification layer: exclude cloud, cirrus, cloud shadow, saturated and no-data classes.
Why doesn't my detection match Hansen?
Different periods (dry season to dry season against calendar year), resolutions, masks and methods.
Can this detect selective logging?
Not reliably; canopy gaps from logging barely change NDVI at 10 m.
What about the rainy season?
Optical images are mostly cloudy then; radar-based methods fill the gap.