How to Map Forest Loss Hotspots with H3 in Python

On this page

Problem statement

A 30 m forest-loss map of a region has tens of millions of pixels, and nobody can read it. Aggregating loss to hexagons of a few square kilometres turns it into a map of where clearing is concentrated — the view that monitoring agencies publish and enforcement teams act on. H3 hexagons are the natural grid: equal-ish area, a hierarchy of resolutions, and stable cell IDs that make maps from different years directly comparable.

"Hotspot" should mean more than "a dark hexagon". A local statistic such as Getis-Ord Gi* tests whether a cluster of high-loss cells is more concentrated than chance, which is the difference between a map and a finding.

Quick answer

Convert forest-loss pixels to points, index them to H3, sum area per cell, and normalise by the cell's forest baseline:

import numpy as np, pandas as pd, geopandas as gpd, h3, rasterio
from shapely.geometry import Polygon

with rasterio.open("rondonia_gfc.tif") as src:          # from the reading guide
    treecover, lossyear, datamask = src.read()
    tf = src.transform
rows, cols = np.nonzero((treecover > 30) & (datamask == 1))
xs, ys = rasterio.transform.xy(tf, rows, cols)
lat = np.array(ys)
px_ha = (abs(tf.a) * 111_320 * np.cos(np.radians(lat))) * (abs(tf.e) * 110_574) / 1e4

RES = 7                                                 # ~5 km² per cell
cells = [h3.latlng_to_cell(y, x, RES) for y, x in zip(ys, xs)]
df = pd.DataFrame({"cell": cells, "ha": px_ha, "year": lossyear[rows, cols]})
g = df.groupby("cell").agg(forest2000_ha=("ha", "sum"))
g["lost_ha"] = df[df.year > 0].groupby("cell").ha.sum()
g["lost_2020_24_ha"] = df[df.year >= 20].groupby("cell").ha.sum()
g = g.fillna(0)
g["share_lost"] = g.lost_ha / g.forest2000_ha
g["recent_share"] = g.lost_2020_24_ha / (g.forest2000_ha - g.lost_ha + g.lost_2020_24_ha)

hexes = gpd.GeoDataFrame(g, geometry=[Polygon([(lng, la) for la, lng in h3.cell_to_boundary(c)])
                                      for c in g.index], crs=4326)
print(f"{len(hexes)} H3 cells at resolution {RES} "
      f"(mean area {h3.average_hexagon_area(RES, 'km^2'):.1f} km²)")
print(hexes[["forest2000_ha", "lost_ha", "share_lost", "recent_share"]].describe(
    percentiles=[.5, .9]).loc[["50%", "90%", "max"]].round(2).to_string())
hexes.to_file("loss_hexes.gpkg")

652 H3 cells at resolution 7 (mean area 5.2 km²)
     forest2000_ha  lost_ha  share_lost  recent_share
50%         201.43    78.71        0.43          0.06
90%         383.79   193.39        0.66          0.22
max         509.07   507.34        1.00          1.00

Two normalisations answer two questions: the share of the 2000 forest ever lost shows where the frontier has been; the share of forest standing in 2020 that was lost since shows where it is now.

A resolution for the audience, index pixels, normalise

1. Choose the resolution for the audience

H3 resolution 7 (about 5 km²) suits regional overviews; 8 (about 0.7 km²) local ones. Too fine and most cells are empty; too coarse and hotspots dissolve into their surroundings.

2. Index pixels, not polygons

Assign every forest pixel's centre to a cell, with the pixel's area. That keeps areas exact and avoids polygon clipping.

3. Normalise by the baseline

Hectares lost per cell reward cells with more forest to lose. Divide by the forest in the cell at the start of the period for a rate.

4. Separate the past from the present

The cumulative share lost since 2000 maps the historical frontier; the share of remaining forest lost in the last few years maps the active one.

5. Test for clustering

Gi* compares each cell's neighbourhood sum with the expectation under randomness. Cells with high z-scores are statistically significant hotspots.

Pixels to cells, cells to rates, rates to hotspots.

Getis-Ord hotspots on H3 cells

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 — Getis-Ord Gi*

from libpysal import graph
from esda import G_Local

hexes = hexes.reset_index(names="cell")
w = graph.Graph.build_contiguity(hexes, rook=False)
gi = G_Local(hexes.recent_share.to_numpy(), w, star=True, permutations=999, seed=0)
hexes["gi_z"] = gi.Zs
hexes["hot"] = (gi.Zs > 1.96) & (gi.p_sim < 0.05)
print(f"significant hotspots of recent loss: {hexes.hot.sum()} of {len(hexes)} cells")
print(f"they hold {hexes.loc[hexes.hot, 'lost_2020_24_ha'].sum() / hexes.lost_2020_24_ha.sum():.0%} "
      f"of 2020–24 loss on {hexes.loc[hexes.hot].shape[0] / len(hexes):.0%} of the cells")

