#Compression, Chunking & Spatial Indexing: The Foundation of Cloud-Native Geospatial Storage

For GIS data engineers, Python backend developers, and cloud architects, the bottleneck in modern spatial analytics is rarely raw compute. It is I/O inefficiency: the wrong codec, oversized or undersized row groups, and storage layouts that force query engines to scan entire files for localized bounding-box queries. When working with terabytes of vector features, raster tiles, or point clouds, traditional monolithic file formats collapse under their own weight.

The architectural solution is a tightly coupled triad — compression, chunking, and spatial indexing — applied together so that each layer amplifies the others. This guide covers the engineering decisions behind each layer, the formats and tools that implement them, production Python patterns, and the trade-offs that determine which configuration fits your workload.

#Four Sub-Problems Every Spatial Storage Layer Must Solve

Geospatial datasets introduce four engineering challenges that interact with each other. Getting one wrong degrades the other three.

Geometry encoding. Raw WKT strings are verbose and non-numeric, blocking vectorized operations. The GeoParquet columnar format, defined by the OGC GeoParquet specification, standardizes Well-Known Binary (WKB) geometry encoding within Parquet column chunks, enabling CPU-vectorized predicate evaluation without a geometry parser in the hot path. Encoding choice feeds directly into compression efficiency: WKB is 30–50% smaller than equivalent WKT before any codec is applied.

Coordinate reference system (CRS) management. Mixed CRS within a dataset produces incorrect bounding-box statistics, breaks spatial join predicates, and violates the GeoParquet spec. CRS must be normalized to a single system — typically WGS 84 / EPSG:4326 for global data — before spatial indexing is applied, because the index’s correctness depends on coordinate comparability.

Spatial indexing. Traditional B-tree indexes fail on two-dimensional coordinates. Space-filling curves map 2D/3D coordinates to a 1D integer sequence that preserves locality, so nearby geographic features land in the same or adjacent chunks. Spatial partitioning with quadtree indexes provides hierarchical space decomposition on top of this ordering, enabling metadata-driven chunk skipping before any decompression occurs.

Compression and layout. Codec selection, compression level, and row group sizing strategies for Parquet determine how many bytes cross the network per query and how many CPU cycles decode them. These parameters must be co-designed: a row group sized for bulk analytics but indexed for point lookups wastes both I/O and memory.

#How Hilbert Curve Sorting Enables Chunk Skipping

The hardest concept to internalize in cloud-native GIS is why physical row order in a Parquet file determines query speed. The following diagram shows the mechanism: a 4×4 grid of geographic cells (representing a small region) is traversed by a Hilbert curve, producing a 1D sequence number for each cell. Features in cells 0–3 land in row group A; features in cells 4–7 land in row group B; and so on. A bounding-box query covering only the top-left corner touches row group A exclusively — the other groups are skipped by reading footer statistics alone.

The key insight: without Hilbert sorting, a bounding-box query covering the top-left quadrant would force the engine to read all 16 cells to find which rows match. With Hilbert sorting, the engine reads the footer bounding-box statistics for each row group, confirms that groups B–D lie entirely outside the query window, and reads only group A. Scan-to-return ratio drops from 16:1 to 1:1 for this example.

#Format and Technology Landscape

Selection criteria: Use GeoParquet when your primary access pattern is analytical SQL or DataFrame operations on tens of millions of features. Use FlatGeobuf when you serve individual features to a web mapping client and need sub-10 ms random access — the GeoParquet vs FlatGeobuf performance comparison quantifies where each format wins. Use Cloud-Optimized GeoTIFF for raster tile pipelines where HTTP range requests must resolve to a single object without intermediate tiling.

#Architecture: The Medallion Storage Pipeline

The medallion pattern organizes spatial data through three physical layers, each with distinct compression and indexing requirements. The diagram below shows how a single source feed flows through bronze, silver, and gold, with the codec, chunk size, and index strategy tightening at each stage to match how the data is read downstream.

Bronze layer — raw ingestion. Ingest source files (Shapefile, GeoJSON, raw CSV with coordinates) exactly as received. Apply only CRS normalization and geometry validation. Do not compress at this stage; you need maximum read flexibility for schema inspection and re-processing.

