⚙️ How-To Forest monitoring Biodiversity & species data Forest Monitoring · step 9 Notebook By · Updated

How to Compare Forest Loss Inside and Outside Protected Areas in Python

On this page

Problem statement

Do protected areas and indigenous territories keep their forest? The first comparison anyone makes — the share of forest lost inside against outside — almost always shows much less loss inside. Part of that is protection working. Part of it is that protected areas are placed where clearing was less likely anyway: far from roads, on steep or poor land, in remote places. A fair comparison looks at land that was similarly exposed, typically a band just outside the boundary, and separates the boundary effect from the location effect.

Quick answer

Rasterise the protected polygons, build distance bands inside and outside the boundary, and compare the share of 2000 forest lost in each:

import numpy as np, pandas as pd, geopandas as gpd, 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               # the northern edge of Uru-Eu-Wau-Wau
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)
treecover, tf = layer("treecover2000")
lossyear, _ = layer("lossyear")
datamask, _ = layer("datamask")
shape = treecover.shape
px_m = abs(tf.a) * 111_320 * np.cos(np.radians((S + N) / 2))
px_ha = px_m * abs(tf.e) * 110_574 / 1e4

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"])]
print("protected polygons:", pa[["name", "boundary"]].dropna(subset=["name"]).drop_duplicates("name")
      .head(6).to_dict("records"))

inside = features.rasterize(((g, 1) for g in pa.geometry), out_shape=shape, transform=tf) == 1
samp = (abs(tf.e) * 110_574, px_m)
d_out = ndimage.distance_transform_edt(~inside, sampling=samp) / 1000     # km outside
d_in = ndimage.distance_transform_edt(inside, sampling=samp) / 1000       # km inside
signed = np.where(inside, -d_in, d_out)                                    # negative = inside

forest = (treecover > 30) & (datamask == 1)
df = pd.DataFrame({"d": signed[forest], "lost": lossyear[forest] > 0})
df["band"] = pd.cut(df.d, [-100, -10, -5, -1, 0, 1, 5, 10, 100])
out = df.groupby("band", observed=True).agg(forest_ha=("lost", "size"), share_lost=("lost", "mean"))
out["forest_ha"] = (out.forest_ha * px_ha).round(-2)
print(out.round({"share_lost": 3}).to_string())
protected polygons: [{'name': 'Parque Nacional de Pacaás Novos', 'boundary': 'protected_area'}, {'name': 'Terra Indígena Uru-Eu-Wau-Wau', 'boundary': 'aboriginal_lands'}]
             forest_ha  share_lost
band                              
(-100, -10]   330700.0       0.012
(-10, -5]      64900.0       0.001
(-5, -1]       47500.0       0.016
(-1, 0]        11300.0       0.077
(0, 1]          7000.0       0.579
(1, 5]         24500.0       0.655
(5, 10]        22000.0       0.625
(10, 100]      24300.0       0.537

Negative distances are inside the protected land, positive outside. The step at zero is the most direct evidence of protection: 7.7% of the forest lost in the kilometre just inside the boundary, 58% in the kilometre just outside. Deeper inside, loss falls to 1–2%; further outside it stays above half. The first kilometre inside is where the boundary is being nibbled.

Protected boundaries, one baseline, signed distance

1. Get the protected boundaries

National registries or the World Database on Protected Areas are authoritative; OSM maps many protected areas (boundary=protected_area) and indigenous territories (boundary=aboriginal_lands) with names and designations.

2. Use the same forest baseline inside and out

Forest in 2000 by a canopy threshold, water masked. Protected areas often contain non-forest vegetation that should not enter the denominator.

3. Build signed distance to the boundary

Distance outside and distance inside, combined into one signed layer. Bands of a few kilometres on each side compare land with similar access.

4. Compare shares, not areas

Share of the 2000 forest lost per band. Areas depend on how much forest each band contains.

5. Report the placement caveat

A lower loss rate inside does not by itself measure the effect of protection. Matching on distance to roads, slope and other covariates is the standard next step.

Signed distance bands inside and outside a protected area boundary, with clearing concentrated outside. signed distance from the boundaryoutside ← boundary → inside+1 km: 58% lost−1 km: 7.7% lostThe step at zero is the most direct evidence; deeper bands add location effects spatialworkflow.io
Signed distance: the boundary effect is the step at zero.

Inside against the outer band, over time

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, since OpenStreetMap is edited daily and a live query would move the counts.

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-forest-loss-protected-areas-python-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 — loss over time, inside against the outer band

near_out = forest & (signed > 0) & (signed <= 5)
near_in = forest & (signed < 0) & (signed >= -5)
for y0, y1 in ((1, 5), (6, 10), (11, 15), (16, 20), (21, 24)):
    lo = lambda m: (m & (lossyear >= y0) & (lossyear <= y1)).sum() / m.sum()
    print(f"{2000 + y0}–{2000 + y1}: lost {lo(near_out):6.2%} of forest 0–5 km outside, "
          f"{lo(near_in):6.2%} 0–5 km inside")
