#ZSTD Compression Levels for Geospatial Data

Geospatial datasets present a distinct compression problem: coordinate arrays carry high spatial autocorrelation, categorical attribute columns repeat values across thousands of rows, and raster bands vary sharply in entropy from one tile to the next. Picking the wrong Zstandard level costs real money — either in storage bills from over-trusting slow but small codecs, or in query latency from decompression becoming the bottleneck. This guide gives GIS data engineers a structured workflow to profile, configure, and validate ZSTD across production vector and raster workloads. It sits within the broader Compression, Chunking & Spatial Indexing topic area.

#Prerequisites

• Python 3.10+ with pyarrow>=14.0, zstandard>=0.22.0, geopandas>=0.14, and rasterio>=1.3
• Representative samples: multi-polygon vector datasets (administrative boundaries, parcel layers) and multi-band raster tiles (Sentinel-2 bands, DEMs)
• Columnar storage familiarity: Parquet row groups, column statistics, and predicate pushdown mechanics
• Monitoring tools: psutil or tracemalloc for memory profiling; time.perf_counter for throughput measurement
• Cloud storage access: an S3-compatible endpoint for validating chunked reads and round-trip decompression latency at network speeds

#Architectural Foundations

Zstandard is a streaming compressor built on a finite-state entropy (FSE/ANS) back-end with a configurable LZ77-style sliding window. Its level parameter trades CPU and memory budget for dictionary search depth: higher levels spend more cycles finding longer back-references, which matters most when the input has long-range repetition.

Geospatial data feeds this well because:

• Coordinate arrays delta-encode naturally. After differencing, consecutive vertices in a polygon boundary differ by small integers, producing runs that a deep ZSTD search can reference across hundreds of kilobytes.
• Categorical attributes — land-cover codes, administrative IDs, status flags — repeat densely. Combining dictionary encoding for categorical GIS attributes at the Parquet column level before ZSTD cuts the symbol alphabet to 16–256 entries, letting even level 3 achieve near-maximum ratio on those columns.
• Raster band data is the exception: adjacent pixels in natural imagery carry high spatial frequency detail, and entropy is already close to the theoretical limit. ZSTD levels above 6 rarely improve raster ratio by more than 2–3% while adding 4–6x more compression time.

The interplay between window size and chunk boundaries matters here. ZSTD operates per-column-chunk inside Parquet, and because Parquet stores GIS data column by column, each geometry, attribute, and category column is compressed as a contiguous byte stream — exactly the long, homogeneous input ZSTD’s deeper levels exploit. ZSTD’s windowLog parameter controls how far back the compressor can search. If your row group size for Parquet is smaller than the effective window, you are leaving ratio on the table; if it is much larger, decompression must load the full group into memory before any query filtering can occur.

The diagram below shows how the level-selection decision flows through a typical ingest pipeline:

#Step-by-Step Workflow

#1. Profile Data Entropy and Geometry Distribution

Geospatial data rarely exhibits uniform entropy across its columns. Before locking in any ZSTD level, scan your dataset for column-level statistics using pyarrow.parquet.read_metadata — row group min/max ranges reveal which columns have high cardinality versus dense repetition.

# pyarrow>=14.0
import pyarrow.parquet as pq

meta = pq.read_metadata("parcels.parquet")
for rg_idx in range(meta.num_row_groups):
    rg = meta.row_group(rg_idx)
    for col_idx in range(rg.num_columns):
        col = rg.column(col_idx)
        stats = col.statistics
        if stats and stats.has_min_max:
            span = stats.max - stats.min if isinstance(stats.max, (int, float)) else None
            print(f"  col={col.path_in_schema}  null_count={stats.null_count}  span={span}")

Pay attention to geometry columns: after WKB encoding, the byte range of coordinate values reflects spatial extent, not repetition. Coordinate arrays with tight bounding boxes (dense urban parcels) compress far better at higher levels than sparse, continent-scale geometries where adjacent vertices may be kilometres apart.

Align ZSTD’s windowLog parameter with your typical geometry extent. The default windowLog=27 (128 MB window) handles most municipal-scale polygon layers. For large continental boundaries, consider windowLog=28 or windowLog=29 when using the zstandard Python library directly.

#2. Map Workload to ZSTD Level

Zstandard supports levels 1–22 (plus 0 as an alias for the default). Not all ranges are practical for geospatial pipelines:

For a detailed breakdown by dataset type, see the optimal ZSTD levels for vector vs raster data reference page.

Raster workloads (Sentinel-2, DEMs, land-cover grids) rarely benefit from going above level 4–6 because natural imagery is already near the entropy ceiling — you pay 3–4x more CPU for under 2% extra compression. Vector workloads, especially polygon boundary layers with long coordinate runs, respond strongly to levels up to 12–14 because the deeper search window finds repeated sub-sequences across entire ring buffers.

