#Spatial Partitioning with Quadtree Indexes

Geospatial workloads at scale demand deterministic query performance, predictable memory footprints, and storage layouts that align with cloud object I/O patterns. The decisive question is whether a bounding-box filter or spatial join can skip irrelevant data at the storage layer — before any geometry is decoded — and that depends entirely on how the dataset is partitioned. A quadtree answers it by recursively subdividing a two-dimensional bounding box into four equal quadrants until a configurable capacity threshold or depth limit is reached, producing partition edges that are known and addressable at write time.

This guide belongs to the Compression, Chunking & Spatial Indexing for Cloud-Native GIS collection, and walks the full workflow from coordinate-system preparation through production tuning. The partitioning decision comes first; the codec and file-layout choices — ZSTD compression levels for geospatial data and row-group sizing strategies for Parquet — are tuned once the quadrant boundaries are fixed.

#Prerequisites

Validate these baseline requirements before building a production quadtree pipeline:

• Python 3.9+ with numpy>=1.24, shapely>=2.0, pyarrow>=12.0, pyproj>=3.5, and zstandard>=0.20.0
• CRS normalization: all geometries must share a single planar coordinate reference system (EPSG:3857 for web-mapping contexts, or a local metric projection to minimize area distortion). Mixed CRS inputs cause quadrant boundary miscalculations that are invisible until query results diverge.
• Validated bounding box: a tight (minx, miny, maxx, maxy) envelope fully containing the dataset. Outliers must be clipped or rejected before ingestion — a single out-of-bounds feature will trigger unbounded recursion.
• Storage target: Parquet or Arrow IPC files on S3-compatible object storage with row-group-level statistics enabled for min/max predicate pushdown.
• Monitoring: tracemalloc or psutil for heap profiling during recursive construction; time.perf_counter for per-partition throughput measurement.

#Architectural Foundations

A quadtree operates on a recursive invariant: a root node covers the entire dataset envelope. When a node accumulates more features than max_features_per_leaf, it divides into four children — NW, NE, SW, SE — and redistributes its geometries. This repeats depth-first until every leaf satisfies the capacity constraint or the tree reaches max_depth.

The deterministic structure of this subdivision delivers three architectural advantages over dynamic index types:

Predictable I/O boundaries. Partition edges are known at write time, so a bounding-box query translates directly into an exact list of leaf files without traversing the tree at query time. Query engines like DuckDB or AWS Athena can prune irrelevant partitions purely from directory metadata.

Cache-aligned reads. Fixed quadrant extents align with CPU cache lines and cloud storage block sizes (typically 4 KB–1 MB), reducing partial-read penalties when engines fetch contiguous byte ranges within a leaf file.

Embarrassingly parallel construction. Once the root splits into its first four quadrants, those branches are fully independent and can be filled concurrently across worker threads or distributed compute nodes.

Where quadtrees underperform is in skewed distributions: a metropolitan area surrounded by rural empty space will produce one deep subtree and three shallow, near-empty quadrants. For datasets with extreme geometric heterogeneity or high polygon overlap, ZSTD compression levels for geospatial data and row-group sizing strategies for Parquet remain critical regardless of index type — but the partitioning decision itself may favor a Hilbert-curve sort or R-tree instead.

The diagram below shows the end-to-end relationship: raw point data enters the quadtree builder, the root node recursively splits into leaf quadrants, and each leaf is serialized as a ZSTD-compressed Parquet file with embedded row-group statistics that enable predicate pushdown in cloud query engines.

#Step-by-Step Workflow

#1. Normalize CRS and Validate the Bounding Box

Project all geometries to a metric planar CRS using pyproj. Validate that no feature centroid falls outside the global envelope before tree construction begins. A strict bounding check at this stage prevents degenerate recursion and ensures every geometry routes to a deterministic quadrant.

# pyproj>=3.5, shapely>=2.0
from __future__ import annotations
import logging
from typing import List, Tuple

from pyproj import Transformer
from shapely.geometry import Point, box

logger = logging.getLogger(__name__)

BBox = Tuple[float, float, float, float]  # (minx, miny, maxx, maxy)


def normalize_to_epsg3857(
    points: List[Tuple[float, float]],
    source_epsg: int,
) -> List[Point]:
    """Re-project (lon, lat) or projected coords to EPSG:3857 Points."""
    transformer = Transformer.from_crs(
        f"EPSG:{source_epsg}", "EPSG:3857", always_xy=True
    )
    result: List[Point] = []
    for x, y in points:
        px, py = transformer.transform(x, y)
        result.append(Point(px, py))
    return result


