How to Measure Deforestation Near Roads in Python

On this page

Problem statement

In the Amazon, clearing follows roads. Official highways open a frontier; settlement roads branch off them at regular intervals, giving Rondônia its "fishbone" pattern; unofficial tracks push further in. Quantifying that relationship — how much of the forest within 1 km of a road was lost, against 10 km away — is one of the most common analyses of deforestation drivers, and a basic input to models that project where clearing will happen next.

The analysis is a distance transform and a grouped sum. What makes it credible is the road layer: official road maps miss most of the tracks that clearing actually uses, and a map of distance to the roads that happen to be mapped measures mapping as much as deforestation.

Quick answer

Rasterise the roads onto the forest grid, compute distance, and summarise loss by distance band:

import numpy as np, pandas as pd, osmnx as ox, rasterio
from rasterio import features
from scipy import ndimage

with rasterio.open("rondonia_gfc.tif") as src:          # from the reading guide
    treecover, lossyear, datamask = src.read()
    tf, shape = src.transform, (src.height, src.width)
    W, S, E, N = src.bounds
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"
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"])]
print("road types:", roads.highway.value_counts().head(6).to_dict())

burn = features.rasterize(((g, 1) for g in roads.geometry), out_shape=shape, transform=tf,
                          all_touched=True)
dist_km = ndimage.distance_transform_edt(burn == 0, sampling=(abs(tf.e) * 110_574, px_m)) / 1000

forest = (treecover > 30) & (datamask == 1)
df = pd.DataFrame({"km": dist_km[forest], "lost": lossyear[forest] > 0})
df["band"] = pd.cut(df.km, [0, 0.5, 1, 2, 5, 10, 50], include_lowest=True)
out = df.groupby("band", observed=True).agg(forest_ha=("lost", "size"), share_lost=("lost", "mean"))
out["forest_ha"] = out.forest_ha * px_ha
print(out.round({"forest_ha": -2, "share_lost": 2}).to_string())

road types: {'residential': 965, 'unclassified': 801, 'service': 313, 'track': 195, 'secondary': 106, 'tertiary': 93}
               forest_ha  share_lost
band                                
(-0.001, 0.5]    37800.0        0.42
(0.5, 1.0]       41000.0        0.46
(1.0, 2.0]       54400.0        0.45
(2.0, 5.0]        6400.0        0.35
(5.0, 10.0]        100.0        0.24

The gradient is flatter than the fishbone pattern suggests: 42–46% of the 2000 forest lost within 2 km of a mapped way, 35% at 2–5 km, 24% in the thin band beyond. The reason is saturation. Almost all of the forest lies within 2 km of some mapped way, so nearly every pixel is "near a road" and distance has little left to discriminate. Distance to roads explains clearing where distance varies — across a frontier and the intact forest beyond it — not inside a settled frontier.

A complete road layer, rasterised, metric distances

1. Assemble the most complete road layer you can

OSM highway=* in the Amazon includes official roads and many track and unclassified ways traced from imagery. Research road datasets add tracks mapped specifically for deforestation studies. Official road maps alone miss most of the network.

2. Rasterise onto the analysis grid

Burn every road into the same grid as the forest data, with all_touched=True so narrow roads are not lost between pixel centres.

3. Compute distance in metres

distance_transform_edt with sampling set to the pixel size in metres (different in x and y on a degree grid) gives true distances.

4. Summarise loss by distance band

Share of the 2000 forest lost in each band, and forest area per band so small bands are not over-interpreted.

5. Beware of the reverse arrow

Many tracks were built to clear the forest, after the fact. Distance to today's roads explains past clearing partly because clearing created the roads. For causal questions, use the road network as of the start of the period.

The fishbone: clearing spreads from a highway along its side roads.

Which roads matter, and distance at the time of 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. Roads come from OpenStreetMap through osmnx, via the maps.mail.ru Overpass mirror while overpass-api.de was refusing connections. The OpenStreetMap answers are read as the Overpass mirror gave them on 10 October 2026.

import io
import os
import requests
import zipfile

import rasterio
from rasterio.windows import from_bounds

# OpenStreetMap as Overpass returned it on 10 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-deforestation-distance-to-roads-python-2026-10-10.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

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)

import osmnx as ox

ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"    # overpass-api.de refused us

3 (2000, 2000) EPSG:4326 (0.00025, 0.00025)

Example 1 — which roads matter

main = roads[roads.highway.isin(["trunk", "primary", "secondary", "tertiary"])]
burn_main = features.rasterize(((g, 1) for g in main.geometry), out_shape=shape, transform=tf,
                               all_touched=True)
