#Space-Filling Curves for Spatial Partitioning

A columnar file only skips data it can prove is irrelevant. In GeoParquet the proof is per-row-group bounding-box statistics: if a group’s bbox does not overlap the query window, the engine never reads it. That mechanism only pays off when the rows inside each group are spatially compact — and physical row order is what makes them compact. Load features in insertion order and every row group spans the whole dataset extent, so no group is ever skipped. This problem sits at the heart of the Compression, Chunking & Spatial Indexing topic: layout, not the codec, decides how many bytes a bounding-box query has to touch.

Space-filling curves solve the ordering problem by mapping two-dimensional coordinates onto a single axis while preserving locality — points that are close on the map stay close in the 1D key. Sort rows by that key before writing and spatially adjacent features land in the same or neighbouring row groups, giving each group a tight bounding box. This page covers the two curves that matter in production, Hilbert and Z-order (Morton), how to encode centroids to 1D keys, how curve ordering couples with row group sizing, and where it complements tree-based indexes rather than replacing them.

#Prerequisites

• Python 3.10+ with geopandas>=0.14, pyarrow>=14.0, numpy>=1.26, and shapely>=2.0
• A single, consistent CRS across the input — mixed coordinate systems corrupt both the curve keys and the resulting bounding-box statistics
• Familiarity with Parquet internals: row groups, footer column statistics, and predicate pushdown on min/max bounds
• Representative query windows: the actual bounding boxes your workload issues, so you can measure scan-to-return ratio rather than guess
• A partitioning target: a rough sense of whether you serve bounding-box scans, point lookups, or nearest-neighbour queries, since that determines whether curve ordering alone is sufficient

#Architectural Foundations

A space-filling curve visits every cell of a discretised grid exactly once, assigning each cell a sequential integer. The defining property for storage is locality preservation: two cells that are neighbours in 2D usually receive close integers. Sorting rows by that integer therefore groups geographic neighbours into contiguous stretches of the file, and contiguous stretches become row groups with small, non-overlapping bounding boxes.

The two curves used in geospatial systems differ in how faithfully they preserve locality:

• Z-order (Morton) interleaves the bits of the x and y coordinates. It is trivial to compute and its keys are monotonic within any power-of-two quadrant, but it contains long jump discontinuities: when the curve finishes a quadrant it can leap across the entire extent to start the next one, so two rows with adjacent keys are occasionally far apart on the map.
• Hilbert uses a recursive rotation of quadrants so the curve never jumps — consecutive keys are always spatial neighbours. Locality is measurably better, at the cost of a more involved encoder. The Hilbert vs Z-order curve comparison quantifies the gap and shows when Z-order’s cheaper encode wins anyway.

The diagram below shows a Hilbert traversal over a 4×4 grid and how it feeds row-group assignment. Because the curve keeps neighbours together, a bounding-box query over one corner reads a single row group and skips the rest on bbox statistics alone.

Because Parquet stores each column as a contiguous byte stream, the sort also improves compression: neighbouring features share similar coordinate prefixes and repeated attribute values, so a spatially clustered geometry column compresses tighter than a randomly ordered one. Curve ordering is therefore a two-for-one — better skipping and a better ratio from the same pass.

#Step-by-Step Workflow

#1. Normalise CRS and Derive Centroids

The curve encoder needs one comparable coordinate per row. Re-project everything to a single CRS first, then reduce each geometry to a representative point. Centroids are the standard choice: cheap to compute and stable for polygons and lines alike.

# geopandas>=0.14, numpy>=1.26
import geopandas as gpd
import numpy as np

gdf = gpd.read_file("parcels.gpkg")
if gdf.crs is None:
    raise ValueError("Input has no CRS; cannot derive comparable keys")
gdf = gdf.to_crs("EPSG:4326")

cx = gdf.geometry.centroid.x.to_numpy(dtype=np.float64)
cy = gdf.geometry.centroid.y.to_numpy(dtype=np.float64)

For very large or elongated geometries the centroid can fall outside the shape, but for ordering purposes that is harmless — the key only needs to place the feature roughly, since the row group’s actual bbox is computed from full geometry extents.

#2. Quantise Coordinates to an Integer Grid

Curve encoders operate on non-negative integers of a fixed bit width, called the curve order. Scale each axis from its geographic range onto [0, 2**order). An order of 16 gives a 65 536×65 536 lattice — far finer than any realistic row group, so quantisation never limits skipping precision.

# numpy>=1.26
def quantise(cx: np.ndarray, cy: np.ndarray, order: int = 16) -> tuple[np.ndarray, np.ndarray]:
    side = (1 << order) - 1
    nx = ((cx + 180.0) / 360.0 * side).astype(np.uint32)
    ny = ((cy + 90.0) / 180.0 * side).astype(np.uint32)
    return nx, ny

Clamp inputs to the valid coordinate range before scaling if your source data may contain out-of-bounds coordinates; a single stray value at longitude 999 will otherwise compress the useful range into a few cells.