Silver layer — cleaned, partitioned. Convert to GeoParquet with ZSTD level 3, explicit row group sizes tuned to analytical workloads (64–128 MB uncompressed), and Hilbert-sorted geometry. Apply dictionary encoding for categorical GIS attributes on land-use codes, sensor classifications, and administrative identifiers. Embed CRS metadata per the GeoParquet spec. This layer feeds all query engines.

Gold layer — query-optimized views. Materialize pre-aggregated or pre-filtered subsets for specific dashboards or services. Row group size shrinks to 8–32 MB for high-selectivity point lookups. Spatial partitioning aligns with the application’s bounding boxes (e.g., administrative regions, map tile grids). Compression levels rise to 5–6 for cold archival views where read frequency is low.

A single Python pipeline manages all three layers:

# Requires: pyarrow>=14.0, geopandas>=0.14, pyproj>=3.6
import geopandas as gpd
import pyarrow as pa
import pyarrow.parquet as pq
from pathlib import Path
from typing import Literal


def write_silver_layer(
    source_path: Path,
    dest_path: Path,
    target_crs: str = "EPSG:4326",
    zstd_level: int = 3,
    row_group_mb: int = 64,
) -> dict[str, object]:
    """
    Convert a source vector file to a silver-layer GeoParquet file.

    Returns a dict with row count, compressed size, and elapsed seconds.
    """
    import time
    import os

    t0 = time.perf_counter()

    gdf = gpd.read_file(source_path)

    # Enforce single CRS — required before spatial indexing
    if gdf.crs is None:
        raise ValueError(f"Source has no CRS: {source_path}")
    if str(gdf.crs) != target_crs:
        gdf = gdf.to_crs(target_crs)

    # Drop null geometries
    gdf = gdf[gdf.geometry.notna()].copy()

    # Hilbert sort for spatial locality
    gdf = gdf.sort_values(by=gdf.geometry.name, key=lambda col: col.apply(
        lambda g: _hilbert_key(g.centroid.x, g.centroid.y)
    )).reset_index(drop=True)

    # Estimate rows per group from uncompressed size
    mem_bytes = gdf.memory_usage(deep=True).sum()
    rows_total = len(gdf)
    target_bytes = row_group_mb * 1024 * 1024
    row_group_size = max(1000, int(rows_total * target_bytes / max(mem_bytes, 1)))

    # In production: use geopandas.io.arrow for full GeoParquet 1.0 compliance
    gdf.to_parquet(
        dest_path,
        compression="zstd",
        compression_level=zstd_level,
        row_group_size=row_group_size,
        schema_version="1.0",
        write_covering_bbox=True,
    )

    elapsed = time.perf_counter() - t0
    compressed_bytes = os.path.getsize(dest_path)

    return {
        "rows": rows_total,
        "compressed_bytes": compressed_bytes,
        "elapsed_s": round(elapsed, 2),
        "row_group_size": row_group_size,
    }


def _hilbert_key(x: float, y: float) -> int:
    """Map WGS-84 lon/lat to a 32-bit Hilbert index."""
    # Normalize to [0, 2^16)
    nx = int((x + 180.0) / 360.0 * 65535)
    ny = int((y + 90.0) / 180.0 * 65535)
    return _xy_to_hilbert(nx, ny, 16)


def _xy_to_hilbert(x: int, y: int, order: int) -> int:
    d = 0
    s = 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)
        x, y = _rotate(s, x, y, rx, ry)
        s >>= 1
    return d


def _rotate(n: int, x: int, y: int, rx: int, ry: int) -> tuple[int, int]:
    if ry == 0:
        if rx == 1:
            x = n - 1 - x
            y = n - 1 - y
        return y, x
    return x, y

#Compression Strategies for Geospatial Payloads

Geospatial data exhibits distinct statistical properties: coordinate arrays are highly correlated (sequential precision-decimal-degree values), categorical attributes repeat frequently (land-use codes, sensor types, administrative boundaries), and geometry bytes compress better than random binary. Generic compression algorithms waste cycles on these predictable patterns.

