Avoided Deforestation Explained

Problem statement

A protected area kept 95% of its forest while the land around it kept 60%. Did protection save 35% of the forest? Only if, without protection, the protected land would have lost as much as its surroundings — and that is exactly what cannot be observed. Protected areas, indigenous territories and conservation projects are placed deliberately, often where clearing was unlikely anyway: far from roads, on poor soils, up steep slopes. Comparing them with everything else credits the policy with what the geography did.

Avoided deforestation is the difference between what happened and a counterfactual — what would have happened without the intervention. Estimating the counterfactual well is the whole problem, and it matters in practice because carbon credits are paid on it.

Quick answer

Three comparisons, in increasing credibility:

1. Inside vs outside — all protected land against all unprotected land. Confounded by placement.
2. Inside vs a boundary band — protected land against land just outside the boundary. Controls for regional conditions.
3. Matched comparison — protected pixels against unprotected pixels with similar covariates (distance to roads, slope, distance to past clearing). The standard in impact evaluation.

import numpy as np, osmnx as ox, rasterio
from rasterio import features
from rasterio.windows import from_bounds
from scipy import ndimage

V, TILE = "GFC-2024-v1.12", "10S_070W"
W, S, E, N = -63.5, -11.3, -62.8, -10.6
def layer(name):
    with rasterio.open(f"/vsicurl/https://storage.googleapis.com/earthenginepartners-hansen/"
                       f"{V}/Hansen_{V}_{name}_{TILE}.tif") as src:
        win = from_bounds(W, S, E, N, src.transform)
        return src.read(1, window=win), src.window_transform(win)
tc, tf = layer("treecover2000")
ly, _ = layer("lossyear")
shape = tc.shape
samp = (abs(tf.e) * 110_574, abs(tf.a) * 111_320 * np.cos(np.radians((S + N) / 2)))
px_ha = samp[0] * samp[1] / 1e4                                 # ~0.075 ha per pixel here

ox.settings.http_user_agent = "spatialworkflow-examples/1.0"
pa = ox.features_from_bbox((W, S, E, N), {"boundary": ["protected_area", "aboriginal_lands"]})
pa = pa[pa.geom_type.isin(["Polygon", "MultiPolygon"])]
inside = features.rasterize(((g, 1) for g in pa.geometry), out_shape=shape, transform=tf) == 1
roads = ox.features_from_bbox((W - 0.1, S - 0.1, E + 0.1, N + 0.1), {"highway": True})
roads = roads[roads.geom_type.isin(["LineString", "MultiLineString"])]
burn = features.rasterize(((g, 1) for g in roads.geometry), out_shape=shape, transform=tf,
                          all_touched=True)
d_road = ndimage.distance_transform_edt(burn == 0, sampling=samp) / 1000
d_out = ndimage.distance_transform_edt(~inside, sampling=samp) / 1000

forest = tc > 30
lost = ly > 0
naive = lost[forest & ~inside].mean() - lost[forest & inside].mean()
band = forest & ~inside & (d_out <= 5)
boundary = lost[band].mean() - lost[forest & inside & (ndimage.distance_transform_edt(inside, sampling=samp) / 1000 <= 5)].mean()
print(f"loss share inside {lost[forest & inside].mean():.1%}, outside {lost[forest & ~inside].mean():.1%}")
print(f"naive difference: {naive:.1%} points; boundary-band difference: {boundary:.1%} points")

loss share inside 1.2%, outside 60.3%
naive difference: 59.1% points; boundary-band difference: 61.1% points

Two comparisons, two estimates of the protection effect — and here the boundary band makes the gap larger, not smaller. Land just outside these boundaries has been cleared at least as heavily as unprotected land in general. That can be leakage, clearing pushed up to the edge, or boundaries that were drawn along an existing frontier; either way the band is not the like-for-like comparison it is meant to be.

Name the counterfactual, measure drivers, compare like with like

1. Name the counterfactual

"What would the protected land have lost without protection?" Everything else is an attempt to approximate it.

2. Measure the covariates that drive clearing

Distance to roads and towns, slope, soil, distance to earlier clearing, land tenure. These are what placement correlates with.

3. Compare like with like

Stratify or match protected and unprotected pixels on those covariates; compare loss within each stratum or pair.

4. Check balance

After matching, the covariate distributions of the two groups should be similar. If they are not, the comparison is still confounded.

5. Watch for leakage

Protection can push clearing next door. If the comparison land is next to the protected area, leakage inflates the apparent effect.

Stratifying by road distance and checking balance

The data used on this page

The Hansen/UMD Global Forest Change layers (version GFC-2024-v1.12, tile 10S 070W, Rondônia, Brazil), read by window straight from Google Cloud Storage, and protected areas, indigenous lands and roads from OpenStreetMap, fetched with OSMnx by the page's own code. The setup points OSMnx at an Overpass mirror and unpacks its answers as saved on 9 October 2026: OpenStreetMap is edited daily, and a live query would move the counts a little.

import io
import os
import osmnx as ox
import requests
import zipfile

# OpenStreetMap as Overpass returned it on 9 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-avoided-deforestation-explained-2026-10-09.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"
ox.settings.overpass_rate_limit = False          # its /status endpoint can hang
ox.settings.requests_timeout = 600
print("Overpass:", ox.settings.overpass_url)