#3. Encode Each Centroid to a 1D Key

Now map the quantised pair to a single integer. Z-order interleaves bits; Hilbert applies a per-level rotation. Both are shown so you can swap them without touching the rest of the pipeline.

# numpy>=1.26
def z_order_key(x: int, y: int, order: int = 16) -> int:
    key = 0
    for i in range(order):
        key |= ((x >> i) & 1) << (2 * i)
        key |= ((y >> i) & 1) << (2 * i + 1)
    return key


def hilbert_key(x: int, y: int, order: int = 16) -> int:
    rx = ry = d = 0
    s = 1 << (order - 1)
    while s > 0:
        rx = 1 if (x & s) > 0 else 0
        ry = 1 if (y & s) > 0 else 0
        d += s * s * ((3 * rx) ^ ry)
        if ry == 0:                      # rotate the quadrant
            if rx == 1:
                x = s - 1 - x
                y = s - 1 - y
            x, y = y, x
        s >>= 1
    return d

Choosing between them is the single biggest quality decision here, and it is covered in depth in Hilbert vs Z-order for spatial partitioning.

#4. Sort Rows and Size Row Groups

Sort the frame by the key, then choose a row group size so each group covers a compact extent. This is where curve ordering couples with the row group sizing strategy for Parquet: the curve sets adjacency, the group size sets how tight each bbox statistic is.

# geopandas>=0.14, numpy>=1.26
keys = np.array([hilbert_key(int(x), int(y)) for x, y in zip(*quantise(cx, cy))])
gdf = gdf.iloc[np.argsort(keys)].reset_index(drop=True)
gdf.to_parquet(
    "parcels_hilbert.parquet",
    compression="zstd",
    compression_level=3,
    row_group_size=100_000,      # tune so each group covers a compact extent
    write_covering_bbox=True,    # embed per-row-group bbox statistics
)

#5. Verify Row-Group Skipping

Ordering is only useful if the engine actually prunes groups. Read the footer statistics back and confirm each group’s bounding box is narrow relative to the dataset extent; a query engine such as DuckDB will then skip non-overlapping groups. Pair aggressive skipping with a modest ZSTD compression level so the groups you do read decompress quickly.

#Partitioning Strategy: From Sort Order to File Layout

Curve ordering operates at two granularities, and production systems use both together. Inside a single file, the sort clusters rows so that per-row-group bounding-box statistics enable skipping. Across many files, the same 1D key becomes a partition boundary: you split the sorted stream into files whose names encode their curve-key range, and a query planner prunes whole files before opening any footer.

The mechanism is a direct extension of the in-file sort. Compute every centroid’s Hilbert or Z-order key, sort, then cut the ordered stream at fixed key intervals — one file per interval. Because the key preserves locality, each file covers a compact, roughly square region of the map. A directory listing then acts as a coarse spatial index: a query with a bounding box computes the key range its window spans and reads only the files overlapping that range. This is cheaper than opening footers because it happens from object-storage metadata alone, and it composes cleanly with the in-file row-group skipping that runs afterwards.

Two parameters govern the strategy. The number of files trades planner overhead against pruning precision: too few files and each spans a wide extent that most queries touch; too many and you pay per-object request latency and metadata bloat on cloud storage. A common target is files of 128–512 MB compressed, sized so a typical query reads a handful rather than hundreds. The key interval should be uniform in data volume, not in key value — dense hotspots would otherwise produce a few enormous files. Split at quantile boundaries of the key distribution so every file holds a similar row count, which keeps parallel readers balanced.

File-level curve partitioning also sets up the tree indexes cleanly. When each file already covers a compact extent, an R-tree over the file bounding boxes is small and shallow, and a quadtree partition aligns naturally with the curve’s quadrant structure. The curve gives you cheap coarse pruning; the tree gives you fine random access to individual features within the surviving files. Combine them rather than choosing between them.

#Production-Ready Implementation

The function below wraps the full workflow: CRS enforcement, centroid quantisation, pluggable curve encoding, curve sort, and a bbox-aware write. It is the minimum viable spatial ordering pass for a production write path.

# Requires: geopandas>=0.14, pyarrow>=14.0, numpy>=1.26, shapely>=2.0
from __future__ import annotations

import time
from pathlib import Path
from typing import Callable, Literal

import geopandas as gpd
import numpy as np


def _z_order(x: int, y: int, order: int) -> int:
    key = 0
    for i in range(order):
        key |= ((x >> i) & 1) << (2 * i)
        key |= ((y >> i) & 1) << (2 * i + 1)
    return key


def _hilbert(x: int, y: int, order: int) -> int:
    d = 0
    s = 1 << (order - 1)
    while s > 0:
        rx = 1 if (x & s) > 0 else 0
        ry = 1 if (y & s) > 0 else 0
        d += s * s * ((3 * rx) ^ ry)
        if ry == 0:
            if rx == 1:
                x = s - 1 - x
                y = s - 1 - y
            x, y = y, x
        s >>= 1
    return d


