⚙️ How-To Forest monitoring Raster & rasterio Forest Monitoring · step 8 Notebook By · Updated

How to Calculate Forest Loss by Year for a Region in Python

On this page

Problem statement

A forest-loss figure is almost always wanted per something: per municipality, per concession, per watershed, per protected area, per year. The Global Forest Change layers give a year of loss for every 30 m pixel; turning that into a table of hectares by polygon by year is a zonal statistics problem with two details that decide whether the numbers are right. Pixel areas vary with latitude on the degree grid, and the baseline — forest in 2000 above some canopy threshold — has to be applied before counting loss, or loss in scrub and pasture is counted as deforestation.

Quick answer

Rasterise the polygons onto the GFC grid, then sum pixel areas by polygon and loss year:

import numpy as np, pandas as pd, geopandas as gpd, osmnx as ox, rasterio
from rasterio import features

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

lat = N - (np.arange(shape[0]) + 0.5) * abs(tf.e)
px_ha = (abs(tf.a) * 111_320 * np.cos(np.radians(lat))) * (abs(tf.e) * 110_574) / 1e4
area = np.broadcast_to(px_ha[:, None], shape)

ox.settings.http_user_agent = "spatialworkflow-examples/1.0"
mun = ox.features_from_bbox((W, S, E, N), {"boundary": "administrative", "admin_level": "8"})
mun = mun[(mun.get("admin_level") == "8") & mun.geom_type.isin(["Polygon", "MultiPolygon"])]
mun = mun[["name", "geometry"]].drop_duplicates("name").reset_index(drop=True)   # OSM can return a boundary twice
zone = features.rasterize(((g, i + 1) for i, g in enumerate(mun.geometry)), out_shape=shape,
                          transform=tf, fill=0)

forest = (treecover > 30) & (datamask == 1)
df = pd.DataFrame({"zone": zone[forest], "year": lossyear[forest].astype(int), "ha": area[forest]})
table = df[df.year > 0].pivot_table(index="zone", columns="year", values="ha", aggfunc="sum",
                                    fill_value=0)
table.index = mun.name.reindex(table.index - 1).values
table.columns = table.columns + 2000
base = df.groupby("zone").ha.sum()
summary = pd.DataFrame({"forest 2000 (ha)": base.values}, index=mun.name.reindex(base.index - 1).values)
summary = summary.drop(index=[np.nan], errors="ignore")
summary["lost 2001–24 (ha)"] = table.sum(axis=1)
summary["share lost"] = summary["lost 2001–24 (ha)"] / summary["forest 2000 (ha)"]
print(summary.round({"forest 2000 (ha)": -2, "lost 2001–24 (ha)": -2, "share lost": 2})
      .sort_values("share lost", ascending=False).head(8).to_string())
                           forest 2000 (ha)  lost 2001–24 (ha)  share lost
Governador Jorge Teixeira           58900.0            30200.0        0.51
Ouro Preto do Oeste                   400.0              200.0        0.51
Jaru                                41100.0            16400.0        0.40
Cacaulândia                         33900.0            12900.0        0.38
Nova União                           3300.0             1200.0        0.36
Theobroma                            2000.0              400.0        0.21

Each municipality's share of its 2000 forest lost by 2024. Because the window only covers part of most municipalities, these are figures for the part inside the window — the next step for a real report is reading enough tiles to cover whole polygons.

Polygons on the raster grid, baseline first, area weights

1. Put the polygons on the raster grid

Rasterising the polygons onto the GFC transform gives an integer zone for every pixel, which is fast and exact to the pixel. For small polygons, all_touched and area-weighted fractions matter.

2. Apply the forest baseline first

Count only pixels that were forest in 2000 by your threshold. Loss outside that mask is disturbance of non-forest cover.

3. Weight every pixel by its area

Sum hectares, not pixels. On the degree grid each row has its own pixel area.

4. Pivot by zone and year

A zone × year table is the standard output. Keep the baseline beside it so shares can be computed.

5. Check coverage of each polygon

If a polygon extends beyond the data you read, its totals are partial. Compare the polygon's area with the area of its pixels.

From the forest change layers and polygons through rasterising, the baseline and area weighting to a zone by year table. zonal loss, the order mattersrasterise polygonszone id per pixelforest 2000 maskthreshold + landpixel areasper row, by latitudepivotzone × year, hectaresBaseline first, area-weighted sums second, coverage check last spatialworkflow.io
Baseline first, area-weighted sums second.

The annual series for one municipality

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-forest-loss-by-year-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 — coverage check

covered = pd.Series(np.bincount(zone.ravel(), weights=area.ravel()))
for i, row in mun.iterrows():
    poly_ha = gpd.GeoSeries([row.geometry], crs=4326).to_crs("ESRI:102033").area.iloc[0] / 1e4      # South America Albers, equal-area
    got = covered.get(i + 1, 0)
    if poly_ha > 0:
        print(f"{row['name']:<24} polygon {poly_ha:>10,.0f} ha, inside the window {got / poly_ha:5.0%}")
