gridpack

Lossless compression for gridded scientific data (ERA5 and friends).

gridpack shrinks NetCDF/HDF5 gridded fields — (time × lat × lon) float arrays — without changing a single bit. decode(encode(x)) == x holds byte-for-byte, including -0.0, NaN and ±Inf. It also exposes an array-level codec so other projects can use it as a drop-in compression layer.

On real ERA5 10 m wind fields it reaches 10.4× vs the raw float bytes while staying bit-exact, and stores every block with a SHA-256 so silent corruption is impossible.

Why it works

Most gridded products are not arbitrary floats. ERA5 fields, for example, sit on an exact binary lattice: 10 m wind is quantised to 2**-4 (0.0625 m/s), total precipitation to 2**-14, 2 m temperature to 2**-4. gridpack detects that lattice exactly and rewrites the field as integers k = v / 2**E. Small integers compress far better than 32-bit float bit patterns.

Detection is not a heuristic. For a non-zero float32 value v = sig · 2**(exp-150) with a 24-bit significand sig, v is a multiple of 2**E iff tz(sig) + exp - 150 >= E, so the largest usable exponent is E_max = min(tz(sig) + exp - 150). The code verifies the rewrite is exact before committing to it, and falls back to an order-preserving bit map otherwise.

-0.0 cannot be represented as an integer, so its positions are stored separately as a sparse bitmask — that is what keeps the lattice path bit-exact.

Install

pip install gridpack                 # core: numpy + h5py (zstandard optional, recommended)
pip install "gridpack[xarray]"       # + xarray backend support

Or from source:

git clone https://github.com/showkeyjar/gridpack && cd gridpack
pip install -e .

Installing the package also registers the xarray backend entry point, so engine="gpz" works immediately -- see xarray.

Command line

python -m gridpack compress   in.nc  out.gpz [options]
python -m gridpack decompress in.gpz out.nc
python -m gridpack verify     in.nc  out.gpz      # byte-for-byte comparison
python -m gridpack info       out.gpz
python -m gridpack bench      in.nc --nt 48

Example:

python -m gridpack compress era5_10u.nc era5_10u.gpz --level 12 --entropy zstd --mode auto
python -m gridpack verify   era5_10u.nc era5_10u.gpz
# EXACT

Parameters

compress

--mode auto is per-variable: the codec probes the first block with every candidate and locks in the winner for that dataset. Different variables routinely pick different models.

Python API

import gridpack

# file level
stats = gridpack.compress_file(src, dst, level=12, chunk_t=16, threads=4,
                               quant_bits=0, mode="auto", entropy="zstd",
                               file_chunk=None, verify=True)
gridpack.verify_file(src, dst)          # -> {"exact": True, "details": [...]}
gridpack.decompress_file(dst, out_nc)

# array level
blob = gridpack.encode_array(arr, level=12, chunk_t=16, threads=1,
                             quant_bits=0, mode="dt", entropy="zstd")
arr2 = gridpack.decode_array(blob, arr.shape, arr.dtype)

# what representation will be used?
gridpack.detect_lattice_e(arr)          # -> -4, or None if no lattice

chunk_t=None encodes the array as a single block. With chunk_t=N the array is split along axis 0 and each chunk is stored independently, so decode_array can rebuild without the whole payload in memory.

xarray

Installing the package registers an xarray backend entry point, so a .gpz file opens like any other dataset:

import xarray as xr
import gridpack                                  # registers backend + accessor

ds = xr.open_dataset("era5_10u.gpz", engine="gpz")
ds = xr.open_dataset("era5_10u.gpz", engine="gpz",
                     chunks={"time0": 24})       # dask-backed, lazy chunks

engine="gpz" is also selected automatically by xr.open_dataset("...gpz"), because the backend's guess_can_open matches the .gpz suffix. Dimensions and coordinates are rebuilt from the dimension-scale links stored in the container, and every value is bit-identical to the source .nc.

xarray has no generic write backend interface, so writing goes through the accessor:

ds.gpz.to_file("out.gpz", level=12, threads=4)    # -> compress_file() stats dict
ds.gpz.info("out.gpz")                            # -> ratio, per-variable model, ...

gridpack.from_xarray(ds, "out.gpz")               # same, without the accessor
gridpack.open_xarray("out.gpz")                   # same as xr.open_dataset(..., engine="gpz")

to_file writes a temporary netCDF (needs h5netcdf or netCDF4) and compresses that, then removes it. netCDF-internal attributes (_NCProperties, CLASS, _Netcdf4Coordinates, ...) are stripped automatically -- if a backend rejects any further reserved name, it is dropped and the write is retried.

Verified on the bundled ERA5 subset (120 time steps, 4 variables): all four variables round-trip bit-exactly through engine="gpz", and a ds.isel(time0=slice(0, 6)).gpz.to_file() write-back re-reads bit-exactly at 10.30x.

python tests/test_xarray.py

Benchmarks

Real ERA5, 10 m eastward wind (eastward_wind_at_10_metres, float32, 120 time steps, 47.53 MB raw). Full reproducible output in docs/BENCHMARKS.md; regenerate with python scripts/bench_era5.py <file.nc>.

Representation × time model (entropy = zstd level 12)

Two things worth noting. The lattice representation roughly doubles the ratio (10.36× vs 5.25×). And with the plain order map, adding prediction makes things worse (5.25× → 4.29×): the order-preserving bit map has a large step at the zero crossing, so a field that changes sign produces enormous residuals. Integer-lattice residuals do not have that problem.

Entropy backend (first 200 000 bytes of the same byte-plane buffer)

