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

How to Measure Forest Fragmentation in Python

On this page

Problem statement

Deforestation totals say how much forest disappeared; they do not say what is left. Two landscapes that each lost 40% of their forest can look completely different — one with a single large block remaining, the other with thousands of small fragments separated by pasture. For biodiversity, carbon and fire risk the second is far worse, because forest edges dry out, burn more easily and lose the species that need interior conditions.

Fragmentation metrics describe that structure: how many patches, how big, how much edge, how much interior far enough from an edge to be unaffected. They are simple to compute from a forest mask and sensitive to choices — connectivity rule, pixel size, edge depth — that have to be reported with them.

Quick answer

Label connected forest patches, then compute patch statistics, edge length and core area:

import numpy as np, rasterio
from scipy import ndimage

with rasterio.open("rondonia_gfc.tif") as src:          # from the reading guide
    treecover, lossyear, datamask = src.read()
    tf = src.transform
px_m = abs(tf.a) * 111_320 * np.cos(np.radians(-10.55))  # ~27 m at this latitude
px_ha = px_m * abs(tf.e) * 110_574 / 1e4

f2000 = (treecover > 30) & (datamask == 1)
f2024 = f2000 & (lossyear == 0)

def fragmentation(forest, edge_depth_m=100):
    labels, n = ndimage.label(forest, structure=np.ones((3, 3)))       # 8-connectivity
    sizes = np.bincount(labels.ravel())[1:] * px_ha
    edge = forest & ~ndimage.binary_erosion(forest)                   # forest pixels on an edge
    core = ndimage.binary_erosion(forest, iterations=int(round(edge_depth_m / px_m)))
    total = forest.sum() * px_ha
    return {"forest ha": total, "patches": n, "mean patch ha": sizes.mean(),
            "largest patch %": 100 * sizes.max() / total,
            "edge m per ha": edge.sum() * px_m / (forest.size * px_ha),
            "core %": 100 * core.sum() * px_ha / total}

for label, f in (("2000", f2000), ("2024", f2024)):
    m = fragmentation(f)
    print(f"{label}: " + ", ".join(f"{k} {v:,.0f}" if k == "patches" else f"{k} {v:,.1f}"
                                   for k, v in m.items()))
2000: forest ha 139,688.6, patches 14,079, mean patch ha 9.9, largest patch % 57.3, edge m per ha 34.0, core % 49.9
2024: forest ha 78,383.1, patches 25,871, mean patch ha 3.0, largest patch % 4.3, edge m per ha 37.0, core % 28.3

The forest lost 44% of its area, but its structure changed faster: nearly twice as many patches, a largest patch that went from 57% of the forest to 4%, and a core — forest more than 100 m from any edge — that shrank by two-thirds.

Forest masks per date, a connectivity rule, patches and edge

1. Build the forest mask for each date

Forest in the baseline year from the canopy threshold; forest in a later year as baseline minus loss up to that year.

2. Choose the connectivity rule

Four-connectivity (edges only) splits diagonal connections that eight-connectivity (edges and corners) joins. Patch counts differ substantially; state the rule.

3. Label patches and summarise

Number of patches, mean and largest patch size, and the largest patch as a share of all forest are the core structure measures.

4. Compute edge

Count forest pixels adjacent to non-forest, or measure boundary length. Express per hectare of landscape to compare areas of different sizes.

5. Compute core area with an edge depth

Erode the forest mask by the depth of edge effects — often 100 m for microclimate, more for some species. Report core as a share of forest.

One forest block and the same forest cut by roads into fragments, with the core area more than 100 m from an edge shaded. the core goes firstintact blockcorecut by two roadslittle forest lost, most core goneRondônia: forest −44% by 2024, core area −68% spatialworkflow.io
The same forest area; the core is what fragmentation takes first.

Fragmentation year by year

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 — year by year

print(f"{'year':<6}{'forest %':>9}{'patches':>9}{'largest %':>11}{'core %':>8}")
for y in (2000, 2005, 2010, 2015, 2020, 2024):
    f = f2000 & ((lossyear == 0) | (lossyear > y - 2000)) if y > 2000 else f2000
    m = fragmentation(f)
    print(f"{y:<6}{100 * m['forest ha'] / (f2000.size * px_ha):>9.1f}{m['patches']:>9,}"
          f"{m['largest patch %']:>11.1f}{m['core %']:>8.1f}")
