Building a quadtree for spatial data indexing

Spatial datasets attach meaning to coordinates: delivery stops belong to addresses, tree samples have locations, and council assets occupy parcels. A full scan is simple but costly as records grow. A quadtree repeatedly divides a two-dimensional region into four smaller regions, letting a lookup ignore most records.

The payoff is narrowing space before examining records. Points are cleanest because each follows one path to a leaf bucket. Bounding boxes and polygons need extra rules when geometry crosses a division. A sound implementation therefore specifies coordinate conventions, boundary behaviour, depth limits, and query semantics rather than hiding them in a library.

Australian workloads make these decisions concrete. Sydney and Melbourne delivery and rideshare systems search nearby drivers, while utilities inspect assets over broad territories. Brisbane and Perth add sparse distributions that can create deep trees. Council datasets and the Australian Geospatial National Address Framework mix addresses, parcels, and service boundaries.

This guide builds a working point index, then explains how to adapt it for objects with extents. It covers an in-memory structure for prototypes, services, and algorithm exercises, with complexity analysis and transferable tests. The broader hello ML community publishes related explanations for the supporting programming ideas.

Why a quadtree suits spatial workloads

A node stores a rectangular region, a small record bucket, and up to four children. A search starts at the root and descends only into regions that can contain results. The structure works well for “everything in this viewport” queries, where spatial locality allows broad areas to be pruned.

Choice of index matters. A k-d tree partitions points along alternating axes and handles nearest-neighbour work well, but does not naturally represent overlapping rectangles. An R-tree groups bounding boxes and often suits polygons and lines. A quadtree is attractive for two-dimensional data, moderate updates, and teams wanting predictable control over memory and latency.

Defining points, bounds, and node invariants

Start with immutable Point and Bounds records. A point holds coordinates and application data; bounds define an axis-aligned rectangle. Keep the coordinate reference system outside these geometry classes, but validate minimum and maximum values so malformed extents cannot corrupt traversal.

from dataclasses import dataclass
from typing import Any


@dataclass(frozen=True)
class Point:
    x: float
    y: float
    data: Any = None


@dataclass(frozen=True)
class Bounds:
    min_x: float
    min_y: float
    max_x: float
    max_y: float

    def __post_init__(self):
        if self.min_x > self.max_x or self.min_y > self.max_y:
            raise ValueError("invalid bounds")

    def contains(self, point: Point) -> bool:
        return (
            self.min_x <= point.x <= self.max_x
            and self.min_y <= point.y <= self.max_y
        )

    def intersects(self, other: "Bounds") -> bool:
        return not (
            self.max_x < other.min_x
            or other.max_x < self.min_x
            or self.max_y < other.min_y
            or other.max_y < self.min_y
        )

A node also needs its bounds, points, child map, depth, and bucket capacity. Its central invariant is ownership: after subdivision, each point belongs to exactly one child. This implementation sends boundary points west or south at each split, making that ownership deterministic.

Closed comparisons count a point on a viewport edge as a match, which is useful for map selection boxes. That convention should be tested explicitly; applications needing half-open regions can change the comparisons without changing the tree’s shape.

Implementing insertion and subdivision

Insertion starts at the root. If a child quadrant contains the point, the node delegates to it; otherwise it stores the point locally. When a bucket exceeds capacity, the node splits, creates four children, and redistributes its points. A maximum depth prevents pathological data from producing an unbounded hierarchy.

The code uses Python programming features such as dataclasses, type hints, and a dictionary-backed child collection. The child map keeps quadrant lookup explicit, while a fixed quadrant order makes subdivision deterministic.

class QuadNode:
    def __init__(
        self,
        bounds: Bounds,
        capacity: int = 4,
        max_depth: int = 12,
        depth: int = 0,
    ):
        if capacity < 1 or max_depth < 0:
            raise ValueError("invalid tree limits")

        self.bounds = bounds
        self.capacity = capacity
        self.max_depth = max_depth
        self.depth = depth
        self.points: list[Point] = []
        self.children: dict[int, "QuadNode"] = {}

    def _quadrant(self, point: Point) -> int:
        mid_x = (self.bounds.min_x + self.bounds.max_x) / 2
        mid_y = (self.bounds.min_y + self.bounds.max_y) / 2

        if point.x <= mid_x:
            return 0 if point.y >= mid_y else 1
        return 2 if point.y >= mid_y else 3

    def insert(self, point: Point) -> None:
        if not self.bounds.contains(point):
            raise ValueError("point lies outside node bounds")

        quadrant = self._quadrant(point)
        child = self.children.get(quadrant)

        if child is not None:
            child.insert(point)
            return

        self.points.append(point)
        if len(self.points) > self.capacity and self.depth < self.max_depth:
            self._subdivide()

    def _subdivide(self) -> None:
        min_x, min_y = self.bounds.min_x, self.bounds.min_y
        max_x, max_y = self.bounds.max_x, self.bounds.max_y
        mid_x = (min_x + max_x) / 2
        mid_y = (min_y + max_y) / 2

        self.children = {
            0: QuadNode(
                Bounds(min_x, mid_y, max_x, max_y),
                self.capacity,
                self.max_depth,
                self.depth + 1,
            ),
            1: QuadNode(
                Bounds(min_x, min_y, mid_x, mid_y),
                self.capacity,
                self.max_depth,
                self.depth + 1,
            ),
            2: QuadNode(
                Bounds(mid_x, mid_y, max_x, max_y),
                self.capacity,
                self.max_depth,
                self.depth + 1,
            ),
            3: QuadNode(
                Bounds(mid_x, min_y, max_x, mid_y),
                self.capacity,
                self.max_depth,
                self.depth + 1,
            ),
        }

        pending, self.points = self.points, []
        for point in pending:
            quadrant = self._quadrant(point)
            self.children[quadrant].insert(point)