Nova União               polygon     80,682 ha, inside the window   10%
Cacaulândia              polygon    196,434 ha, inside the window   32%
Ouro Preto do Oeste      polygon    197,324 ha, inside the window    0%
Jaru                     polygon    294,542 ha, inside the window   39%
Theobroma                polygon    219,705 ha, inside the window    3%
Governador Jorge Teixeira polygon    506,438 ha, inside the window   22%

Municipalities cut by the window edge are reported for their covered part only. The check prevents presenting a partial total as a municipal one.

Example 2 — the annual series for one municipality

top = summary["lost 2001–24 (ha)"].idxmax()
s = table.loc[top]
print(f"{top}: loss per year (ha)")
print(" ".join(f"{y}:{v:,.0f}" for y, v in s.items() if y % 3 == 1 or y == 2024))
print(f"five-year totals: 2001–05 {s.loc[2001:2005].sum():,.0f}, "
      f"2011–15 {s.loc[2011:2015].sum():,.0f}, 2020–24 {s.loc[2020:2024].sum():,.0f}")
Governador Jorge Teixeira: loss per year (ha)
2002:3,723 2005:1,429 2008:796 2011:886 2014:464 2017:836 2020:767 2023:439 2024:992
five-year totals: 2001–05 12,915, 2011–15 2,251, 2020–24 5,216

Annual numbers are noisy; multi-year totals show the trajectory more reliably. Comparisons across the whole 2001–2024 record carry the caveat that detection methods changed.

Example 3 — threshold sensitivity per municipality

for t in (10, 30, 50):
    f = (treecover > t) & (datamask == 1)
    lost = (f & (lossyear > 0))
    print(f"threshold {t:>2}%: forest {area[f].sum():>9,.0f} ha, lost {area[lost].sum():>9,.0f} ha, "
          f"share {area[lost].sum() / area[f].sum():.1%}")
threshold 10%: forest   151,295 ha, lost    62,185 ha, share 41.1%
threshold 30%: forest   139,684 ha, lost    61,301 ha, share 43.9%
threshold 50%: forest   129,746 ha, lost    59,897 ha, share 46.2%

The threshold changes both the baseline and the loss. Shares move less than absolute areas, but they move.

Forest loss in the Rondônia window in selected years from 2001 to 2024. loss on 2000 forest, km², selected years200177.62002 (peak)97200436.9200814.02012 (lowest)4.6201619.8202117.4202416.5050100The record spans several algorithm versions; compare periods, not single years spatialworkflow.io
Annual loss in the window; the record spans several algorithm versions.

Rasterise once, baseline first

Why rasterise the polygons rather than mask per polygon

Rasterising all polygons once gives a zone array; every statistic is then a grouped sum, which is fast for hundreds of polygons. Masking the raster separately for each polygon repeats work and is slower for large sets.

Why the baseline has to come first

GFC loss is any stand-replacing disturbance of tree cover, including in areas with little tree cover to begin with. Without a baseline threshold, clearing of a 15% canopy pasture counts the same as clearing closed forest, and "deforestation" totals inflate.

Why partial coverage matters

Official statistics are reported per whole administrative unit. A window that clips a municipality produces a correct figure for an arbitrary part of it — useful for the window, misleading if labelled with the municipality's name.

Forest in 2000 and loss by 2024 for the parts of Rondônia municipalities inside the window. by municipality, inside the window onlymunicipalityforest 2000lost 2001–24sharepolygon insideGov. Jorge Teixeira58,900 ha30,200 ha51%22%Jaru41,100 ha16,400 ha40%39%Cacaulândia33,900 ha12,900 ha38%32%Nova União3,300 ha1,200 ha36%10%Partial polygons: the next step is reading enough tiles to cover them spatialworkflow.io
A municipal table is only municipal if the whole polygon was read.

Equal-area CRS, small and overlapping polygons

  • Equal-area CRS: for national-scale work, reproject once to an equal-area CRS instead of using row areas.
  • Small polygons: pixel-centre rasterisation mis-assigns edges; compute coverage fractions for polygons of a few hundred hectares or less.
  • Overlapping polygons (e.g. protected areas overlapping municipalities) need separate zone arrays.
  • Loss attribution: GFC does not say why trees were lost; drivers need a separate dataset.
  • Version: store the GFC version with the table.

FAQ

How do I calculate forest loss per municipality?

Rasterise the municipalities onto the GFC grid, mask to forest in 2000, and sum area-weighted pixels by zone and loss year.

Why do I need a forest baseline?

GFC loss includes disturbance of sparse tree cover; the baseline limits it to what you define as forest.

How do I get pixel areas right?

Multiply by each row's pixel area on the degree grid, or reproject to an equal-area CRS.

Can I sum annual losses to get total loss?

Yes, within one version of the dataset; comparisons of early and late years carry method caveats.

What if my polygon crosses tiles?

Read and mosaic all intersecting tiles first; report coverage per polygon.

Does GFC say why the forest was lost?

No. Drivers come from separate classifications.