_ENCODERS: dict[str, Callable[[int, int, int], int]] = {
    "hilbert": _hilbert,
    "zorder": _z_order,
}


def write_curve_ordered_geoparquet(
    source: Path,
    output: Path,
    *,
    curve: Literal["hilbert", "zorder"] = "hilbert",
    crs: str = "EPSG:4326",
    order: int = 16,
    zstd_level: int = 3,
    row_group_size: int = 100_000,
) -> dict[str, object]:
    """Sort a vector dataset by a space-filling curve and write GeoParquet.

    Args:
        source: Input vector file (Shapefile, GeoJSON, GeoPackage, FlatGeobuf).
        output: Destination .parquet path.
        curve: Ordering curve — 'hilbert' for best locality, 'zorder' for
            cheapest encode.
        crs: Target CRS applied to every geometry before keying.
        order: Curve bit-order; 16 yields a 65_536 x 65_536 lattice.
        zstd_level: ZSTD level applied to all columns.
        row_group_size: Rows per Parquet row group.

    Returns:
        Dict with row count, row group count, file size, and elapsed seconds.
    """
    if curve not in _ENCODERS:
        raise ValueError(f"Unknown curve {curve!r}; expected one of {list(_ENCODERS)}")

    t0 = time.perf_counter()
    gdf: gpd.GeoDataFrame = gpd.read_file(source)

    if gdf.crs is None:
        raise ValueError(f"Input has no CRS; cannot key safely: {source}")
    if gdf.crs.to_epsg() != int(crs.split(":")[1]):
        gdf = gdf.to_crs(crs)

    gdf = gdf[gdf.geometry.notna() & gdf.geometry.is_valid].copy()
    if gdf.empty:
        raise ValueError("No valid geometries remain after cleaning")

    cx = gdf.geometry.centroid.x.to_numpy(dtype=np.float64)
    cy = gdf.geometry.centroid.y.to_numpy(dtype=np.float64)
    side = (1 << order) - 1
    nx = np.clip((cx + 180.0) / 360.0 * side, 0, side).astype(np.uint32)
    ny = np.clip((cy + 90.0) / 180.0 * side, 0, side).astype(np.uint32)

    encode = _ENCODERS[curve]
    keys = np.fromiter(
        (encode(int(a), int(b), order) for a, b in zip(nx, ny)),
        dtype=np.uint64,
        count=len(gdf),
    )
    gdf = gdf.iloc[np.argsort(keys)].reset_index(drop=True)

    gdf.to_parquet(
        output,
        compression="zstd",
        compression_level=zstd_level,
        row_group_size=row_group_size,
        write_covering_bbox=True,
    )

    import pyarrow.parquet as pq
    meta = pq.read_metadata(output)
    return {
        "rows": len(gdf),
        "row_groups": meta.num_row_groups,
        "file_bytes": Path(output).stat().st_size,
        "curve": curve,
        "elapsed_s": round(time.perf_counter() - t0, 2),
    }

#Reference: Ordering Strategies Compared

Representative behaviour on a 20-million-row parcel layer written to GeoParquet on a 4-core cloud VM. Scan fraction is the share of row groups a typical municipal-scale bounding-box query must read.

For point-in-polygon and nearest-neighbour access, curve ordering alone is not enough — pair it with spatial partitioning using quadtree indexes, and see quadtree vs R-tree index build time for how those trees are built and maintained.

#Failure Modes and Gotchas

#Frequently Asked Questions

Do I need a Hilbert curve, or is a simple sort by latitude enough?

A one-dimensional sort by latitude clusters rows along a single axis, so a bounding-box query that is narrow in longitude still touches almost every row group. A space-filling curve interleaves both axes, keeping features that are close in two dimensions close in the 1D key. In practice, Hilbert ordering cuts the number of row groups scanned by a bounding-box query by three to ten times compared with a single-axis sort on the same data.

How does curve ordering interact with row group size?

The curve decides which rows are adjacent; the row group size decides how many adjacent rows share one bounding-box statistic. Smaller groups produce tighter bounding boxes and more precise skipping, but they add footer overhead and can under-fill the ZSTD window. Target 32 to 128 MB per group and tune it together with the sort, never in isolation.

Does a space-filling curve replace a quadtree or R-tree?

No. Curve ordering is a physical row layout that makes coarse bounding-box pruning effective inside a columnar file. Quadtree and R-tree structures are explicit indexes that accelerate point-in-polygon and nearest-neighbour lookups. They are complementary: sort rows by a curve for scan-heavy analytics, and build a tree index when you need fine-grained random access to individual features.

#Related

• Hilbert vs Z-Order Curve for Spatial Partitioning
• Quadtree vs R-Tree Index Build Time
• Spatial Partitioning with Quadtree Indexes
• Row Group Sizing Strategies for Parquet
• ZSTD Compression Levels for Geospatial Data

← Back to Compression, Chunking & Spatial Indexing