Optimal ZSTD Levels for Vector vs Raster Data

For production geospatial workloads, the optimal ZSTD level for vector data is 3–6 and for raster data is 1–3. Vector formats benefit from higher levels because coordinate sequences, attribute dictionaries, and CRS identifiers repeat heavily across features — ZSTD’s sliding-window entropy model captures this structure efficiently. Raster formats require lower levels because they depend on chunk- and tile-based random access over HTTP range requests; decompression time scales nearly linearly with level, so pushing above level 3 inflates tile latency without meaningfully shrinking tile size. Levels above 6 for vector yield under 2% additional ratio while increasing CPU time by 30–50% and peak memory during encoding. This page works through the reasoning, gives you ready-to-run Python, and maps the edge cases where these defaults do not apply. For the full algorithmic context and format-level reference matrix, see ZSTD Compression Levels for Geospatial Data.

Quick-Reference Table

Data Type Recommended Level Typical Ratio Gain vs Uncompressed Decode Latency Primary Use Case
Vector — GeoParquet, FlatGeobuf, GeoPackage (interactive) 3 60–75% < 2 ms per row group Feature servers, WFS, map tile generation
Vector — GeoParquet, FlatGeobuf (archival / batch ETL) 5–6 75–85% 3–6 ms per row group Data lakes, ML training, cold analytical scans
Raster — Cloud-Optimized GeoTIFF, Zarr (interactive) 1–2 40–55% < 3 ms per tile Web mapping, tile servers, real-time dashboards
Raster — Zarr, NetCDF4 (batch analysis) 3 50–60% < 5 ms per tile Batch post-processing, scientific workflows

Why These Levels Work

ZSTD Level vs Compression Ratio and Decode Latency Chart showing that vector data ratio gains flatten above level 6, while raster decode latency climbs steeply from level 3 onward. Two shaded bands mark the recommended ranges: levels 1 to 3 for raster, levels 3 to 6 for vector. Relative scale (%) Raster: 1–3 Vector: 3–6 1 2 3 4 5 6 7 8 9 ZSTD compression level 0 25 50 75 100 Vector ratio Raster latency Shaded bands mark the recommended level ranges per data type. Curves are illustrative; actual values vary with dataset entropy.

Vector geometries and attribute data exhibit high structural locality. Coordinate deltas between adjacent features are small integers; categorical attribute columns (land-use codes, administrative names, sensor IDs) repeat thousands of times within a single row group. ZSTD’s sliding window and dictionary mechanisms are well-suited to this structure: at level 3 you capture roughly 85% of the theoretical compression ceiling with near-zero latency overhead, and pushing to level 5–6 adds approximately 5–8% ratio for archival or batch ETL pipelines where encode time is amortised across scheduled runs. Modern columnar formats like GeoParquet default to ZSTD level 3 for precisely this reason — it balances chunk compression with fast predicate pushdown during analytical scans. When you also apply dictionary encoding for categorical GIS attributes, the two techniques compound: the dictionary removes redundancy within the column before ZSTD sees the data, letting ZSTD concentrate its window on cross-column patterns.

Raster datasets are structurally different. Cloud-Optimized GeoTIFFs and Zarr arrays are divided into spatial tiles that clients fetch individually over HTTP range requests. Each tile must be fully decompressed before any pixel value can be read. ZSTD’s compression ratio scales with level sub-linearly for typical 32-bit float or 8-bit integer band data — but decompression time scales nearly linearly. At level 4+, tile decompression routinely crosses 10–15 ms, introducing visible lag in web mapping clients and throttling concurrent tile request throughput against object storage. At level 1–3, tiles decompress in under 5 ms while still achieving 40–60% size reduction. Pairing these levels with matching tile sizes (256×256 for web maps, 1024×1024 for batch analysis) ensures I/O parallelism rather than CPU saturation drives your throughput. This aligns with the principles covered in Compression, Chunking & Spatial Indexing, where physical layout governs which compression trade-offs are viable.

Production Implementation

The snippet below implements format-aware ZSTD levels for both pipelines. It uses typed signatures, explicit level control, tile/row-group alignment, and context-managed I/O — the minimum required for production use. Comments reference GDAL and library versions that introduced ZSTD support.

python
# Requires: pyarrow>=14.0, rasterio>=1.3, GDAL>=2.3
from __future__ import annotations

import logging
from pathlib import Path

import pyarrow as pa
import pyarrow.parquet as pq
import rasterio
from rasterio.enums import Resampling

logger = logging.getLogger(__name__)