#3. Configure Chunk Boundaries and Row Group Size

Compression efficiency scales with chunk size. Smaller chunks limit ZSTD’s search window and shrink ratio potential; larger chunks improve ratio but increase peak memory during decompression and can stall concurrent query workers.

The interaction is concrete: at level 9, a 50 MB row group typically achieves 60–65% of the ratio you would get at 500 MB, because the compressor exhausts its back-reference window midway through the smaller group. However, a 500 MB row group means your query engine must decompress the entire group before any predicate filtering takes effect.

For most GeoParquet workloads, 128–256 MB row groups hit the best tradeoff: the ZSTD back-reference window is well-populated, and decompression memory stays manageable. Coordinate this with your overall row group sizing strategy rather than treating chunk size and compression level as independent knobs.

When co-locating data with spatial partitions — for instance, files partitioned by Hilbert curve range — smaller row groups (32–64 MB) may be warranted because each file already covers a compact spatial extent, and predicate pushdown at the file level compensates for the reduced per-group ratio.

#4. Benchmark Compression Ratio vs Decompression Throughput

Run controlled benchmarks on your actual data before committing to a level. See the production-ready implementation section below for a complete, reusable benchmarking function. When testing categorical attribute columns (land-cover codes, administrative classification codes), also measure whether applying dictionary encoding for categorical GIS attributes before ZSTD reduces both file size and decompression time — dictionary-encoded columns compress to near-theoretical limits even at level 3.

#5. Validate Spatial Query Performance

Compression tuning is only meaningful if spatial operations remain fast. After writing your benchmarked files, run representative queries and measure end-to-end latency:

• Point-in-polygon: latency against administrative boundary layers at level 3 vs level 9
• Raster band extraction: tile read time for multi-band GeoTIFF encoded at different ZSTD levels
• Spatial join: memory-peak tracking during large geometry intersection scans

Watch predicate pushdown behaviour closely. ZSTD compresses entire row groups, so a query engine must decompress the full group before any spatial predicate can filter rows. If your workload is highly selective (returning less than 1% of rows per group), pair lower ZSTD levels (3–6) with tighter spatial partitioning with quadtree indexes to minimise the volume of data decompressed per query.

# Quick DuckDB query-latency check — duckdb>=0.10.0
import time
import duckdb

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

for level in [3, 6, 9, 12]:
    path = f"parcels_level_{level}.parquet"
    t0 = time.perf_counter()
    result = conn.execute(
        f"""
        SELECT COUNT(*) FROM read_parquet('{path}')
        WHERE ST_Within(
            ST_GeomFromWKB(geometry),
            ST_GeomFromText('POLYGON((...))')
        )
        """
    ).fetchone()
    elapsed = time.perf_counter() - t0
    print(f"Level {level}: {result[0]} rows in {elapsed:.3f}s")

#Production-Ready Implementation

The function below benchmarks multiple ZSTD levels against a pyarrow.Table, tracks both write and read throughput, and cleans up scratch files on every exit path. Drop it into your ingest pipeline to establish a level baseline before deploying to production.

# pyarrow>=14.0, zstandard>=0.22.0
from __future__ import annotations

import os
import time
import tracemalloc
from pathlib import Path
from typing import Any

import pyarrow as pa
import pyarrow.parquet as pq


def benchmark_zstd_levels(
    table: pa.Table,
    levels: list[int] | None = None,
    row_group_size: int = 250_000,
    scratch_dir: Path = Path("/tmp/zstd_bench"),
) -> list[dict[str, Any]]:
    """Benchmark ZSTD compression at multiple levels on a PyArrow Table.

    Args:
        table: The PyArrow Table to benchmark.
        levels: ZSTD levels to test. Defaults to [3, 6, 9, 12].
        row_group_size: Row group size passed to write_table.
        scratch_dir: Temporary directory for benchmark files.

    Returns:
        List of dicts with write/read timings, file size, and peak memory.
    """
    if levels is None:
        levels = [3, 6, 9, 12]

    scratch_dir.mkdir(parents=True, exist_ok=True)
    results: list[dict[str, Any]] = []

    for level in levels:
        path = scratch_dir / f"bench_level_{level}.parquet"
        try:
            # --- Write ---
            tracemalloc.start()
            t0 = time.perf_counter()
            pq.write_table(
                table,
                str(path),
                compression="zstd",
                compression_level=level,
                row_group_size=row_group_size,
                write_statistics=True,
            )
            write_s = time.perf_counter() - t0
            _, peak_write = tracemalloc.get_traced_memory()
            tracemalloc.stop()

            file_mb = path.stat().st_size / 1024**2

            # --- Read ---
            tracemalloc.start()
            t0 = time.perf_counter()
            _ = pq.read_table(str(path))
            read_s = time.perf_counter() - t0
            _, peak_read = tracemalloc.get_traced_memory()
            tracemalloc.stop()

            results.append({
                "level": level,
                "write_s": round(write_s, 3),
                "read_s": round(read_s, 3),
                "file_mb": round(file_mb, 2),
                "ratio": round(
                    (table.nbytes / 1024**2) / file_mb, 2
                ) if file_mb > 0 else None,
                "peak_write_mb": round(peak_write / 1024**2, 2),
                "peak_read_mb": round(peak_read / 1024**2, 2),
            })

        except Exception as exc:
            results.append({"level": level, "error": str(exc)})
        finally:
            if path.exists():
                path.unlink()

    return results