significant hotspots of recent loss: 10 of 652 cells
they hold 21% of 2020–24 loss on 2% of the cells

Ten cells, one and a half percent of the grid, hold a fifth of the loss since 2020. That concentration — not the colour of individual cells — is what makes a hotspot map useful for targeting.

Example 2 — resolution sensitivity

for res in (6, 7, 8):
    c = [h3.latlng_to_cell(y, x, res) for y, x in zip(ys, xs)]
    d = pd.DataFrame({"cell": c, "ha": px_ha, "recent": lossyear[rows, cols] >= 20})
    share = d.groupby("cell").apply(lambda x: (x.ha * x.recent).sum() / x.ha.sum(),
                                    include_groups=False)
    top = share.sort_values(ascending=False)
    k = max(1, int(len(top) * 0.1))
    print(f"res {res}: {len(top):>5} cells, top 10% of cells have mean recent share "
          f"{top.iloc[:k].mean():.1%}")

res 6:   108 cells, top 10% of cells have mean recent share 19.3%
res 7:   652 cells, top 10% of cells have mean recent share 26.2%
res 8:  4300 cells, top 10% of cells have mean recent share 36.0%

Finer cells find more extreme values because they average over less; coarser cells smooth them. A hotspot threshold set at one resolution does not transfer to another.

Example 3 — cells that stopped being frontier

early = df[(df.year > 0) & (df.year <= 10)].groupby("cell").ha.sum() / g.forest2000_ha
late = g.lost_2020_24_ha / g.forest2000_ha
both = pd.DataFrame({"early": early, "late": late}).fillna(0)
moved = both[(both.early > both.early.quantile(0.9)) & (both.late < both.late.median())]
print(f"cells in the top 10% for 2001–10 loss but below median for 2020–24: {len(moved)}")

cells in the top 10% for 2001–10 loss but below median for 2020–24: 47

Forty-seven cells were in the top tenth for early loss and are below the median now. A quiet cell has two readings: clearing stopped, or there is little forest left to clear and the frontier has moved on. Check each cell's remaining forest before reading its quiet as success.

A few hotspot cells, a large share of the recent loss.

Hexagons, normalisation and significance

Why hexagons

Hexagons have six equidistant neighbours, which makes neighbourhood statistics like Gi* cleaner than on squares, and H3's global index lets maps of different years and regions share cell IDs. Their areas vary slightly across the globe, which is why sums of pixel areas, not cell counts, carry the numbers.

Why normalise by the remaining forest

A cell cleared long ago has no forest left to lose; its low recent loss says nothing about pressure. Dividing recent loss by the forest still standing at the start of the recent period measures the rate at which what is left is disappearing.

Why test for significance

Random variation produces some high cells anywhere. Gi* asks whether a cell and its neighbours together are higher than expected, with permutation tests giving a p-value. It turns "these look high" into "these are unlikely by chance", with the usual caveats about multiple testing.

Finer cells find more extreme values because they average over less.

Window edges, multiple testing and units

Window edges: cells cut by the data window have partial baselines; drop cells whose forest area is well below the cell area without water.
Multiple testing: with hundreds of cells, some pass p < 0.05 by chance; use a false discovery rate correction for formal claims.
Pixel centres near cell edges can fall in either cell; at 30 m against 5 km² cells the effect is negligible.
H3 versions: the h3 Python API changed names in version 4 (latlng_to_cell, cell_to_boundary).
Area units: h3.average_hexagon_area is an average; use actual sums for statistics.

Internal links

How to aggregate points into H3 hexagons — H3 aggregation in general
How to calculate forest loss by year for a region — administrative units instead
Hot spot analysis in Python — Gi* in more depth
Deforestation rates explained — normalising loss
How to read Hansen Global Forest Change data — the pixel source
Spatial autocorrelation explained — the statistics behind hotspots

FAQ

How do I aggregate forest loss to hexagons?

Index each forest pixel's centre to an H3 cell and sum pixel areas by cell and loss year.

Which H3 resolution should I use?

Resolution 7 (about 5 km²) for regional overviews; 8 for local detail.

Should I map hectares lost or a rate?

A rate — loss divided by the forest present at the start of the period — unless the question is about absolute area.

What makes a cell a hotspot?

A significantly high local Gi* statistic, not just a high value.

Why did an old hotspot go quiet?

Often because its forest is gone; normalising by remaining forest shows whether pressure fell.

Does the resolution change the hotspots?

Yes. Finer cells produce more extreme values; compare maps only at the same resolution.