year   forest %  patches  largest %  core %
2000       46.2   14,079       57.3    49.9
2005       36.8   18,432        9.4    39.5
2010       33.3   20,185        7.6    36.6
2015       31.8   20,933        7.2    35.1
2020       28.8   23,825        5.0    32.7
2024       25.9   25,871        4.3    28.3

Core area falls faster than forest area, and the largest patch collapses early — from 57% to 9% of the forest within the first five years, while the forest share fell by only nine points. The timing of that collapse, not only the rate of loss, is what a fragmentation analysis adds.

Example 2 — four- or eight-connectivity

for name, st in (("4-connected", ndimage.generate_binary_structure(2, 1)),
                 ("8-connected", np.ones((3, 3)))):
    _, n = ndimage.label(f2024, structure=st)
    print(f"{name}: {n:,} patches in 2024")
4-connected: 40,547 patches in 2024
8-connected: 25,871 patches in 2024

The same forest, very different patch counts. Diagonal single-pixel links are common in a 30 m classification; whether they connect two patches is a modelling decision.

Example 3 — the edge depth

for depth in (30, 100, 300, 500):
    core = ndimage.binary_erosion(f2024, iterations=max(1, int(round(depth / px_m))))
    print(f"core at {depth:>3} m from an edge: {100 * core.sum() / f2024.sum():5.1f}% of the forest")
core at  30 m from an edge:  60.6% of the forest
core at 100 m from an edge:  28.3% of the forest
core at 300 m from an edge:   9.3% of the forest
core at 500 m from an edge:   3.7% of the forest

Core area depends steeply on the edge depth. Microclimate effects reach about 100 m into tropical forest; fire and some species effects reach much further, and at those depths a fragmented landscape has almost no core left.

Forest cover, patches, largest patch share and core share in the Rondônia window from 2000 to 2024. Rondônia window, 8-connected, 100 m edgeyearforestpatcheslargestcore200046.2%14,07957.3%49.9%200536.8%18,4329.4%39.5%201033.3%20,1857.6%36.6%201531.8%20,9337.2%35.1%202028.8%23,8255.0%32.7%202425.9%25,8714.3%28.3%The largest patch collapsed in the first five years spatialworkflow.io
Area, the largest patch and the core decline at different speeds.

Structure changes faster than area

Why structure changes faster than area

Clearing follows roads. Each new road and its side branches cut through large blocks long before most of the forest is gone, so the landscape passes from one connected forest to many fragments early. Area statistics lag the change in structure.

Why core area is the sensitive measure

Every hectare lost at an edge exposes new edge behind it. Core area therefore shrinks by the lost area plus a band of newly exposed forest, which makes it fall disproportionately fast — and it is the part of the forest that retains interior conditions.

Why the choices must be reported

Connectivity, edge depth, minimum patch size and pixel size all change the numbers. Fragmentation statistics are comparable only when those choices match, which is why landscape-ecology software such as FRAGSTATS asks for them explicitly.

Core share of the 2024 forest at four edge depths, in the Rondônia window. core share of the 2024 forest30 m from an edge60.6%100 m28.3%300 m9.3%500 m3.7%035%70%Four-connected patches: 40,547; eight-connected: 25,871 — state both parameters spatialworkflow.io
The edge depth decides whether a fragmented forest still has an interior.

Window edges, pylandstats and resolution

  • Window edges cut patches artificially; exclude patches touching the window border from size statistics, or read a larger area.
  • pylandstats implements FRAGSTATS-style metrics in Python for a fuller set.
  • Resolution: aggregating to coarser pixels merges patches and removes small gaps; see the fix guide.
  • Linear features narrower than a pixel (small roads, rivers) may not break patches in the map even where they do on the ground.
  • Regrowth: GFC does not map forest recovery after 2012; later forest masks treat all lost pixels as non-forest.

FAQ

How do I measure forest fragmentation in Python?

Label connected patches of a forest mask with scipy.ndimage.label, then compute patch counts and sizes, edge length and core area.

What is core area?

Forest further than a chosen depth — often 100 m — from any edge.

Should I use 4- or 8-connectivity?

Either, stated. Eight-connectivity joins diagonal pixels and gives fewer patches.

What is the largest patch index?

The largest patch as a share of the total forest (or landscape) area.

Why does core area fall faster than forest area?

Each clearing exposes a new band of edge in the remaining forest.

Is there a Python version of FRAGSTATS?

pylandstats implements many FRAGSTATS metrics.