def write_vector_zstd(
    table: pa.Table,
    output_path: str | Path,
    *,
    level: int = 3,
    row_group_size: int = 100_000,
) -> None:
    """Write vector data with format-aware ZSTD compression.

    Args:
        table: Arrow table containing geometry (WKB) and attribute columns.
        output_path: Destination Parquet file path (local or object-storage URI).
        level: ZSTD compression level. Use 3 for interactive workloads,
               5–6 for archival or batch ETL. Avoid >6 unless storage
               egress cost explicitly justifies the encode overhead.
        row_group_size: Rows per row group. Keep between 50k–200k to
                        balance dictionary effectiveness and scan granularity.
    """
    if level > 6:
        logger.warning(
            "ZSTD level %d for vector data yields <2%% extra ratio vs level 6 "
            "but increases encode CPU by 30–50%%. Consider whether the trade-off "
            "is justified for this workload.", level
        )

    pq.write_table(
        table,
        str(output_path),
        compression="zstd",
        compression_level=level,
        use_dictionary=True,     # Critical for categorical attribute columns
        write_statistics=True,   # Required for predicate pushdown
        row_group_size=row_group_size,
    )
    logger.info("Wrote vector Parquet at ZSTD level %d → %s", level, output_path)


def write_raster_zstd(
    src_path: str | Path,
    output_path: str | Path,
    *,
    level: int = 2,
    tile_size: int = 256,
) -> None:
    """Write a Cloud-Optimized GeoTIFF with raster-appropriate ZSTD compression.

    Args:
        src_path: Source raster file readable by GDAL/rasterio.
        output_path: Destination COG path. Must end in .tif or .tiff.
        level: ZSTD compression level. Keep at 1–3 to stay below 5 ms
               tile decompression. Level 1 for real-time tile servers;
               level 3 for batch post-processing.
        tile_size: Tile width and height in pixels. 256 for web maps,
                   512–1024 for batch analysis pipelines.
    """
    if level > 3:
        logger.warning(
            "ZSTD level %d for raster data risks tile decompression latency "
            "> 10 ms. Use level 1–3 for interactive workloads.", level
        )

    with rasterio.open(src_path) as src:
        profile = src.profile.copy()
        profile.update(
            driver="GTiff",
            compress="zstd",
            zstd_level=level,         # GDAL ZSTD encoder level (GDAL >= 2.3)
            tiled=True,
            blockxsize=tile_size,
            blockysize=tile_size,
            BIGTIFF="IF_SAFER",       # Auto-upgrade to BigTIFF for files > 4 GB
            NUM_THREADS="ALL_CPUS",   # Parallel tile encode
        )

        with rasterio.open(output_path, "w", **profile) as dst:
            for band_idx in range(1, src.count + 1):
                dst.write(src.read(band_idx), band_idx)

            # Build overviews so tile servers can serve zoom-out requests
            # without reading full-resolution tiles
            dst.build_overviews(
                [2, 4, 8, 16, 32],
                Resampling.average,
            )
            dst.update_tags(ns="rio_overview", resampling="average")

    logger.info(
        "Wrote COG at ZSTD level %d, tile size %dx%d → %s",
        level, tile_size, tile_size, output_path,
    )

Key notes:

  • compression_level in PyArrow maps directly to ZSTD levels 1–22. Setting use_dictionary=True alongside ZSTD is critical for categorical dictionary-encoded GIS attributes — the two techniques compound.
  • zstd_level in Rasterio sets the GDAL ZSTD encoder. ZSTD support in GeoTIFF requires GDAL ≥ 2.3; always confirm with gdal-config --version before deployment.
  • ZSTD’s internal window size grows with level. For levels > 4, monitor resident set size (RSS) during encoding of large single-chunk raster arrays.

Validation and Verification

Before committing a level to your pipeline, measure the actual ratio and decode latency against your own dataset distribution. The snippet below uses a Python microbenchmark that logs both metrics side by side:

python
# Requires: zstandard>=0.22, pyarrow>=14.0
import io
import time

import pyarrow.parquet as pq
import zstandard as zstd


def benchmark_zstd_levels(
    parquet_path: str,
    levels: list[int] | None = None,
) -> list[dict]:
    """Benchmark ZSTD levels against a real Parquet file's row groups.

    Returns a list of dicts with keys: level, compressed_bytes,
    ratio_pct, decode_ms.
    """
    if levels is None:
        levels = [1, 2, 3, 4, 5, 6, 9]

    # Read the first row group as raw Arrow record batch bytes
    pf = pq.ParquetFile(parquet_path)
    batch = pf.read_row_group(0)
    sink = io.BytesIO()
    with pq.ParquetWriter(sink, batch.schema, compression="none") as w:
        w.write_table(batch)
    raw_bytes = sink.getvalue()
    raw_size = len(raw_bytes)

    results = []
    for level in levels:
        cctx = zstd.ZstdCompressor(level=level)
        compressed = cctx.compress(raw_bytes)
        comp_size = len(compressed)
        ratio_pct = round((1 - comp_size / raw_size) * 100, 1)

        # Measure decompression (20 iterations for stable average)
        dctx = zstd.ZstdDecompressor()
        t0 = time.perf_counter()
        for _ in range(20):
            dctx.decompress(compressed)
        decode_ms = round((time.perf_counter() - t0) / 20 * 1000, 2)

        results.append({
            "level": level,
            "compressed_bytes": comp_size,
            "ratio_pct": ratio_pct,
            "decode_ms": decode_ms,
        })
        print(
            f"Level {level:2d}: {ratio_pct:5.1f}% ratio, "
            f"{decode_ms:6.2f} ms decode, "
            f"{comp_size // 1024} KB"
        )

    return results