def validate_envelope(points: List[Point], envelope: BBox) -> List[Point]:
    """Reject points outside *envelope*; warn on any exclusions."""
    region = box(*envelope)
    valid = [p for p in points if region.contains(p)]
    dropped = len(points) - len(valid)
    if dropped:
        logger.warning("Dropped %d out-of-envelope points before indexing", dropped)
    return valid

#2. Initialize the Root Node and Configure Capacity Thresholds

Instantiate a root node covering the full dataset envelope. Two parameters govern the entire tree structure:

• max_features_per_leaf: triggers subdivision when exceeded. Cloud workloads typically perform well at 100–500 features; lower values improve selectivity, higher values reduce file count.
• max_depth: prevents runaway recursion in extreme density clusters. Values of 8–12 cover virtually all real-world spatial distributions.

The interaction between these two parameters determines partition granularity. A dataset with 10 million points at max_features_per_leaf=200 and max_depth=10 produces at most 4^10 ≈ 1 million possible leaves, but in practice far fewer because most quadrants exhaust the data long before maximum depth.

#3. Insert Features and Trigger Recursive Splits

Route each point to its containing quadrant. When a leaf exceeds max_features_per_leaf and has not reached max_depth, split it into four children and redistribute its accumulated features.

# shapely>=2.0, Python 3.9+
from __future__ import annotations
from dataclasses import dataclass, field
from typing import Iterator, List, Tuple

from shapely.geometry import Point, box as shapely_box

BBox = Tuple[float, float, float, float]


@dataclass
class QuadTreeNode:
    bounds: BBox
    depth: int = 0
    points: List[Point] = field(default_factory=list)
    children: List["QuadTreeNode"] = field(default_factory=list)

    def insert(
        self,
        point: Point,
        max_features: int,
        max_depth: int,
    ) -> None:
        if not shapely_box(*self.bounds).contains(point):
            return
        if not self.children:
            self.points.append(point)
            if len(self.points) > max_features and self.depth < max_depth:
                self._split(max_features, max_depth)
        else:
            for child in self.children:
                child.insert(point, max_features, max_depth)

    def _split(self, max_features: int, max_depth: int) -> None:
        minx, miny, maxx, maxy = self.bounds
        midx = (minx + maxx) / 2.0
        midy = (miny + maxy) / 2.0
        quadrants: List[BBox] = [
            (minx, midy, midx, maxy),  # NW
            (midx, midy, maxx, maxy),  # NE
            (minx, miny, midx, midy),  # SW
            (midx, miny, maxx, midy),  # SE
        ]
        self.children = [
            QuadTreeNode(bounds=q, depth=self.depth + 1) for q in quadrants
        ]
        for p in self.points:
            for child in self.children:
                child.insert(p, max_features, max_depth)
        self.points.clear()

    def leaves(self) -> Iterator["QuadTreeNode"]:
        """Yield all non-empty leaf nodes."""
        if not self.children:
            if self.points:
                yield self
        else:
            for child in self.children:
                yield from child.leaves()

For polygonal or line geometries, replace contains with intersects when routing features to quadrants. Track geometries that fall into more than one quadrant and decide whether to duplicate them or route them to a parent overflow bucket (see the failure modes section below).

#4. Serialize Leaf Partitions to Columnar Layout

Convert each leaf node into a discrete Parquet file. Store the quadrant bounding box as metadata columns alongside geometry coordinates so downstream engines can filter at the row-group level using min/max statistics. Keeping a single row group per leaf file means the engine’s min/max filter covers the entire file, maximizing skip efficiency — see row-group sizing strategies for Parquet for detailed guidance on this trade-off. If your features carry repeated categorical attributes — land-use codes, administrative names, sensor classes — pair this layout with dictionary encoding for categorical GIS attributes so those columns shrink independently of the coordinate columns.

# pyarrow>=12.0, numpy>=1.24
from __future__ import annotations
import hashlib
from pathlib import Path
from typing import TYPE_CHECKING

import numpy as np
import pyarrow as pa
import pyarrow.parquet as pq

if TYPE_CHECKING:
    from .quadtree import QuadTreeNode  # local import for typing only


def stable_quad_id(node: "QuadTreeNode") -> str:
    """Return a deterministic, human-readable ID for this leaf."""
    coords = f"d{node.depth}_{node.bounds[0]:.6f}_{node.bounds[1]:.6f}"
    short_hash = hashlib.sha1(coords.encode()).hexdigest()[:8]
    return f"q{node.depth}_{short_hash}"


