#Geospatial Storage Fundamentals & Format Comparison

Modern geospatial pipelines have outgrown legacy file formats and monolithic databases. As spatial datasets scale into terabytes and real-time analytics become standard, platform teams and GIS data engineers must treat geometry as first-class data within cloud-native storage architectures. Four interconnected problems drive every storage decision: how coordinates are encoded into bytes, how projections are tracked across systems, how spatial indexes accelerate queries, and how physical layout is aligned to access patterns. Getting any one of these wrong cascades into degraded query performance, incorrect spatial joins, or brittle schema evolution.

This guide — one of the core reference tracks linked from the data-compression.net home page — establishes the foundational principles, evaluates modern format trade-offs, and provides actionable implementation patterns for production systems, from raw ingestion to query-optimised gold layers. Teams building the ingest side of these pipelines should pair this material with the data conversion and migration pipelines reference, which covers batch conversion patterns, schema mapping, and metadata preservation in depth; teams tuning the read side should pair it with the compression, chunking and spatial indexing track.

#1. Core Concept Map: Four Interconnected Sub-problems

Every geospatial storage decision reduces to four inseparable engineering constraints. Solving them in isolation introduces hidden failure modes; solving them together produces pipelines that scale predictably.

#Geometry Encoding

Geometry encoding converts coordinates and topological relationships into bytes. Well-Known Binary (WKB) is the production standard: compact, unambiguous, and supported by every major GIS library. Well-Known Text (WKT) is reserved for debugging and human-readable interchange. Production systems store WKB internally and transcode to WKT only at API boundaries. The encoding choice directly affects compression efficiency — WKB’s fixed-width coordinate sequences compress more predictably than JSON-encoded coordinate arrays.

#CRS Management

Spatial data is semantically void without projection context. Modern formats embed EPSG codes or PROJ strings alongside geometry columns, enabling automatic transformation during ingestion. The most common failure mode is implicit CRS assumptions: two datasets in different projections loaded into the same index produce silently incorrect spatial joins. Enforce explicit CRS validation at the silver layer of every ingestion pipeline, before any spatial index is built. CRS errors are most disruptive during schema mapping from legacy formats, where source files rarely carry authoritative projection metadata.

#Spatial Indexing

Sequential scans are prohibitive for polygon-heavy workloads. Index structures — R-trees, quadtrees, and space-filling curves — enable bounding-box filtering and spatial joins without full dataset reads. In distributed environments, global indexing strategies such as spatial partitioning with quadtree indexes are critical for co-locating spatially adjacent records in the same storage partitions, minimising cross-partition shuffle during joins.

#Compression and Layout

Physical byte arrangement governs I/O efficiency. Columnar layouts isolate geometry and attribute columns, enabling predicate pushdown and compressed reads that skip irrelevant data. Row-oriented or interleaved layouts suit streaming and web delivery where entire feature records are needed at once. ZSTD compression levels configured per column type — higher levels for coordinate columns with repetitive deltas, lower for densely packed polygon arrays — can reduce storage footprints by 60–80 % versus uncompressed WKB. Categorical attribute columns (land-use class, road type, jurisdiction code) benefit additionally from dictionary encoding, which can compress these columns by 10–100× before any codec is applied.

#2. Format Landscape: Trade-offs and Selection Criteria

The geospatial ecosystem has evolved from single-purpose interchange formats to specialised, cloud-native alternatives. Each format solves a different access pattern; selecting the wrong one creates bottlenecks that cannot be patched at the query layer.

#Legacy Interchange: Shapefile

Despite its age, the ESRI Shapefile remains ubiquitous in source data from government agencies and legacy GIS vendors. However, its architecture is fundamentally misaligned with distributed systems. As detailed in Shapefile limitations in modern data stacks, the mandatory multi-file structure (.shp, .shx, .dbf, .prj) creates atomicity issues in object storage — partial uploads leave datasets in an unreadable state. The 2 GB file size cap and 10-character field-name restriction break modern schema evolution patterns, and the format carries no native support for null geometry or multi-part feature types.

For new pipelines, treat Shapefiles strictly as ingestion artifacts: land them in the bronze layer and convert immediately upon arrival.

#Web-Native Interchange: GeoJSON

GeoJSON became the de facto standard for web mapping and RESTful APIs due to its JSON-native structure and human readability. While excellent for frontend consumption and lightweight feature exchange, it carries significant overhead for analytical workloads. Text-based coordinate serialisation inflates storage footprints by 3–5× compared to binary equivalents, and JSON parsing overhead scales poorly with large feature collections.

The GeoJSON overhead and serialisation costs analysis quantifies CPU and memory penalties across feature sizes. GeoJSON remains optimal for client-side rendering, API responses under a few hundred features, and configuration-driven geometry, but it should not persist in data lake storage layers.