2001–2005: lost 26.04% of forest 0–5 km outside,  0.86% 0–5 km inside
2006–2010: lost 10.69% of forest 0–5 km outside,  0.41% 0–5 km inside
2011–2015: lost  6.02% of forest 0–5 km outside,  0.23% 0–5 km inside
2016–2020: lost  9.88% of forest 0–5 km outside,  0.50% 0–5 km inside
2021–2024: lost 11.19% of forest 0–5 km outside,  0.73% 0–5 km inside

The contrast between the two sides holds in every period: clearing just outside runs 15 to 30 times the rate just inside, though the ratio has narrowed since 2011. Both sides slowed to 2015 and rose again after; inside, the share cleared per period has tripled from 0.23% to 0.73%. A rising inside share is the pattern monitoring programmes look for.

Example 2 — are the two sides comparable?

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)
dist_road = ndimage.distance_transform_edt(burn == 0, sampling=samp) / 1000
for label, m in (("0–5 km outside", near_out), ("0–5 km inside", near_in)):
    print(f"{label:<15} median distance to a mapped road {np.median(dist_road[m]):5.1f} km")
0–5 km outside  median distance to a mapped road   1.0 km
0–5 km inside   median distance to a mapped road   3.5 km

Even within 5 km of the boundary, the two sides differ in access. Part of the difference in loss is roads, not protection — which is why careful studies match pixels on distance to roads before comparing.

Example 3 — a matched comparison on road distance

bins = [0, 1, 2, 5, 10, 20, 100]
rows = []
for lo, hi in zip(bins[:-1], bins[1:]):
    band = (dist_road >= lo) & (dist_road < hi)
    a, b = forest & band & near_out, forest & band & near_in
    if a.sum() > 1000 and b.sum() > 1000:
        rows.append((f"{lo}–{hi} km", (lossyear[a] > 0).mean(), (lossyear[b] > 0).mean()))
print(f"{'road distance':<14}{'outside':>9}{'inside':>9}")
for r in rows:
    print(f"{r[0]:<14}{r[1]:>9.1%}{r[2]:>9.1%}")
road distance   outside   inside
0–1 km            64.2%    12.7%
1–2 km            66.2%     6.4%
2–5 km            57.9%     1.6%

Within the same road-distance class, the boundary still makes a large difference — 64% against 13% within a kilometre of a road — in the three classes where both sides have enough forest to compare. Beyond 5 km from a road one side or the other has too little forest near the boundary to compare. That residual gap is closer to the effect of protection than the raw comparison.

Share of the 2000 forest lost by 2024 in distance bands inside and outside the protected boundary. share of 2000 forest lost, by signed distanceinside, > 10 km1.2%inside, 5–10 km0.1%inside, 1–5 km1.6%inside, 0–1 km7.7%outside, 0–1 km57.9%outside, 1–5 km65.5%outside, 5–10 km62.5%outside, > 10 km53.7%035%70%Pacaás Novos National Park and the Uru-Eu-Wau-Wau territory spatialworkflow.io
Loss share by signed distance from the boundary.

Why raw comparisons overstate protection

Why the raw inside/outside comparison overstates protection

Protected areas are rarely placed at random. Governments protect land that is remote, steep or otherwise unattractive to clear, partly because that is where protection is cheapest. Comparing them with all unprotected land compares remote with accessible land.

Why the boundary band helps

Land just outside a boundary is similar to land just inside in climate, soil and remoteness. A sharp change in clearing at the boundary is hard to explain by anything other than the boundary itself. It is the spatial equivalent of a regression discontinuity.

Why matching is the standard

Covariate matching pairs protected pixels with unprotected pixels of similar distance to roads, slope, distance to past clearing and so on. The difference in loss between matched pairs estimates avoided deforestation with far less placement bias than a simple comparison.

Share of forest lost per period within 5 km outside and inside the protected boundary. both sides of the boundary, by periodperiod0–5 km outside0–5 km inside2001–0526.0%0.86%2006–1010.7%0.41%2011–156.0%0.23%2016–209.9%0.50%2021–2411.2%0.73%Clearing inside has tripled since 2011–15: the pattern monitoring looks for spatialworkflow.io
Outside runs 15 to 30 times the inside rate; the inside rate is rising again.

Leakage, designation dates and OSM completeness

  • Leakage: protection can push clearing into neighbouring land, inflating loss just outside the boundary.
  • Designation dates: areas protected after 2001 should be compared from their designation year.
  • Different designations (strict reserves, sustainable-use areas, indigenous lands) have different rules and effects.
  • OSM completeness for protected areas varies; use the national registry for official analyses.
  • Loss inside protected areas includes fire and natural disturbance, not only illegal clearing.

FAQ

How do I measure deforestation inside protected areas?

Rasterise the protected polygons onto the forest grid and compute the share of baseline forest lost inside and in comparison bands outside.

Why compare with a band just outside the boundary?

Land on both sides of the boundary is similar in access and environment, so a difference there is more likely due to protection.

Does lower loss inside prove protection works?

No — protected areas are often placed where clearing was less likely. Match on covariates to estimate the effect.

What is leakage?

Clearing displaced from protected land into nearby unprotected land.

Where do I get protected area boundaries?

National registries, the World Database on Protected Areas, or OSM tags boundary=protected_area and boundary=aboriginal_lands.

How do I handle areas protected after 2001?

Compare loss only from the year of designation onward.