ZSTD has become the default codec for cloud-native geospatial formats. Unlike gzip or Snappy, ZSTD supports dictionary pre-training, fine-grained level tuning, and multi-threaded decode. Selecting the appropriate ZSTD compression levels for geospatial data directly determines query latency and storage cost: levels 1–3 prioritize decompression throughput (>1 GB/s per core), levels 4–9 increase ratio at mild CPU cost, and levels 10+ are suitable only for cold archival where write latency is unconstrained.

Dictionary encoding handles categorical columns before the byte-level codec sees them. By mapping repeated string values to integer IDs, dictionary encoding for categorical GIS attributes reduces payload size by 60–90% for land-cover codes, sensor classifications, and administrative boundaries. PyArrow applies this automatically when column cardinality is below 50% of row count, but you should verify it is active with pq.read_schema(path).field('column_name').metadata.

Per-column codec selection. Parquet permits different codecs per column. Apply ZSTD to geometry (WKB) and numeric coordinate columns. Apply dictionary encoding to categoricals before ZSTD. For timestamp columns, use DELTA_BINARY_PACKED encoding (built into Parquet) before ZSTD to exploit temporal monotonicity in sensor and satellite data ingestion.

#Performance and Trade-off Analysis

#Compression Ratio vs Query Latency

#Indexing Strategy Comparison

#Decision Flowchart

#Production Implementation

The pattern below demonstrates a complete chunked write pipeline with Hilbert sorting, per-column compression, and PyArrow metadata validation — the minimum viable implementation for a silver-layer write in production.

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

import time
from pathlib import Path

import geopandas as gpd
import numpy as np
import pyarrow as pa
import pyarrow.parquet as pq
from shapely import to_wkb


def build_geoparquet_silver(
    source: Path,
    output: Path,
    *,
    crs: str = "EPSG:4326",
    zstd_level: int = 3,
    row_group_rows: int = 131_072,  # ~64 MB uncompressed for typical vector
    write_stats: bool = True,
) -> dict[str, object]:
    """
    Read a vector file, normalize CRS, Hilbert-sort, and write a GeoParquet
    silver-layer file with per-column compression and embedded bbox statistics.

    Parameters
    ----------
    source : Path
        Input vector file (Shapefile, GeoJSON, FlatGeobuf, GeoPackage).
    output : Path
        Destination .parquet file path.
    crs : str
        Target CRS for all geometry. Default WGS 84.
    zstd_level : int
        ZSTD compression level applied to all columns. 1–3 for hot paths.
    row_group_rows : int
        Number of rows per Parquet row group.
    write_stats : bool
        Whether to embed per-column min/max statistics in the footer.

    Returns
    -------
    dict with 'rows', 'row_groups', 'file_bytes', 'elapsed_s'.
    """
    t0 = time.perf_counter()

    gdf: gpd.GeoDataFrame = gpd.read_file(source)

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

    # Drop null/invalid geometries before indexing
    valid_mask = gdf.geometry.notna() & gdf.geometry.is_valid
    dropped = len(gdf) - valid_mask.sum()
    if dropped:
        print(f"[warn] Dropped {dropped} null/invalid geometries")
    gdf = gdf[valid_mask].copy()

    # Hilbert sort: map centroid lon/lat to 32-bit key
    cx = gdf.geometry.centroid.x.to_numpy(dtype=np.float64)
    cy = gdf.geometry.centroid.y.to_numpy(dtype=np.float64)
    hkeys = _batch_hilbert(cx, cy, order=16)
    gdf = gdf.iloc[np.argsort(hkeys)].reset_index(drop=True)

    # Convert geometry to WKB bytes for Parquet storage
    wkb_array = pa.array(to_wkb(gdf.geometry.values), type=pa.large_binary())

    # Build PyArrow table from non-geometry columns
    non_geom = [c for c in gdf.columns if c != gdf.geometry.name]
    tbl: pa.Table = pa.Table.from_pandas(gdf[non_geom], preserve_index=False)
    tbl = tbl.append_column(
        pa.field(gdf.geometry.name, pa.large_binary()),
        wkb_array,
    )

    # GeoParquet 1.0 geometry column metadata
    geo_meta = {
        "version": "1.0",
        "primary_column": gdf.geometry.name,
        "columns": {
            gdf.geometry.name: {
                "encoding": "WKB",
                "crs": gdf.crs.to_wkt(),
                "bbox": list(gdf.total_bounds),
            }
        },
    }
    import json
    existing_meta = tbl.schema.metadata or {}
    tbl = tbl.replace_schema_metadata({
        **existing_meta,
        b"geo": json.dumps(geo_meta).encode(),
    })

    parquet_opts = pq.ParquetWriter(
        str(output),
        tbl.schema,
        compression="zstd",
        compression_level=zstd_level,
        write_statistics=write_stats,
        data_page_version="2.0",
    )

    n_groups = 0
    for start in range(0, len(tbl), row_group_rows):
        parquet_opts.write_table(tbl.slice(start, row_group_rows))
        n_groups += 1
    parquet_opts.close()

    elapsed = round(time.perf_counter() - t0, 2)
    file_bytes = output.stat().st_size

    return {
        "rows": len(gdf),
        "row_groups": n_groups,
        "file_bytes": file_bytes,
        "elapsed_s": elapsed,
    }


