Forestry GIS Python · Forestry · Ecology

Generating a Hexagonal Plot Grid with GeoPandas

A systematic sampling lattice has to answer one question well: where do I put n plots so that no part of the ownership is over- or under-represented? A square grid is the obvious answer and a slightly worse one, because a square’s corners are much further from its centre than its edges are, so coverage is uneven in a way that shows up in variance estimates. A hexagonal lattice fixes that with the same amount of code. This guide builds one, as a concrete design inside Spatial Plot Sampling Design, part of Ecological GIS Data Foundations in Python.

When to use this approach

Layout Coverage evenness Bias risk Use when
Hexagonal lattice Best — every point is closest to a near-equidistant centre Aliasing against periodic features Systematic inventory, monitoring networks, habitat surveys
Square grid Good, but corner-biased Aliasing against roads, rows, and drainage Compatibility with an existing grid or a national frame
Simple random Even in expectation, clumpy in practice None, but high variance Small samples where design-based inference is essential
Stratified random Even within strata Depends on the strata Strata that genuinely differ, as in Stratified Random Sampling for Forest Plots
Hexagonal with random start Best, plus randomisation Low The practical default for a repeatable inventory

The hexagonal lattice with a random origin is the design most national inventories converge on, because it combines even coverage with a randomisation that supports design-based estimation.

Spacing from a target plot count

The geometry is fixed. A regular hexagon of circumradius — centre to vertex — has area

and hexagons tile without gaps, so for a target of plots over an area the spacing follows directly:

For 200 plots over a 12,000 ha ownership that gives m and a centre-to-centre spacing of about 263 m. Working from the target count rather than picking a round spacing is what keeps the sample size predictable when the ownership boundary is irregular.

Coverage evenness of a square lattice against a hexagonal lattice Two lattices of equal cell area are drawn side by side with a sample point at every cell centre and a circle of equal radius around each point. In the square lattice the circles leave conspicuous four-cornered gaps where the corners of each cell are much further from its centre than its edges. In the hexagonal lattice the same circles leave smaller, more uniform gaps because a hexagon is closer to a circle than a square is. Square lattice Hexagonal lattice corner distance = 1.41 × edge distance corner distance = 1.15 × edge distance The same plot count, more evenly spread — which is what a variance estimate rewards.

Minimal reproducible example

Build the lattice in the ownership’s own projected CRS, generate centres row by row with alternate rows offset, and clip to the boundary.

import math

import geopandas as gpd
import numpy as np
from shapely.geometry import Point, Polygon
from shapely.ops import unary_union


def hex_spacing_for_count(area_m2: float, n_plots: int) -> tuple[float, float]:
    """Circumradius R and centre spacing for n hexagons tiling an area."""
    R = math.sqrt(2.0 * area_m2 / (3.0 * math.sqrt(3.0) * n_plots))
    return R, math.sqrt(3.0) * R


def hex_centres(bounds, R: float, rng: np.random.Generator) -> list[Point]:
    """Hexagon centres over a bounding box, with a random origin offset."""
    minx, miny, maxx, maxy = bounds
    dx = math.sqrt(3.0) * R          # horizontal centre spacing
    dy = 1.5 * R                     # vertical spacing between rows
    # Random start so the lattice is not aligned to the bounding box corner.
    x0 = minx - dx + rng.uniform(0, dx)
    y0 = miny - dy + rng.uniform(0, dy)

    pts, row = [], 0
    y = y0
    while y <= maxy + dy:
        offset = (dx / 2.0) if row % 2 else 0.0
        x = x0 + offset
        while x <= maxx + dx:
            pts.append(Point(x, y))
            x += dx
        y += dy
        row += 1
    return pts


def hexagon(centre: Point, R: float) -> Polygon:
    """Pointy-top regular hexagon of circumradius R around a centre."""
    return Polygon([
        (centre.x + R * math.sin(math.pi / 3 * k),
         centre.y + R * math.cos(math.pi / 3 * k))
        for k in range(6)
    ])

Putting them together, clipping to the ownership and keeping only cells whose centre lands on forest:

def hex_plot_grid(boundary_path: str, n_plots: int, seed: int = 20260809,
                  min_overlap: float = 0.5) -> gpd.GeoDataFrame:
    """A hexagonal plot lattice clipped to an ownership boundary."""
    boundary = gpd.read_file(boundary_path)
    if boundary.crs is None or boundary.crs.is_geographic:
        raise ValueError("reproject to a projected, metre-based CRS first")

    area_geom = unary_union(boundary.geometry.values)
    R, spacing = hex_spacing_for_count(area_geom.area, n_plots)
    rng = np.random.default_rng(seed)

    rows = []
    for i, c in enumerate(hex_centres(area_geom.bounds, R, rng)):
        cell = hexagon(c, R)
        inter = cell.intersection(area_geom)
        if inter.is_empty or inter.area < min_overlap * cell.area:
            continue
        rows.append({
            "plot_id": len(rows) + 1,
            "centre_x": c.x,
            "centre_y": c.y,
            "cell_area_ha": cell.area / 10_000.0,
            "inside_frac": inter.area / cell.area,
            "on_boundary": inter.area < 0.999 * cell.area,
            "geometry": cell,
        })
    gdf = gpd.GeoDataFrame(rows, crs=boundary.crs)
    gdf.attrs["circumradius_m"] = R
    gdf.attrs["centre_spacing_m"] = spacing
    gdf.attrs["seed"] = seed
    return gdf