Overpass: https://maps.mail.ru/osm/tools/overpass/api

Example 1 — stratify by distance to roads

bins = [0, 1, 2, 5, 10, 20, 100]
tot_w, eff = 0, 0
print(f"{'road distance':<14}{'inside':>9}{'outside':>9}{'difference':>12}{'protected ha':>14}")
for lo, hi in zip(bins[:-1], bins[1:]):
    s = forest & (d_road >= lo) & (d_road < hi)
    a, b = s & inside, s & ~inside
    if a.sum() > 5000 and b.sum() > 5000:
        diff = lost[b].mean() - lost[a].mean()
        tot_w += a.sum(); eff += diff * a.sum()
        print(f"{lo:>3}–{hi:<3} km{lost[a].mean():>11.1%}{lost[b].mean():>9.1%}{diff:>12.1%}"
              f"{a.sum() * px_ha:>14,.0f}")
print(f"stratified effect, weighted by protected area: {eff / tot_w:.1%} points")

road distance    inside  outside  difference  protected ha
  0–1   km      11.6%    60.7%       49.1%         4,240
  1–2   km       5.9%    60.9%       55.1%         9,724
  2–5   km       1.4%    55.1%       53.8%        37,825
stratified effect, weighted by protected area: 53.6% points

Within each road-distance class the difference between protected and unprotected land is smaller than the naive gap, because the naive gap mixed remote protected land with accessible unprotected land. But only three classes have enough of both to compare. Beyond 5 km from a road there is almost no unprotected forest left to compare with, so the stratified estimate speaks for the protected forest within 5 km of a road — about 52,000 ha, a tenth of the total. That lack of overlap is itself a finding: the rest of the protected forest has no counterpart outside.

Example 2 — covariate balance

print(f"median distance to a road: inside {np.median(d_road[forest & inside]):.1f} km, "
      f"outside {np.median(d_road[forest & ~inside]):.1f} km, "
      f"boundary band {np.median(d_road[band]):.1f} km")

median distance to a road: inside 17.6 km, outside 1.0 km, boundary band 1.0 km

The unprotected forest is a median 1 km from a road; the protected forest 17.6 km. The boundary band does not close that gap at all — it sits at 1 km too — which is why it did not reduce the estimate. Balance on the covariates is what makes a comparison credible, and it has to be shown, not assumed.

Example 3 — the size of the counterfactual

prot_forest_ha = (forest & inside).sum() * px_ha
for label, effect in (("naive", naive), ("boundary band", boundary), ("stratified", eff / tot_w)):
    print(f"{label:<14} avoided loss ≈ {effect * prot_forest_ha:>9,.0f} ha of the protected forest")

naive          avoided loss ≈   268,308 ha of the protected forest
boundary band  avoided loss ≈   277,544 ha of the protected forest
stratified     avoided loss ≈   243,714 ha of the protected forest

The same protected forest, three estimates of what protection saved, spanning 34,000 ha — and the lowest still applies an effect measured on the accessible tenth of the protected forest to the remote rest, which faced little threat. For a carbon project, the gap between estimates like these is the gap between credits claimed and credits earned.

Placement bias, matching and carbon markets

Why placement bias is large

Governments and communities establish protected areas where they can: land without competing claims, often remote and less suited to agriculture. Those are also the places least likely to be cleared. The naive comparison measures that difference in threat as much as the effect of protection.

Why matching helps but does not solve everything

Matching controls for the covariates you measured. If protection was placed according to something unmeasured — political influence, community presence — bias remains. Matching makes the comparison more credible; it does not make it an experiment.

Why this matters for carbon markets

Forest carbon credits are issued for avoided emissions: counterfactual loss minus actual loss. Inflated baselines — counterfactuals built from naive comparisons or exaggerated regional rates — produce credits for forest that was never at risk. Independent studies have found many projects' baselines far above what matched comparisons support.

Designation dates, spillover and synthetic controls

• Designation dates matter: compare from the year of protection, not from 2001.
• Indigenous lands often show strong effects that partly reflect active defence by communities, not only legal status.
• Spillover: protection can also reduce clearing nearby (positive spillover), biasing boundary comparisons the other way.
• Spatial autocorrelation makes pixel-level standard errors far too small; use blocks or areas as units.
• Synthetic control and difference-in-differences add time to the comparison when pre-protection data exist.

Internal links

• How to compare forest loss inside and outside protected areas — the full workflow
• How to measure deforestation near roads — the main covariate
• Protected area data explained — sources and designations
• Deforestation rates explained — rates for the counterfactual
• Spatial autocorrelation explained — why pixels are not independent
• Global Forest Change explained — what counts as loss

FAQ

What is avoided deforestation?

The forest that would have been lost without an intervention but was not: counterfactual loss minus actual loss.

Why is comparing inside and outside protected areas misleading?

Protected areas are often placed where clearing was unlikely, so the comparison mixes placement with protection.

What is matching?

Pairing protected and unprotected units with similar characteristics — road distance, slope, and so on — before comparing outcomes.

What is leakage?

Clearing displaced from protected land to nearby unprotected land.

Why does this matter for carbon credits?

Credits are based on the counterfactual. An inflated counterfactual issues credits for forest that was not threatened.

Are pixel-level statistics valid here?

Point estimates can be, but standard errors are not; nearby pixels are not independent.