def _batch_hilbert(lons: np.ndarray, lats: np.ndarray, order: int) -> np.ndarray:
    """Vectorized Hilbert index for WGS-84 lon/lat arrays."""
    scale = (1 << order) - 1
    nx = ((lons + 180.0) / 360.0 * scale).astype(np.int32).clip(0, scale)
    ny = ((lats + 90.0) / 180.0 * scale).astype(np.int32).clip(0, scale)
    return np.array([_xy_to_hilbert_scalar(int(x), int(y), order) for x, y in zip(nx, ny)], dtype=np.int64)


def _xy_to_hilbert_scalar(x: int, y: int, order: int) -> int:
    d = 0
    s = 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

#Schema Governance and Compliance

CRS standardization. Establish a single organization-wide default CRS and enforce it at ingestion time, not at query time. Re-projecting at query time doubles I/O for cross-dataset joins and breaks bounding-box filter statistics. Use pyproj.CRS.from_user_input() to parse both EPSG codes and WKT strings from legacy sources.

Schema evolution. Parquet is schema-on-read, which sounds permissive but creates silent errors when upstream sources add or rename columns. Version your GeoParquet schemas using a schema_version key in file-level Parquet metadata. Pin PyArrow schemas explicitly in your pipeline config; never infer schema from the first batch of a streaming source.