The min_overlap rule is the edge policy, and it needs to be stated rather than defaulted. Keeping every cell that touches the boundary over-samples the edge; keeping only fully interior cells under-samples it and can lose a thin ownership entirely. Half coverage is a defensible middle, and recording inside_frac lets an estimator weight edge plots properly afterwards.

Parameter reference

Parameter Type Default Range Rationale
n_plots int — 30–1000 Drives spacing; set it from a variance target, not from habit
seed int fixed any A recorded seed makes the layout reproducible and auditable
min_overlap float 0.5 0.3–0.7 Edge policy; below 0.3 the sample is dominated by slivers
orientation — pointy-top pointy or flat Rotate the lattice if the ownership has a strong directional grain
R float (m) derived — Circumradius; centre spacing is √3 R
plot radius float (m) 11.3 8–20 The measurement plot inside each cell, unrelated to cell size

Expected output and verification

Three properties define a correct lattice: the plot count is close to target, the nearest-neighbour distances are tight around the intended spacing, and the plots cover the ownership without a systematic hole.

import numpy as np
from scipy.spatial import cKDTree


def audit_grid(gdf, target_n: int) -> dict:
    """Check count, spacing regularity, and coverage of a plot lattice."""
    xy = np.column_stack([gdf["centre_x"], gdf["centre_y"]])
    tree = cKDTree(xy)
    d, _ = tree.query(xy, k=2)          # k=1 is the point itself
    nn = d[:, 1]
    return {
        "n": len(gdf),
        "target": target_n,
        "count_error_pct": round(100.0 * (len(gdf) - target_n) / target_n, 1),
        "nn_mean_m": round(float(nn.mean()), 1),
        "nn_cv": round(float(nn.std() / nn.mean()), 3),
        "edge_plots": int(gdf["on_boundary"].sum()),
    }


report = audit_grid(grid, target_n=200)
assert abs(report["count_error_pct"]) < 15, f"plot count far from target: {report}"
assert report["nn_cv"] < 0.08, f"lattice is irregular — check the row offset: {report}"

The coefficient of variation of nearest-neighbour distance is the sharpest diagnostic. A correct hexagonal lattice has every interior centre at exactly the same distance from its six neighbours, so the CV is dominated by edge cells and stays below about 0.05. A CV above 0.2 means the alternate-row offset is missing and you have built a rectangular grid with extra steps.

Nearest-neighbour distance as a lattice diagnostic Two histograms of nearest-neighbour distance between plot centres. The correct hexagonal lattice produces a single narrow spike at the intended spacing with a small tail from edge plots. The lattice built without the alternate-row offset produces two separate peaks, one at the row spacing and one at the column spacing, which is the signature of a rectangular grid. Correct hexagonal lattice Missing row offset 200 263 340 nearest-neighbour distance (m) CV = 0.04 200 263 340 nearest-neighbour distance (m) CV = 0.24 — two peaks, not one Three edge policies at the ownership boundary A hexagonal lattice meets a curved ownership boundary. Under the keep-everything policy, slivers barely inside the boundary become full plots and the edge is over-sampled. Under the interior-only policy, a whole band of genuine forest along the boundary is never visited. Under the half-coverage policy, cells above fifty percent inside are kept and each records its inside fraction, so an estimator can weight them correctly. keep all touched interior only ≥ 50% inside edge over-sampled boundary forest never visited kept, with inside fraction recorded The policy is a design decision and belongs in the method, not in a default.

Common pitfalls

  • Building the lattice in degrees. Spacing in a geographic CRS is not a distance, and the cells become progressively wider towards the equator. Reproject first.
  • Anchoring the origin to the bounding box. Without a random offset the lattice is deterministic given the extent, which forfeits the randomisation that design-based estimators rely on.
  • No edge policy. Silently keeping every touched cell over-weights the perimeter, which in a fragmented ownership can be most of the sample.
  • Confusing cell size with plot size. The hexagon is the sampling frame; the measurement plot is a small circle at its centre. Setting the plot radius to the hexagon radius means adjacent plots overlap.
  • Regenerating the grid without the seed. A monitoring network that relocates between visits is no longer a time series.

Frequently Asked Questions

Does hexagonal really beat square in practice?

For coverage evenness, measurably: a hexagon’s furthest interior point is 1.15 times its edge distance against 1.41 for a square, so the worst-case gap is smaller. Whether that improves your specific estimator depends on the spatial structure of what you are measuring, but it never makes it worse and it costs nothing.

How do I add a second, denser phase later?

Keep the same origin and seed, halve the circumradius, and the new lattice nests inside the old one — every original centre remains a centre in the finer lattice. That preserves the earlier plots as a subsample and makes the two phases directly comparable.

Should the hexagons be pointy-top or flat-top?

Geometrically identical, just rotated. Choose the orientation whose rows run across the dominant grain of the landscape — perpendicular to contour, planting rows, or drainage — so the lattice does not alias against a periodic feature.

Up: Spatial Plot Sampling Design · Ecological GIS Data Foundations in Python