class Quadtree:
    def __init__(
        self,
        bounds: Bounds,
        capacity: int = 4,
        max_depth: int = 12,
    ):
        self.root = QuadNode(bounds, capacity, max_depth)

    def insert(self, point: Point) -> None:
        self.root.insert(point)

This remains a point quadtree, not a hidden rectangle index. For a box record, descend only when one child wholly contains it; otherwise retain it at the current node. Rectangle queries must test the true extent, preventing a record that crosses a split from disappearing from a neighbouring quadrant.

Answering rectangle and point queries

A range query visits the root and prunes nodes whose bounds do not intersect the requested rectangle. Surviving leaf points are tested against the query. An explicit stack avoids Python’s recursion limit during traversal.

def query_range(root: QuadNode, area: Bounds) -> list[Point]:
    matches: list[Point] = []
    stack = [root]

    while stack:
        node = stack.pop()
        if not node.bounds.intersects(area):
            continue

        matches.extend(
            point for point in node.points
            if area.contains(point)
        )
        stack.extend(node.children.values())

    return matches


def query_point(root: QuadNode, x: float, y: float) -> list[Point]:
    area = Bounds(x, y, x, y)
    return [
        point for point in query_range(root, area)
        if point.x == x and point.y == y
    ]


tree = Quadtree(Bounds(0, 0, 100, 100))
tree.insert(Point(20, 30, "depot"))
tree.insert(Point(80, 70, "vehicle"))

hits = query_range(tree.root, Bounds(0, 0, 50, 50))

For polygons, parcels, or segments, replace the stored point with a record containing a box and test its extent instead. A record straddling quadrants can stay at their common ancestor or be inserted into every intersecting child; the latter duplicates data, while the former can concentrate large features near the root. An R-tree is often safer for this pattern.

Exact nearest-neighbour search needs another traversal. The branch nearest the query may not contain the nearest point, so use best-first traversal with a priority queue ordered by minimum squared distance to each cell, pruning cells whose bound exceeds the best candidate. This suits delivery fleets in Melbourne and Sydney, where road distance may still be calculated separately.

Measuring performance and planning production use

For a balanced point tree, insertion is typically O(log n), while a range query costs returned records plus visited boundary cells. Both have O(n) worst cases: dense clusters can force repeated splits, and a root-sized query inspects everything. Memory is O(n) plus node overhead.

Real data breaks synthetic assumptions. Duplicate addresses, exact ties, boundary points, and sudden clusters affect occupancy. Measure p50, p95, and p99 latency, allocation rate, and memory on representative data; capacity four is a starting point. Deletions need invariant-preserving removal or tombstones plus rebuilding, while concurrent writes require locking or snapshots.

Coordinate systems need care in Australia. GDA2020 underpins much of the national spatial framework, state services often use projected GDA2020/MGA zones such as 55 or 56, and web maps commonly display WGS84. Mislabelling coordinates can scatter records across distant cells. A compact service can reproduce the design using C programming, but explicit ownership, bounds checks, and profiling matter more than language choice.

Before owning a custom index, compare it with a spatial database or maintained R-tree library. Australian councils, utilities, logistics firms, and environmental-monitoring programs often need durability and standard query tools as well as speed. A custom structure earns its cost when latency targets or embedded constraints justify the maintenance burden.

A production-ready adoption checklist

The design is ready for a controlled trial when its data model, invariants, and failure behaviour are explicit. Use the following checks before moving beyond a prototype:

• Fix the spatial model. Record the coordinate reference system, axis order, units, allowed precision, and whether queries use inclusive or exclusive edges.
• Set structural limits. Choose bucket capacity and maximum depth from representative distributions, and define duplicate-coordinate behaviour at the depth limit.
• Prove correctness. Compare results with a brute-force scan using empty windows, skewed clusters, duplicate points, boundary points, rectangles, and split-crossing records.
• Benchmark real workloads. Track p50, p95, and p99 latency, throughput, allocations, memory, and rebuild time using sanitised Sydney, Melbourne, Brisbane, or council datasets.
• Plan operations. Specify persistence, concurrency, deletion, tombstones, compaction, versioning, and when an R-tree or spatial database becomes preferable.

A focused quadtree can give predictable response times for Australian fleet, asset, and environmental datasets without a large service. Its strongest case is stable bounds, frequent locality-based queries, and a team prepared to test edge cases. Under those conditions, the hierarchy turns spatial locality into less work per request.