Generating Hierarchical Administrative Boundaries

A synthetic boundary set has three levels — region, district, ward — and a spatial join between them returns wards that straddle two districts, districts whose area does not equal the sum of their wards, and a handful of wards belonging to no district at all.

Part of Polygon Tessellation Algorithms: a single tessellation has to partition its window, and a hierarchy has to do that at every level and have the levels agree with each other. The second requirement is the one that is usually discovered late.

Root Cause: Independent Generation Cannot Nest

The natural implementation generates each level separately — tessellate the extent into regions, tessellate it again into districts, again into wards — and then hopes the levels align. They never do, and the reason is arithmetic rather than tolerance.

Three independent tessellations of the same window produce three independent sets of edges. A district boundary and a ward boundary that ought to be the same line are two different lines, computed from different seeds, agreeing nowhere except by accident. Snapping them together afterwards is the sliver problem at every level boundary simultaneously, and it does not converge: snapping wards to districts moves ward edges, which breaks the ward-level partition, which needs re-snapping.

The resolution is to invert the construction. Generate the finest level only, and derive every coarser level by aggregating it. Then a district is by definition the union of its wards, its boundary is by definition made of ward edges, and the containment relationships hold by construction rather than by check.

Three hierarchy construction strategies against four nesting invariants A matrix with three construction strategies as rows and four invariants as columns. The first strategy generates each level independently: it satisfies the per-level partition invariant, because each level was tessellated properly, but it fails containment, area conservation and shared-edge geometry, because the levels were computed from different seeds and their edges agree nowhere. The second strategy generates the coarsest level first and subdivides downward: it satisfies containment and shared edges, since every child is cut from its parent, but the cutting introduces slivers along parent boundaries so area conservation holds only to a tolerance rather than exactly. The third generates the finest level once and unions upward: all four invariants hold by construction, because a parent is literally the union of its children and its boundary is literally made of their edges. The note underneath draws the conclusion: the third strategy is the only one where the checks are assertions rather than repairs, and the other two require a reconciliation pass that does not converge. Only one of the three makes the invariants true by construction strategy level partitions child in one parent areas conserve shared edges independent per level yes each level valid no straddles + orphans no gaps at every level no different seeds top-down subdivision yes cut from parent yes by construction to tolerance slivers on cuts yes parent edge reused bottom-up aggregation yes one tessellation yes by construction exactly union of children yes child edges only Bottom-up is the only strategy where the containment and area checks are assertions rather than repairs. The other two need a reconciliation pass, and that pass does not converge: snapping children to a parent edge moves child edges, which breaks the child-level partition, which needs snapping again.
Three construction strategies against the four invariants a hierarchy must satisfy — only bottom-up aggregation satisfies all of them by construction rather than by repair.

Prerequisite Check: State the Invariants Before Building Anything

python
def containment(child_gdf, parent_gdf, child_key: str, parent_key: str) -> list[str]:
    """Every child must lie inside exactly one parent."""
    joined = child_gdf.sjoin(parent_gdf, predicate="within", how="left")
    orphans = joined[joined[parent_key].isna()][child_key].tolist()
    straddlers = (joined.groupby(child_key)[parent_key].nunique()
                  .loc[lambda s: s > 1].index.tolist())
    return [f"orphan: {o}" for o in orphans] + [f"straddles: {s}" for s in straddlers]


def area_conservation(child_gdf, parent_gdf, parent_key: str, tol=1e-9) -> list[str]:
    """A parent's area must equal the sum of its children's."""
    summed = child_gdf.groupby(parent_key).geometry.apply(lambda g: g.area.sum())
    out = []
    for key, parent_area in parent_gdf.set_index(parent_key).geometry.area.items():
        got = summed.get(key, 0.0)
        if abs(got - parent_area) / max(parent_area, 1e-12) > tol:
            out.append(f"{key}: children sum {got:.4f} vs parent {parent_area:.4f}")
    return out

Both of these should return empty on a correctly constructed hierarchy, and both are cheap. Writing them before the generator is worth doing, because they are also the specification: a construction that cannot satisfy them is a construction to abandon rather than to debug.

The tolerance in area_conservation is deliberately near zero. On a bottom-up hierarchy the parent is literally the union of the children, so the areas agree to floating-point precision; a tolerance loose enough to accommodate a top-down construction is a tolerance that hides real gaps.

Fix: Generate the Finest Level, Aggregate Upward

1 — Tessellate once, at the finest level

python
from shapely import set_precision