def serialize_leaf(
    node: "QuadTreeNode",
    output_dir: Path,
    zstd_level: int = 3,
) -> Path:
    """
    Write a leaf node's points to a ZSTD-compressed Parquet file.

    Row-group statistics (min/max x, y) are written automatically by
    PyArrow and enable predicate pushdown in DuckDB, Athena, and Trino.
    """
    if not node.points:
        raise ValueError(
            f"Leaf node at depth {node.depth} is empty — prune before serializing"
        )

    coords = np.array([(p.x, p.y) for p in node.points], dtype=np.float64)
    quad_id = stable_quad_id(node)

    table = pa.table(
        {
            "x": pa.array(coords[:, 0], type=pa.float64()),
            "y": pa.array(coords[:, 1], type=pa.float64()),
            "quad_id": pa.array([quad_id] * len(coords), type=pa.string()),
            "bbox_minx": pa.array([node.bounds[0]] * len(coords), type=pa.float64()),
            "bbox_miny": pa.array([node.bounds[1]] * len(coords), type=pa.float64()),
            "bbox_maxx": pa.array([node.bounds[2]] * len(coords), type=pa.float64()),
            "bbox_maxy": pa.array([node.bounds[3]] * len(coords), type=pa.float64()),
        }
    )

    output_path = output_dir / f"{quad_id}.parquet"
    pq.write_table(
        table,
        output_path,
        compression="zstd",
        compression_level=zstd_level,
        # Keep one row group per leaf so min/max stats cover the full file
        row_group_size=len(node.points),
    )
    return output_path

For codec selection guidance, including the spatial correlation argument for choosing ZSTD compression levels for geospatial data over alternatives, see that dedicated page.

#5. Benchmark and Validate Predicate Pushdown

After writing leaf files, confirm that query engines skip irrelevant partitions. The fastest validation path uses DuckDB’s spatial extension with EXPLAIN ANALYZE:

# duckdb>=0.9.0 — install with: pip install duckdb
import duckdb

con = duckdb.connect()
con.execute("INSTALL spatial; LOAD spatial;")

# Point the engine at the partition directory
con.execute("""
    CREATE VIEW spatial_parts AS
    SELECT * FROM read_parquet('output_dir/*.parquet', hive_partitioning=false)
""")

# Bounding-box query — check 'Parquet Files' rows_read in EXPLAIN output
result = con.execute("""
    EXPLAIN ANALYZE
    SELECT COUNT(*) FROM spatial_parts
    WHERE x BETWEEN 1_800_000 AND 1_900_000
      AND y BETWEEN 6_100_000 AND 6_200_000
""").fetchall()

for row in result:
    print(row[1])
# Look for "Parquet Files": files_read << total_files to confirm skipping

If files_read equals total_files, your row-group statistics are not being used. Likely causes: the query engine lacks pushdown support for the chosen file layout, or row groups span too many distinct x/y ranges (fix by keeping one row group per leaf).

#Production-Ready Full Pipeline

# numpy>=1.24, shapely>=2.0, pyarrow>=12.0, pyproj>=3.5, zstandard>=0.20.0
from __future__ import annotations
import logging
import tracemalloc
from pathlib import Path
from typing import List, Tuple

from shapely.geometry import Point

from .normalize import normalize_to_epsg3857, validate_envelope
from .quadtree import QuadTreeNode
from .serialize import serialize_leaf

logger = logging.getLogger(__name__)

BBox = Tuple[float, float, float, float]


def build_quadtree_partitions(
    raw_points: List[Tuple[float, float]],
    source_epsg: int,
    envelope: BBox,
    output_dir: Path,
    max_features_per_leaf: int = 200,
    max_depth: int = 10,
    zstd_level: int = 3,
) -> dict[str, int | float]:
    """
    Full pipeline: normalize → index → serialize → return stats.

    Returns a dict with leaf_count, total_points, and peak_memory_mb.
    """
    output_dir.mkdir(parents=True, exist_ok=True)
    tracemalloc.start()

    # 1. Normalize CRS
    points: List[Point] = normalize_to_epsg3857(raw_points, source_epsg)
    # 2. Validate bounding box
    points = validate_envelope(points, envelope)

    # 3. Build quadtree
    root = QuadTreeNode(bounds=envelope, depth=0)
    for pt in points:
        root.insert(pt, max_features=max_features_per_leaf, max_depth=max_depth)

    # 4. Serialize non-empty leaves
    leaf_count = 0
    for leaf in root.leaves():
        try:
            path = serialize_leaf(leaf, output_dir, zstd_level=zstd_level)
            leaf_count += 1
            logger.debug("Wrote %s (%d points)", path.name, len(leaf.points))
        except ValueError as exc:
            logger.warning("Skipping empty leaf: %s", exc)

    _, peak_mem = tracemalloc.get_traced_memory()
    tracemalloc.stop()

    return {
        "leaf_count": leaf_count,
        "total_points": len(points),
        "peak_memory_mb": round(peak_mem / 1_048_576, 2),
    }

#Benchmark Reference Matrix

The table below shows measured outcomes for a 10 M-point uniform distribution across a 1,000 km × 1,000 km EPSG:3857 envelope, written to a local SSD then queried with DuckDB 0.10.