Access control. Cloud object storage IAM policies should grant query engines read access at the prefix level (e.g., s3://bucket/silver/vectors/) rather than per-file. This enables partition-aligned access control — a spatial region’s data is isolated in its own prefix — without per-file policy maintenance. Avoid presigned URLs for bulk analytical access; they bypass audit logging.

Specification compliance. The OGC GeoParquet specification requires that the geo metadata key in the Parquet file footer contain at minimum: version, primary_column, columns[name].encoding, and columns[name].crs. Validate compliance with geopandas.io.arrow.geopandas_to_arrow() and inspect output with pq.read_schema(path).metadata[b'geo'].

#Strategic Selection Framework

#Operational Resilience and Maintenance

Cloud-native spatial storage introduces failure modes absent from traditional file-based GIS. Corrupted row groups, mismatched spatial metadata, and codec incompatibilities can silently degrade query accuracy without raising exceptions.

Automated integrity validation. After every write, verify that embedded bbox values in the geo metadata match gdf.total_bounds on the written data. Check row group counts against expected output using pq.ParquetFile(path).metadata.num_row_groups. Any mismatch indicates a truncated write or serialization error.

Three production metrics to monitor:

• Scan-to-return ratio. Bytes read from storage divided by bytes in the result set. A ratio above 10:1 indicates poor spatial sorting or row group sizes that don’t match query selectivity. Target below 5:1 for dashboard workloads.
• Decompression latency p99. CPU time spent in ZSTD decode per query. A rising p99 with stable query volume suggests dictionary drift in categorical columns — retrain or rebuild affected files.
• Chunk cache hit rate. Fraction of row groups served from in-process cache (Arrow BufferReader) versus network. Below 60% for a recurring dashboard indicates the cache capacity is insufficient for the working set.

Maintenance cycles. Rebuild spatial sorting and update bbox footer statistics when the dataset grows by more than 20% through appends — Hilbert locality degrades with unsorted appends. Retrain ZSTD dictionaries for categorical columns quarterly when attribute vocabularies evolve (e.g., new land-cover classification codes from updated surveys).

#Frequently Asked Questions

When should I keep Shapefiles instead of migrating to GeoParquet?

Keep Shapefiles only when interoperability with desktop GIS software that cannot read GeoParquet is a hard requirement, or when datasets are smaller than 50 MB and are never queried programmatically. For any cloud pipeline, analytical query, or dataset over 100 MB, GeoParquet with ZSTD compression outperforms Shapefiles on every measurable dimension: file size, query latency, and schema flexibility. The limitations of Shapefiles in modern data stacks are well-documented — the 2 GB per-.dbf file cap alone is disqualifying for most production datasets.

What row group size should I use?

Target 32–64 MB uncompressed for spatially filtered workloads, and 128–512 MB for bulk aggregations. Smaller groups increase selectivity precision — fewer false-positive rows are read after the bounding-box filter — but add file-footer overhead per group. Benchmark your actual query bounding boxes against candidate group sizes before committing. See row group sizing strategies for Parquet for a full benchmark methodology.

Does ZSTD compression hurt query latency?

At levels 1–3, ZSTD decompression throughput exceeds 1 GB/s per core on modern hardware, which outpaces the read bandwidth of S3 Standard and GCS Standard. Latency impact is negligible for remote reads. It only becomes measurable when processing locally cached data on NVMe storage, where Snappy’s simpler algorithm has a slight edge on raw decompression speed with no network bottleneck.

Can I query GeoParquet files directly in S3 with DuckDB?

Yes. Run INSTALL httpfs; LOAD httpfs; INSTALL spatial; LOAD spatial;, then execute SELECT * FROM read_parquet('s3://bucket/prefix/*.parquet') WHERE ST_Within(geometry, ST_GeomFromText('...')). DuckDB reads Parquet footers first, applies predicate pushdown on min/max bounding-box statistics per row group, and transfers only relevant chunks. For best performance, files must be Hilbert-sorted before upload.

What is the difference between a quadtree index and a Hilbert curve sort?

A quadtree is a hierarchical tree structure that partitions space into nested quadrants, used for point-in-polygon and range queries with a persistent external index file. A Hilbert curve is a space-filling curve that maps 2D coordinates to a 1D integer, used to physically sort rows in a Parquet file so geographically nearby features reside in adjacent chunks. In GeoParquet, Hilbert sorting is the indexing primitive embedded in the file itself. Spatial partitioning with quadtree indexes describes both approaches and when to combine them.

How do I handle mixed CRS datasets in a single pipeline?

Never mix CRS within a single GeoParquet file. Re-project all source data to a single target CRS before writing. Mixed CRS produces incorrect bounding-box statistics embedded in the file footer — which query engines use for predicate pushdown — causing queries to silently return wrong or incomplete results. Use pyproj.CRS.from_user_input() and GeoPandas .to_crs() to standardize upstream in your pipeline.

#Related

• ZSTD Compression Levels for Geospatial Data — codec level selection by workload and data type
• Row Group Sizing Strategies for Parquet — benchmark-driven chunk size selection
• Dictionary Encoding for Categorical GIS Attributes — attribute compression for high-cardinality columns
• Spatial Partitioning with Quadtree Indexes — hierarchical spatial partitioning strategies
• GeoParquet vs FlatGeobuf Performance — when analytical columnar storage beats streaming feature access
• Geospatial Storage Fundamentals & Format Comparison — format selection across the full GIS storage landscape

← Back to Home