# Example call — replace with your own GeoParquet table:
if __name__ == "__main__":
    import geopandas as gpd  # geopandas>=0.14

    gdf = gpd.read_file("admin_boundaries.gpkg")
    tbl = gdf.to_arrow()  # requires geopandas>=0.14 + pyarrow>=14.0
    for row in benchmark_zstd_levels(tbl, levels=[3, 6, 9, 12, 15]):
        print(row)

#Reference Matrix: Level vs Measured Outcomes

The values below are representative ranges from polygon boundary layers (50 M vertices, GeoParquet) on a 4-core cloud VM. Your numbers will vary by geometry complexity and hardware.

#Advanced Configuration for Production Deployment

#Cloud Storage and Cold-Tier Optimisation

When migrating historical datasets to object storage, network egress and retrieval costs often dominate over compute. For infrequently accessed layers, levels 12–15 can reduce storage bills by 30–40% compared to default Snappy, with decompression still measured in seconds for typical 256 MB row groups. However, cold-storage retrieval latency compounds with decompression time. Implement tiered compression policies aligned with your S3 lifecycle rules: levels 3–6 for standard-tier, 12–15 for infrequent-access, and only evaluate levels above 15 for Glacier-style deep archive where restore time is already measured in hours.

#Production Checklist

1. Validate codec compatibility: modern engines (DuckDB 0.9+, Trino 420+, Athena, Spark 3.3+) all support ZSTD in Parquet metadata. Only legacy tools (Hive 1.x) need a fallback to Snappy.
2. Cap decompression memory: configure your runtime and PyArrow memory pool to prevent OOM during concurrent spatial scans. A rule of thumb: allow 3–4× the uncompressed row group size per concurrent thread.
3. Implement fallback codec routing: route failed decompression attempts to a secondary codec pool (LZ4 or uncompressed) to maintain pipeline resilience during index rebuilds or schema migrations.
4. Monitor compression drift: re-benchmark quarterly as dataset characteristics evolve (new geometry types, updated raster resolutions, schema additions). A single new high-cardinality column can shift the optimal level by 2–3 tiers.

#Failure Modes and Gotchas

#Frequently Asked Questions

Which ZSTD level should I use for interactive tile serving?

Levels 3–6 offer the best balance for interactive workloads. They decompress at 400–600 MB/s on a single core, keeping tile-server latency well under 50 ms, while still achieving 30–45% better ratio than Snappy on coordinate arrays.

Does ZSTD level affect GeoParquet compatibility with DuckDB or Athena?

No. The Parquet codec metadata records ZSTD regardless of level, and all modern engines (DuckDB 0.9+, Trino, AWS Athena, Spark 3.3+) support it. Only very old engines (Hive 1.x, legacy Impala builds) need a fallback to Snappy or uncompressed.

Should I apply dictionary encoding before or after ZSTD?

Apply dictionary encoding first, at the Parquet column level. ZSTD then compresses the already-reduced representation, typically yielding an additional 10–20% reduction on top of dictionary encoding alone. Reversing the order wastes CPU without improving ratio.

What ZSTD level is best for cold-tier archival of large raster datasets?

Levels 12–15 reduce storage by 30–40% compared to level 3, with decompression still completing in seconds for typical row groups. Avoid levels 16–22 in production archival pipelines: the marginal ratio gain (1–3%) rarely justifies the 4–10x memory overhead during restore.

#Related

• Optimal ZSTD Levels for Vector vs Raster Data
• Row Group Sizing Strategies for Parquet
• Dictionary Encoding for Categorical GIS Attributes
• Spatial Partitioning with Quadtree Indexes
• Tuning Row Group Size for Cloud Query Performance

← Back to Compression, Chunking & Spatial Indexing