Expected output on a typical GeoParquet vector file (100k rows, mixed geometry + categoricals):

markup
Level  1: 61.3% ratio,   0.82 ms decode, 1847 KB
Level  2: 63.1% ratio,   0.89 ms decode, 1771 KB
Level  3: 66.8% ratio,   0.94 ms decode, 1649 KB
Level  4: 68.2% ratio,   1.12 ms decode, 1601 KB
Level  5: 70.4% ratio,   1.38 ms decode, 1543 KB
Level  6: 71.9% ratio,   1.67 ms decode, 1507 KB
Level  9: 73.1% ratio,   2.05 ms decode, 1462 KB

The ratio curve flattens sharply after level 6 (under 2 percentage points gained from level 6 to 9) while decode time climbs 23%. For raster tile data, run the same benchmark against a representative tile’s raw bytes — expect decode latency to be the dominant constraint rather than ratio. You can also verify row group sizing is set correctly at the same time, since both parameters interact: undersized row groups limit ZSTD’s sliding-window context, artificially inflating measured ratios at higher levels.

Edge Cases and Caveats

Large single-chunk raster arrays. Some scientific workflows write entire Zarr arrays as a single chunk (for example, a 10 GB temperature grid with no tile split). In this pattern, ZSTD level 3 may actually perform better than level 1 because the longer input gives ZSTD’s window more context for dictionary matches, partially reversing the usual raster logic. Always benchmark with your actual chunk configuration before assuming the default recommendation applies.

GPU-accelerated vector pipelines. If your stack runs RAPIDS/cuDF for vector feature engineering, ZSTD is entirely CPU-bound, meaning it cannot benefit from GPU memory bandwidth. In GPU-heavy workflows, consider switching to LZ4 or Snappy for the hot in-memory representation and reserving ZSTD level 3 only for the final write to object storage. The CPU-bound encode step will not saturate your GPU, but it will compete with CPU-bound data preprocessing.

Mixed streaming and batch consumers. When a single GeoParquet dataset must serve both real-time API requests and scheduled analytical batch jobs, level 3 is the safest single choice. It keeps streaming decode under 2 ms per row group while still delivering 65–70% compression — neither the streaming path nor the batch path is materially handicapped. Avoid the temptation to split into two separate files at different levels unless your storage and ETL budget explicitly justifies the duplication overhead. The spatial partitioning decisions you make upstream also matter here: finer spatial partitions reduce the number of row groups touched per query, making the per-row-group decode cost at any level less significant.

Frequently Asked Questions

Can I use the same ZSTD level for both vector and raster data?

You can, but you should not. Vector and raster data have fundamentally different access patterns. Raster tile servers suffer visible latency increases above level 3 because decompression time scales nearly linearly with level. Vector analytical queries tolerate higher levels because they read entire row groups sequentially. Applying a single level to both data types will either under-compress your vector store or over-compress your raster tiles.

Why does ZSTD level 9–12 only make sense for archival vector data?

Levels 9–12 offer diminishing ratio gains — typically under 3% over level 6 — while increasing encode time by 5–10× and peak memory by 2–4×. For interactive or streaming workloads the added ratio simply is not worth the CPU cost. The only justified use is cold-storage archival of vector datasets where the data is written once, read infrequently, and storage egress costs dominate the budget.

Does ZSTD level affect GeoParquet predicate pushdown performance?

Yes, indirectly. Higher ZSTD levels do not affect the Parquet metadata or statistics blocks — pushdown filtering happens before decompression. However, once the filter narrows to specific row groups, those groups must be decompressed for column reads. A higher level means more CPU time per row group. For queries that regularly touch many row groups (low selectivity), this overhead accumulates. For highly selective queries (under 5% of row groups touched), the difference is negligible.

What ZSTD level should I use in a Zarr store for climate model output?

Level 1–2 for interactive analysis or dashboard workloads; level 3 if you also serve the store for batch post-processing. Climate model output (float32 temperature, pressure, humidity grids) typically has moderate entropy after delta encoding, so level 3 rarely outperforms level 1 by more than 8–10% in ratio. The decode latency penalty from jumping to level 4+ almost always outweighs the marginal storage saving on chunked climate arrays.


← Back to ZSTD Compression Levels for Geospatial Data