The adaptive context range coder (rc1, context = previous byte × bit index) beats zstd by 12.6 % on this data, and rc0 by 4.7 %. That is a real gain but it costs ~15× the time in pure Python, so zstd stays the production default and rc* are there for offline re-encoding or archival runs where CPU time is cheap.

Whole file

Ratio depends heavily on how the source was written. Fields already quantised (most modern ERA5 monthly files) land at ~10×; unquantised, full-precision months land at ~4×; a .nc that is itself zlib-compressed by netCDF will show a smaller gain versus that .nc, so compare against the raw float bytes to get a like-for-like number.

Throughput: 0.4.0 -> 0.5.0

Lattice detection dominated compression time, and inside it the helper that counts trailing zeros. 0.5.0 replaces that helper's shift-and-mask loop with a two-level 16-bit lookup table.

Measured on an ERA5 2 m temperature file (483.1 MB on disk, 3 089.8 MB of raw float32, 6 x 24 time steps; level=3, chunk_t=16, threads=4):

detect_lattice fell from 13 % of the wall clock to 8.6 %.

The compressed output is unchanged: same 296.1 MB container, same 1.631x ratio, and verify_file still reports byte-exact. This is a pure speed-up — no format change, no API change, CODEC_VERSION untouched, and .gpz files written by 0.4.0 still read.

The lookup table was validated against the original implementation over every integer from 1 to 220 and over 2 x 105 random 32-bit values (identical results), and the detected lattice exponent on real ERA5 data is unchanged (E = -4).

Correctness

tests/test_roundtrip.py covers every path and asserts byte equality, not approximate equality:

python tests/test_roundtrip.py                 # unit + array level
GRIDPACK_TEST_NC=file.nc python tests/test_roundtrip.py   # + NetCDF end-to-end

It checks the order map on ±0.0 / ±Inf / NaN / subnormals, lattice detection and round-trip for float32 and float64, negative-zero preservation, all five entropy backends (including single-symbol and empty inputs), all time models on random / lattice / constant / NaN-bearing / integer / 1-D arrays, chunked encoding, and a full NetCDF compress → verify → decompress cycle.

Container guarantees:

• every block carries a SHA-256 of its raw bytes; a mismatch raises instead of returning wrong data,
• the restored .nc keeps variable names, shapes, dtypes, dimension scales, chunking, compression settings, _FillValue and all attributes,
• peak memory is bounded by the file chunk size, not by the file size.

Container format

.gpz (MAGIC = b"GPZ1"): JSON metadata header + a payload index + the compressed blocks. Metadata records the source name/size, codec version, and per-dataset shape, dtype, chunking, attributes, dimension links, chosen mode, and the detected lattice exponent.

Each block starts with a tag byte (0 raw, 1 predicted, 2 byte-planes only, 3 lattice + byte-planes, 4 lattice + predicted), an entropy-backend id, and — for the lattice tags — the exponent and the length of the negative-zero bitmap.

Limitations

• The lattice path only engages when every value shares one exponent and the resulting integers fit comfortably in int32; otherwise the codec falls back to the order map.
• .gpz needs this library to read. It is not a NetCDF file. If you need .nc compatibility, decompress back.
• rc0/rc1/huff are pure Python and unsuitable for large files.
• Throughput is tens of MB/s of raw float data per process; see Throughput: 0.4.0 -> 0.5.0 above for measured numbers. Parallelise across files.
• Only quant_bits=0 is lossless; anything else deliberately discards information.

中文说明

gridpack 是面向 ERA5 这类格点科学数据的严格无损压缩与读写工具包：解码后与原始 数据逐字节相同（含 -0.0、NaN、±Inf），实测在真实 ERA5 10m 风场上达到 10.37 倍压缩（相对原始 float 字节），并支持按需读取。

核心思路：ERA5 各要素的值并不是任意浮点数，而是精确落在同一个二进制格点 2**E 上（10m 风量化到 2**-4，总降水 2**-14，2m 气温 2**-4）。工具会精确判定这个 格点并把浮点无损改写成整数 k = v / 2**E，再做时间/空间一阶差分、zig-zag、字节平面重排 与熵编码。负零单独用稀疏位图记录，保证位级无损。

关键参数：--mode（none/spatial/dt/lagN/auto，按变量自动选）、 --entropy（zstd/zlib/rc0/rc1/huff）、--level（zstd 1–22，默认 12）、 --chunk-t（内部时间分块，默认 16）、--file-chunk（文件级分块，0 = 自适应约 96 MB）、 --threads、--quant-bits（0 = 严格无损，>0 为有容差模式）。

实测（真实 ERA5 10m 风场，120 时次，47.53 MB 原始）：格点整数 + dt 达 10.36×， 而保序映射最好只有 5.25×——格点表示带来近一倍收益；熵编码上 rc1 比 zstd 再小 12.6%，但纯 Python 慢约 15 倍，故生产默认仍用 zstd。整文件 475.28 MB → 45.83 MB， 10.37×，耗时 47.6 s，逐字节校验通过。

xarray 支持：pip install gridpack 后 xarray 后端即注册完成，.gpz 可直接 xr.open_dataset("x.gpz", engine="gpz") 打开（.gpz 后缀也会被自动识别）；维度与坐标 从容器里记录的维标关联还原，数值与原始 .nc 逐位一致。写回用 ds.gpz.to_file("out.gpz", level=12, threads=4)，也可用 gridpack.from_xarray(ds, path)； 查看某个 .gpz 的压缩率与逐变量选中的模型用 ds.gpz.info(path)。实测 120 时次 ERA5 子集 的 4 个变量经 engine="gpz" 读出全部逐位通过，写回再读回同样逐位通过（10.30×）。

License

MIT — see LICENSE.