#Analytical Columnar: GeoParquet

GeoParquet extends the Apache Parquet specification to natively support spatial types, bridging the gap between GIS and big data ecosystems. By storing geometries as WKB within Parquet’s columnar structure, it inherits robust compression (ZSTD, Snappy), schema evolution, and predicate pushdown. The format is governed by the OGC GeoParquet Specification, ensuring cross-platform compatibility across DuckDB, Apache Arrow, and cloud data warehouses.

Understanding Parquet columnar storage for GIS explains how row-group partitioning lets query engines skip irrelevant blocks before deserialising geometry. For teams building analytical data lakes or running heavy spatial aggregations, GeoParquet is the current industry baseline.

#Cloud-Optimised Streaming: FlatGeobuf

FlatGeobuf (FGB) was engineered specifically for HTTP range requests and streaming consumption. Unlike traditional vector formats that require full-file downloads, FGB embeds a packed Hilbert R-tree at the file header, enabling clients to fetch only the bounding boxes or features relevant to the current viewport. This makes it exceptionally fast for web tile servers, mobile applications, and serverless APIs that cannot afford full dataset materialisation.

When evaluating streaming performance against analytical columnar formats, comparing GeoParquet vs FlatGeobuf performance provides benchmark-driven guidance for matching workload characteristics to format strengths.

#Format Selection Criteria Table

#3. Architecture: Medallion Pipeline and Ingestion Patterns

Production geospatial pipelines follow a medallion architecture: raw data lands in a bronze layer untouched, is validated and enriched in silver, and is materialised as query-optimised formats in gold. Each transition enforces a specific set of invariants.

Bronze layer preserves source data exactly as received — no transformation, no re-encoding. This ensures raw data is always reproducible and auditable. Shapefiles, GeoJSON files, and vendor exports land here unchanged. The batch conversion pipeline that moves data from bronze to silver handles format detection, driver selection, and per-file error routing.

Silver layer is where spatial invariants are enforced. ETL steps at this stage include: CRS normalisation (all geometries projected to EPSG:4326 or a project-standard CRS), geometry validation and repair (closing rings, removing self-intersections, filling null geometries with sentinel values), field-name normalisation, and schema enforcement against a registered schema registry. Any row failing spatial validation is routed to a quarantine partition with a structured error code — not silently dropped. The metadata preservation guide details how to carry source attribution, lineage fields, and original CRS strings through this transformation without loss.

Gold layer materialises query-optimised formats for distinct consumer profiles. A single source-of-truth GeoParquet partition table serves analytical engines; FlatGeobuf files are generated on-demand for web APIs. This polyglot approach avoids storage duplication while optimising access paths for different consumers.

#4. Performance and Trade-off Analysis

Storage format selection is only half the equation. Query performance hinges on how spatial data is indexed, compressed, and accessed at runtime.

#Compression Ratio vs Latency Matrix

Geometry columns compress differently than scalar attributes. Coordinate deltas and repeated topology patterns benefit from dictionary encoding and delta compression; dense polygon arrays respond well to ZSTD. The trade-off is CPU decompression cost versus storage and network egress savings.

#Indexing Strategy Comparison

Row-oriented formats rely on external index files or database-managed structures (e.g., PostGIS GiST indexes). Cloud-native formats embed indexes directly into the file layout:

• FlatGeobuf: packed Hilbert R-tree stored at the file header, enabling O(log n) spatial lookups via HTTP Range headers. No external index required; the client fetches only the bytes covering the requested bounding box.
• GeoParquet: row-group-level bounding-box statistics. Query engines use min/max metadata to prune irrelevant row groups before reading geometry columns. Combining row group sizing with Hilbert-curve spatial sorting of features maximises pruning effectiveness.
• PostGIS: GiST R-tree indexes built and maintained by the database engine. Optimal for transactional workloads with frequent inserts and updates, but impractical for read-optimised object storage.
• Distributed partitioning: for datasets exceeding tens of millions of features, partitioning by S2, H3, or Geohash cells distributes spatial locality across storage shards, reducing cross-partition reads during spatial joins.

#Decision Flowchart: Format Selection

#5. Production Implementation

This Python snippet demonstrates the core medallion pipeline pattern: reading a mixed-source bronze input, validating and reprojecting in the silver stage, and writing a spatially sorted GeoParquet gold partition with ZSTD compression.

# Dependencies (pin in requirements.txt):
#   geopandas==1.0.1  pyarrow==16.1.0  pyogrio==0.9.0  shapely==2.0.5

from __future__ import annotations

import logging
from pathlib import Path