d_main = ndimage.distance_transform_edt(burn_main == 0, sampling=(abs(tf.e) * 110_574, px_m)) / 1000
lost = forest & (lossyear > 0)
for label, d in (("all mapped ways", dist_km), ("main roads only", d_main)):
    print(f"{label:<17} median distance of lost forest {np.median(d[lost]):5.1f} km, "
          f"of remaining forest {np.median(d[forest & ~lost]):5.1f} km")

all mapped ways   median distance of lost forest   0.9 km, of remaining forest   0.9 km
main roads only   median distance of lost forest   3.0 km, of remaining forest   3.3 km

With all mapped ways, lost and remaining forest sit at the same median distance, 0.9 km: the side roads reach everything, cleared or not. With main roads only, both are around 3 km away, and lost forest is slightly closer (3.0 against 3.3 km). Neither layer separates cleared from standing forest here; which layer is used changes the distances by a factor of three.

Example 2 — distance at the time of loss

# Loss in 2001–2005 against distance to all roads mapped today
early = forest & (lossyear >= 1) & (lossyear <= 5)
late = forest & (lossyear >= 20)
print(f"median distance, loss 2001–05: {np.median(dist_km[early]):.2f} km")
print(f"median distance, loss 2020–24: {np.median(dist_km[late]):.2f} km")
print(f"remaining forest in 2024: {np.median(dist_km[forest & (lossyear == 0)]):.2f} km")

median distance, loss 2001–05: 0.90 km
median distance, loss 2020–24: 0.88 km
remaining forest in 2024: 0.89 km

Early loss, recent loss and the forest still standing all sit a median 0.9 km from today's roads. By 2024 the network reaches the whole window, so today's roads cannot tell which forest was cleared first. The road layer is a snapshot of the end of the period; a causal question needs the network as it was at the start, which a current OSM extract cannot supply.

Example 3 — save the distance layer

profile = {"driver": "GTiff", "height": shape[0], "width": shape[1], "count": 1,
           "dtype": "float32", "crs": "EPSG:4326", "transform": tf, "compress": "deflate"}
with rasterio.open("dist_to_road_km.tif", "w", **profile) as dst:
    dst.write(dist_km.astype("float32"), 1)
print(f"distance to nearest mapped road: max {dist_km.max():.1f} km, "
      f"share of window within 2 km {np.mean(dist_km <= 2):.0%}")

distance to nearest mapped road: max 5.5 km, share of window within 2 km 97%

The distance layer is a standard covariate for deforestation risk models, alongside slope, soils, protected status and distance to past clearing.

The share of the 2000 forest lost, by distance to the nearest mapped way.

Why clearing follows roads

Roads lower the cost of reaching land, removing timber and bringing cattle and crops to market. In Rondônia, the colonisation schemes of the 1970s and 80s laid out settlement roads at regular intervals perpendicular to the BR-364 highway, and lots along them were cleared progressively — the fishbone visible from space.

Why road completeness drives the result

If the tracks used for clearing are not in the road layer, clearings appear to be far from roads and the relationship weakens. The more complete the road layer, the stronger and more local the relationship looks. Two studies with different road data can reach different conclusions about the same forest.

Why causality is hard

Roads cause clearing and clearing causes roads. Using a road network from before the period, or modelling road building explicitly, is how studies separate the two. A distance-band table from today's roads describes association, not cause.

Lost and standing forest sit at the same distance from today’s network: the saturation problem.

Rivers, padded queries and mapped tracks

Rivers are transport routes in much of the Amazon; add them as a second distance layer.
Padding the road query beyond the window avoids edge cells being far from roads that exist just outside.
Degree grids need anisotropic sampling in the distance transform; equal-area projection avoids it.
Mapped tracks in OSM are unevenly complete; check a sample against imagery.
Protected areas often have fewer roads by design, confounding road and protection effects.

Internal links

How to read Hansen Global Forest Change data — the forest stack
How to measure forest fragmentation — what roads do to structure
How to compare forest loss inside and outside protected areas — another driver
How to measure distance to the grid — the same distance transform in energy siting
OpenStreetMap data model explained — highway tags
Deforestation rates explained — rates rather than shares

FAQ

How do I calculate distance to roads for a raster?

Rasterise the roads onto the grid and run scipy.ndimage.distance_transform_edt on the inverted mask, with sampling in metres.

Which roads should I include?

All mapped ways that vehicles can use, including tracks; main roads alone miss where clearing happens.

Why is deforestation concentrated near roads?

Access lowers the cost of clearing and of getting products to market.

Do roads cause deforestation?

They enable it, and clearing also creates roads. Use a road network from before the period for causal questions.

Is OSM complete enough in the Amazon?

Partly. Many tracks are mapped from imagery; check completeness against imagery for your area.

How do I handle the degree grid?

Pass pixel sizes in metres to sampling, or reproject to an equal-area CRS first.