GRID = 1e-7


def finest_level(seeds, envelope, grid: float = GRID):
    """One tessellation, snapped to the precision grid, covering the envelope exactly."""
    cells = voronoi_clipped(seeds, envelope)
    return [set_precision(c, grid) for c in cells]

Everything above this level is derived, so this is the only place the partition invariants have to be established — and it is the only place they can be violated.

2 — Aggregate with a declared assignment, not a spatial join

python
from shapely.ops import unary_union


def aggregate(children, assignment: dict, key_of) -> dict:
    """Union children by a declared parent key.

    The assignment is data — a lookup from child key to parent key — not a spatial
    predicate. A spatial join here would re-introduce the ambiguity the whole
    construction exists to avoid.
    """
    groups: dict = {}
    for child in children:
        groups.setdefault(assignment[key_of(child)], []).append(child)
    return {parent: unary_union(members) for parent, members in groups.items()}

The distinction in that docstring is the one that matters. Deriving the parent by unioning declared members gives a boundary made of exactly the children’s edges. Deriving it by a spatial join against an independently generated parent geometry gives back every problem the construction was meant to remove.

3 — Generate the assignment so the groups are spatially coherent

Children have to be grouped into contiguous, plausibly-shaped parents rather than arbitrary sets. Region-growing over the adjacency graph does this and is straightforward:

python
def grow_regions(adjacency: dict, keys: list, n_parents: int, target_size: int, rng) -> dict:
    """Grow contiguous groups from seeds until each reaches its target size."""
    seeds = [keys[rng.integers(len(keys))] for _ in range(n_parents)]
    assignment = {s: f"P{i:03d}" for i, s in enumerate(seeds)}
    frontier = {f"P{i:03d}": set(adjacency[s]) for i, s in enumerate(seeds)}
    sizes = {f"P{i:03d}": 1 for i in range(n_parents)}
    while any(f for f in frontier.values()):
        for parent in sorted(frontier):
            if sizes[parent] >= target_size or not frontier[parent]:
                continue
            candidates = sorted(k for k in frontier[parent] if k not in assignment)
            if not candidates:
                frontier[parent] = set()
                continue
            pick = candidates[rng.integers(len(candidates))]
            assignment[pick] = parent
            sizes[parent] += 1
            frontier[parent] |= set(adjacency[pick]) - set(assignment)
    # anything unassigned joins the neighbouring parent it shares the most edge with
    for k in keys:
        if k not in assignment:
            assignment[k] = dominant_neighbour_parent(k, adjacency, assignment)
    return assignment

The sorted calls are not cosmetic. Iterating a set of candidate cells in arbitrary order makes the assignment depend on hash ordering, which makes the whole hierarchy non-reproducible — the same defect as everywhere else in this area, arriving through a different door.

Bottom-up aggregation from fine cells to parent boundaries Three panels read left to right. The first shows the finest level: a grid of small cells that partitions the extent exactly, produced by a single tessellation. The second shows the same cells shaded by the contiguous group each was grown into, starting from four seed cells and expanding over the adjacency graph until every cell belongs to exactly one group. The third shows the result of unioning each group: four parent polygons whose outlines follow cell edges precisely, because they are cell edges. Nothing was snapped, cut or reconciled between the panels. The note underneath states what this buys: containment and area conservation hold to floating-point precision rather than to a tolerance, so the verification checks are assertions about a construction rather than measurements of how badly a repair went. One tessellation, grown into groups, unioned into parents 1 — finest level, one tessellation 2 — cells grown into contiguous groups 3 — groups unioned into parents Nothing is snapped, cut or reconciled between the panels. The parent outlines in the third panel are cell edges from the first, so containment and area conservation hold to floating-point precision — the verification checks become assertions about a construction rather than measurements of how badly a repair went.
The construction in one direction: cells are grown into contiguous groups, groups are unioned into parents, and the parents' edges are by definition cell edges.

Verification Step: Assert Nesting at Every Level Pair

python
def test_each_level_partitions_the_envelope(levels, envelope, tol=1e-9):
    for name, gdf in levels.items():
        union = gdf.geometry.union_all()
        assert union.symmetric_difference(envelope).area < tol, f"{name} is not a partition"


def test_every_child_nests_in_one_parent(levels, hierarchy):
    for child, parent in hierarchy:
        issues = containment(levels[child], levels[parent], f"{child}_id", f"{parent}_id")
        assert not issues, f"{child} → {parent}: {issues[:5]}"