import geopandas as gpd
import pyarrow as pa
import pyarrow.parquet as pq
from shapely.validation import make_valid

logger = logging.getLogger(__name__)

TARGET_CRS = "EPSG:4326"
COMPRESSION = "zstd"
COMPRESSION_LEVEL = 3
ROW_GROUP_SIZE = 200_000


def silver_transform(
    source_path: Path,
    quarantine_path: Path,
) -> gpd.GeoDataFrame:
    """
    Load a bronze-layer file, validate CRS, repair geometry,
    and return a clean GeoDataFrame ready for gold-layer writing.
    """
    gdf = gpd.read_file(source_path, engine="pyogrio")

    if gdf.crs is None:
        raise ValueError(f"No CRS detected in {source_path}. Assign an EPSG code at ingest.")

    if gdf.crs.to_epsg() != 4326:
        logger.info("Reprojecting from %s to %s", gdf.crs, TARGET_CRS)
        gdf = gdf.to_crs(TARGET_CRS)

    invalid_mask = ~gdf.geometry.is_valid
    if invalid_mask.any():
        n_invalid = invalid_mask.sum()
        logger.warning("Repairing %d invalid geometries", n_invalid)
        gdf.loc[invalid_mask, "geometry"] = gdf.loc[invalid_mask, "geometry"].apply(make_valid)

    null_mask = gdf.geometry.is_empty | gdf.geometry.isna()
    if null_mask.any():
        quarantine = gdf[null_mask].copy()
        quarantine["_quarantine_reason"] = "null_or_empty_geometry"
        quarantine.to_file(quarantine_path, driver="GeoJSON")
        logger.warning("Quarantined %d null/empty geometries → %s", null_mask.sum(), quarantine_path)
        gdf = gdf[~null_mask].copy()

    return gdf


def write_gold_geoparquet(
    gdf: gpd.GeoDataFrame,
    output_path: Path,
) -> None:
    """
    Write a Hilbert-sorted GeoParquet file with ZSTD compression.
    Sorting by Hilbert index co-locates spatially adjacent rows in the
    same row groups, maximising predicate pushdown effectiveness.
    """
    # Sort by Hilbert index for spatial locality
    centroids = gdf.geometry.centroid
    gdf = gdf.copy()
    gdf["_hilbert_x"] = centroids.x
    gdf["_hilbert_y"] = centroids.y
    gdf = gdf.sort_values(["_hilbert_y", "_hilbert_x"]).drop(
        columns=["_hilbert_x", "_hilbert_y"]
    )

    gdf.to_parquet(
        output_path,
        compression=COMPRESSION,
        compression_level=COMPRESSION_LEVEL,
        row_group_size=ROW_GROUP_SIZE,
        write_covering_bbox=True,   # embeds bbox statistics for predicate pushdown
        schema_version="1.1.0",
    )
    logger.info("Gold GeoParquet written → %s (%d features)", output_path, len(gdf))


def run_pipeline(
    source: Path,
    gold_dir: Path,
    quarantine_dir: Path,
) -> None:
    gold_dir.mkdir(parents=True, exist_ok=True)
    quarantine_dir.mkdir(parents=True, exist_ok=True)

    silver_gdf = silver_transform(
        source_path=source,
        quarantine_path=quarantine_dir / f"{source.stem}_quarantine.geojson",
    )
    write_gold_geoparquet(
        gdf=silver_gdf,
        output_path=gold_dir / f"{source.stem}.parquet",
    )

Key design decisions in this implementation:

• Explicit CRS check before any spatial operation: reprojecting after indexing produces incorrect results.
• Quarantine partition instead of silent drops: invalid geometries are traceable, not lost.
• Hilbert-order sort before writing: spatially adjacent features land in the same row groups, improving bounding-box predicate pushdown at query time.
• write_covering_bbox=True: embeds per-row-group bounding box statistics that DuckDB, Trino, and Athena use for partition pruning.

#6. Schema Governance and Compliance

#Schema Evolution

Geospatial schemas drift as new attributes are collected or coordinate precision requirements change. GeoParquet handles additive schema changes gracefully — adding new columns is non-breaking for downstream readers that use column projection. However, dropping or renaming geometry columns breaks all downstream consumers. Register schemas in a catalog (Apache Atlas, AWS Glue Data Catalog, or a custom metadata store) and enforce backward-compatibility rules through CI/CD validation hooks at the pipeline boundary.

Track the following per dataset: EPSG code, geometry type (point/polygon/multipolygon), coordinate precision (decimal places), encoding standard (WKB ISO vs Extended WKB), and GeoParquet spec version. Schema drift without version tracking is the primary cause of silent corruption in long-running pipelines.

#CRS Governance in Distributed Indexes