max_features=200, max_depth=10 is the recommended starting point for interactive analytics on cloud storage. Values below 100 fragment storage and can overwhelm S3’s per-prefix request quota during full scans. Target leaf sizes that produce Parquet files of 64 MB–256 MB after compression — this aligns with S3 multipart thresholds and DuckDB’s parallel read scheduler.

#Failure Modes and Gotchas

CRS mismatch before indexing. Mixing geographic (latitude/longitude in degrees) with projected (meters) coordinates causes quadrant boundaries to misalign by orders of magnitude. A bounding box in meters wrapped around geographic coordinates produces a tree where every point falls into the same leaf. Always validate CRS before calling insert.

Skewed point distributions. Dense urban clusters trigger deep subtrees while surrounding rural quadrants remain shallow. When leaf imbalance exceeds 5:1 (deepest leaf vs. shallowest leaf at similar depth), consider a hybrid strategy: apply a Hilbert-curve sort and fixed-size chunk approach for the dense area, and use the quadtree for the sparse remainder.

Empty quadrants consuming metadata overhead. After splitting, three quadrants may be empty while one absorbs all features. Prune zero-feature leaves before serialization with the if node.points guard in leaves(). Unpruned empty Parquet files still consume list entries in the directory manifest and add unnecessary overhead to partition discovery.

Geometry boundary overlap. Polygons and linestrings crossing quadrant edges must either be duplicated or stored in a parent overflow bucket. Duplication is safe for features with a unique ID that consumers can deduplicate; overflow bucketing requires all queries to always include the overflow partition. Track the duplication ratio: if more than 10 % of features are duplicated, the geometry extent is too large relative to leaf size and max_features should be increased.

Unbounded recursion from identical coordinates. Perfectly coincident points grouped at the same coordinate will exhaust max_depth immediately. The max_depth guard prevents infinite recursion, but the resulting leaves will still exceed max_features. Detect and deduplicate exact-coordinate duplicates upstream.

ZSTD level mismatch with column type. Coordinate columns (x, y) compress well at level 3–5 because values within a leaf are spatially correlated. Applying level 12+ to coordinate columns adds substantial CPU overhead with negligible ratio improvement. See ZSTD compression levels for geospatial data for level selection guidance.

#FAQ

When should I use a quadtree instead of an R-tree for spatial indexing?

Quadtrees work best for uniformly distributed point datasets and web-mapping tile generation, where fixed grid boundaries align naturally with display tiles and cloud storage block sizes. R-trees outperform quadtrees when geometries are irregular or heavily overlapping — administrative polygons, road networks — because R-trees fit tight bounding envelopes around each feature rather than imposing a fixed grid. Choose a quadtree when deterministic partition boundaries and parallelizable construction matter; choose an R-tree when geometry heterogeneity demands flexible envelope fitting.

How do I choose max_features_per_leaf and max_depth?

Start with max_features_per_leaf between 100 and 500 and max_depth between 8 and 12, then benchmark read latency and partition file count. Target leaf sizes that produce Parquet files of 64 MB to 256 MB after compression — this aligns with S3 multipart thresholds and DuckDB’s parallel read scheduler. Smaller leaves improve selectivity for tight bounding-box queries but fragment storage and inflate metadata overhead.

What ZSTD compression level should I apply to quadtree leaf partitions?

Level 3–5 is the correct range for hot quadtree partitions. Coordinate arrays within a single quadrant are spatially correlated, so the compressor achieves high ratios without deep search. Level 3 minimizes decompression latency for interactive query engines; level 5 is appropriate for archival partitions. Avoid levels above 7 for hot partitions — the CPU cost rarely justifies the marginal ratio improvement over already-correlated spatial data.

How do I handle geometries that span multiple quadrant boundaries?

The two standard strategies are duplication and overflow bucketing. Duplication copies the geometry into every intersecting quadrant — simple but inflates storage and risks result duplication unless consumers deduplicate on a unique feature ID. Overflow bucketing stores boundary-spanning geometries in a parent node or dedicated overflow partition and always queries it alongside leaf partitions. Overflow bucketing is preferable for large polygons where duplication overhead is high; duplication is acceptable for point and small-polygon datasets where overlap is rare.

#Related

• ZSTD Compression Levels for Geospatial Data — codec selection and level tuning for coordinate and raster columns
• Row-Group Sizing Strategies for Parquet — aligning row group size with quadtree leaf capacity for predicate pushdown
• Dictionary Encoding for Categorical GIS Attributes — compressing repeated string attributes alongside spatial coordinates
• Comparing GeoParquet vs FlatGeobuf Performance — choosing the right output format for partitioned spatial data

← Back to Compression, Chunking & Spatial Indexing for Cloud-Native GIS