def test_areas_conserve_up_the_hierarchy(levels, hierarchy):
    for child, parent in hierarchy:
        issues = area_conservation(levels[child], levels[parent], f"{parent}_id")
        assert not issues, issues[:5]


def test_parents_are_contiguous(levels, adjacency, hierarchy):
    """A region-grown parent should be one piece; a disjoint one is a bug in the growth."""
    for _, parent in hierarchy:
        multi = levels[parent][levels[parent].geometry.geom_type == "MultiPolygon"]
        assert multi.empty, f"{parent} has {len(multi)} disjoint parents"
Straddled child area under independent generation versus aggregation The horizontal axis is the number of parent units the extent is divided into, from six up to thirty, with the child level held fixed at ninety units. The vertical axis is the percentage of each child's area that falls outside the parent it is assigned to — the parent containing the majority of it — averaged over all children. When the two levels are tessellated independently from different seeds, this rises steadily as the parent count grows, because more parent boundaries means more opportunities for one to cut through the middle of a child. Even at the coarsest division, a noticeable share of child area is on the wrong side of a parent line. The flat line along the bottom is bottom-up aggregation, which is exactly zero at every parent count, because a parent is the union of its children and cannot cut through one. The note underneath makes the practical point: the straddle rate is not a tolerance to tighten, it is a structural consequence of generating the levels apart, and no snapping tolerance drives it to zero. Independently generated levels straddle more as the hierarchy gets finer 6 12 20 30 0% 10% 20% 30% parent units the extent is divided into child area outside its parent (%) 8.5% 13.5% 18.4% 24.6% bottom-up aggregation — exactly 0% at every level count independent tessellation per level finest level unioned upward 90 child units on the unit square against 6–30 independently seeded parents, measured on a 100×100 sample grid with a fixed seed. The straddle rate is not a tolerance to tighten. It is a structural consequence of generating the levels apart, and no snapping tolerance drives it to zero — snapping only moves the error into slivers.
Child area falling outside its assigned parent as the hierarchy gets finer: independent generation drifts steadily upward, aggregation stays at exactly zero.

The contiguity test is the one that catches a defect the others do not. A parent assembled from cells that are not adjacent is geometrically valid, nests correctly and conserves area — and it is a district in two pieces on opposite sides of the region, which no administrative geography has and which any consumer will notice.

Edge Cases & Gotchas

Genuinely disjoint administrative units. They exist — enclaves, offshore islands, historical anomalies. If the release models them, the contiguity assertion needs an allow-list rather than a relaxation, because relaxing it hides the accidental cases too.

Levels that do not nest in reality. Postal geographies and administrative geographies famously cross. If the release contains both, they are two separate hierarchies that share a finest level, not one hierarchy with five levels — and modelling them as the latter guarantees failures that are not defects.

Population-weighted targets. Administrative units are usually sized by population rather than by area, so the growth target should be a population sum drawn from the covariate rather than a cell count. That is a one-line change to the growth loop and it makes the output far more plausible.

The finest level is too fine to aggregate quickly. Region growing over a million cells is slow. Growing over a coarser intermediate and then refining is faster and equivalent, provided the intermediate is itself derived from the finest level rather than generated independently — which is the same rule again.

Choosing the Finest Level

Because everything is derived from it, the finest level is the only real design decision in the construction, and it is worth making deliberately rather than defaulting to whatever the first consumer asked for.

It has to be at least as fine as the finest level anyone will publish. Aggregation only goes one way. A hierarchy built on wards can produce districts and regions, and can never produce anything below a ward — so if a consumer later needs output areas, the whole set has to be regenerated from a finer base, and every previously published boundary changes with it.

It should be finer than that, but not by much. A base an order of magnitude below the finest published level gives room to add a level later without regenerating, and costs only storage. A base three orders below multiplies the region-growing cost for no benefit anyone will ever observe, and makes the adjacency graph large enough that the growth pass becomes the slowest step in the pipeline.

Its cell count sets the granularity of every parent’s size. A parent grown from twenty cells can only ever have twenty possible sizes, so a population target expressed to three significant figures cannot be hit. If the release is calibrated to population totals, the base needs enough cells per parent that the growth can land close to a target rather than stepping past it.

The practical shape of this is that the base level is a published artifact in its own right, not an implementation detail, and it deserves the same version identity as the levels above it — a consumer who wants to build their own aggregation needs it, and a producer who wants to add a level in a year’s time needs to know exactly which base the existing levels came from.