Spatial indexing assumes a consistent coordinate space across all records. When ingesting multi-CRS data into a single indexed file, the index must normalise coordinates first or store projection metadata per feature. Mixing projections silently corrupts bounding-box metadata and produces incorrect spatial joins — a failure mode that only surfaces during cross-dataset operations, often months after ingest.

Enforce a single canonical CRS at the silver layer boundary and document deviations explicitly in the metadata catalog.

#Access Control and Compliance

Spatial data frequently contains sensitive location information: infrastructure positions, asset tracking, population health data. Cloud object storage provides bucket-level IAM policies, but fine-grained feature-level masking requires application-layer enforcement. For regulated workloads, implement column-level encryption for geometry columns before writing to the gold layer, and maintain audit logs of all spatial query patterns.

Government and enterprise contracts frequently mandate specific retention periods, encryption standards (AES-256 at rest, TLS 1.3 in transit), and provenance tracking. Embed pipeline run IDs and source data lineage as GeoParquet file metadata so every gold-layer file carries a complete audit trail.

#7. Strategic Selection Framework

Choosing a storage format should be driven by workload characteristics rather than ecosystem familiarity. Use this decision matrix as a starting point, then validate against your measured latency and throughput requirements.

When building hybrid systems, maintain a single source of truth in GeoParquet for analytics while generating FlatGeobuf derivatives on demand for web services. This polyglot approach minimises storage duplication while optimising access paths for distinct consumer profiles. At scale, FlatGeobuf generation can be triggered as an async job when the gold-layer partition is written, with outputs cached in a CDN-backed object bucket for tile server consumption.

#Frequently Asked Questions

When should I keep Shapefiles in my pipeline?

Shapefiles should only persist as ingestion artifacts in the bronze layer. Convert them to GeoParquet or FlatGeobuf immediately after landing. The multi-file structure, 2 GB size cap, and 10-character field-name limit make Shapefiles unsuitable as a persistent storage format in any modern cloud or distributed system. The conversion cost is trivially low compared to the operational risk of building workflows on top of the format’s limitations.

GeoParquet or FlatGeobuf for tile serving?

FlatGeobuf is the better choice for tile and feature APIs that serve dynamic bounding-box queries. Its packed Hilbert R-tree header enables HTTP Range requests so clients fetch only the relevant spatial subset without downloading the full file. GeoParquet is optimal for server-side analytics that scan and aggregate large row groups — use it as the source of truth and generate FlatGeobuf derivatives on demand from the gold layer.

What compression codec should I use for GeoParquet geometry columns?

ZSTD at level 3–6 delivers the best analytical pipeline trade-off: high compression ratios on coordinate-heavy geometry columns with acceptable decompression throughput. For interactive APIs where latency dominates, Snappy or LZ4 reduce CPU overhead at the cost of slightly larger files. Avoid GZIP for geometry columns — decompression is CPU-expensive and parallelises poorly in columnar scan workloads.

How do I handle mixed CRS sources in a single GeoParquet file?

Normalise all geometries to a single CRS — typically EPSG:4326 for interchange or a projected CRS like EPSG:3857 for metric calculations — before writing the file. Mixing projections silently corrupts spatial indexes and bounding-box filters. Embed the target EPSG code in the GeoParquet metadata column and add a CRS validation step at the silver layer of your pipeline. Never defer this normalisation to query time.

Can DuckDB query GeoParquet files directly in object storage?

Yes. DuckDB’s spatial extension reads GeoParquet files from S3, GCS, or Azure Blob Storage via its httpfs extension. It applies predicate pushdown and row-group pruning so only relevant column chunks are fetched over the network. This makes DuckDB an excellent serverless query engine for ad-hoc spatial analysis without a persistent cluster.

What is the right row-group size for spatially partitioned GeoParquet?

Target 100,000–500,000 rows per row group for most analytical workloads. Smaller row groups improve predicate pushdown selectivity for point queries but increase per-file metadata overhead. For polygon-heavy datasets with large geometry columns, lean toward 100,000 rows to keep decompressed memory pressure manageable. Always benchmark against your actual query engine’s read-ahead and prefetch settings — the optimal size varies between DuckDB, Trino, and Spark.

#Related

• Understanding Parquet Columnar Storage for GIS
• Comparing GeoParquet vs FlatGeobuf Performance
• Shapefile Limitations in Modern Data Stacks
• GeoJSON Overhead and Serialization Costs
• ZSTD Compression Levels for Geospatial Data
• Dictionary Encoding for Categorical GIS Attributes
• Spatial Partitioning with Quadtree Indexes
• Data Conversion & Migration Pipelines
• Compression, Chunking & Spatial Indexing

← Back to Home