ERA5 is already gzipped. I still made it 47% smaller — losslessly.

ERA5 is already gzipped. I still made it 47% smaller — losslessly. First, the two baselines The easiest mistake to make when you write a compressor is to quote a flattering baseline while talking about a different one. So both go on the table first: Baseline Result vs the .nc on disk (ERA5 ships gzip-compressed) 1.86–1.88×, 46–47% smaller vs uncompressed float32 bytes median 10.27× Same files, two rulers. The orange bar is what you actually get back on disk. Both numbers are true, but they describe different things. The .nc files are already gzip-compressed — so the disk you actually reclaim is 47%, not 90%. And 10.27× matters because it is the same ruler everyone else uses. fpzip, zfp's reversible mode, and zstd + byte shuffle typically land at 1.5–3× on uncompressed scientific fields. Same ruler, an order of magnitude apart. Why you can compress already-compressed data Because ERA5's floats are not arbitrary floats. Every ERA5 variable is quantized to a fixed precision when it is generated — 10 m wind to 0.0625 (= 2⁻⁴), total precipitation to 2⁻¹⁴. Which means every value in a file is an integer multiple of the same 2^E. gzip cannot see this. It looks at 32-bit patterns, and the long run of zeros in the low mantissa bits is indistinguishable from random bytes to it: LZ77 finds no repeated pattern, and Huffman cannot give zeros a shorter code, because as far as gzip is concerned the high bytes are the content. To us, that run of zeros is where all the information is. Finding E Detection is purely algebraic, with no approximation: float32 non-zero value = sig × 2^(exp−150) # sig is the 24-bit significand v is a multiple of 2^E ⟺ tz(sig) + exp − 150 ≥ E largest E valid for every value = min(tz(sig) + exp − 150) (tz is the count of trailing zero bits. One trap: exp must have the sign bit masked off first, or negative values get their exponent over-counted by 256 — we hit this bug, and it silently turns lossless compression into lossy compression.) Given E, the float is rewritten losslessly as an integer k = v / 2^E. Because it is only a power-of-two scaling, this step has no rounding error. Where the saving comes from A concrete example. Take a grid point whose value is 273.1875 in a file with E = −4: raw float32: 273.1875 → k = 273.1875 / 0.0625 = 4371 right neighbour k = 4373 longitude diff → 2 zigzag → 4 = 0x00000004 byte planes → top 3 bytes all zero, only the lowest byte is 4 Across a whole run of values like that, the high byte planes are almost entirely zero and get eaten whole by the compressor; only the low planes actually cost bits. The real bit pattern of 248.375: 13 trailing zeros in the mantissa are the gap quantization left behind. gzip looks straight past them. The byte planes, on real data This is a slice of a real ERA5 2 m temperature field (201001_2t.nc, 24×32, detected E = −4). Left is the raw float32 byte planes, right is after lattice representation + differencing + zigzag: Byte plane raw float32, zero share after encoding lowest byte 100.0% 88.4% byte 1 5.2% 100% byte 2 0.0% 100% highest byte 0.0% 100% Note that on the left the lowest byte is already 100% zero — that is the quantization footprint — but the high bytes are pure entropy, so gzip has nothing to work with. After encoding it flips: the top three planes go to zero, and only 89 of the 768 differenced residuals are still non-zero. The pipeline float array → lattice representation (or order-preserving fallback) → 3-D separable Lorenzo differencing (time + lat + lon) → zigzag → byte-plane shuffle → zstd Lattice representation: turns 32-bit floats into integers with a far smaller dynamic range. This is where most of the win is. Differencing: assumes the field is continuous in time and space. First-order differences along time, then latitude, then longitude. The time lag is configurable — 24 to align with the diurnal cycle on hourly data. zigzag: differences are signed, and negatives are all-ones in two's complement. zigzag interleaves ± into non-negative values. Byte-plane shuffle: split each int32 into its 4 bytes and group like bytes together. zstd: finish up. How losslessness is guaranteed With quant_bits=0, decode(encode(x)) is byte-for-byte equal to x, including -0.0, NaN, and Inf. The guarantee is not "detect the odd values and special-case them". It is that both encoding paths are bijections on bit patterns: the lattice path relies on "power-of-two scaling is exact + per-block round-trip check", and the order-preserving path is a pure bit-pattern permutation. NaN/Inf are rejected by lattice detection outright and take the order-preserving path, so their bit patterns survive untouched. -0.0 is a separate trap: in the integer representation, -0.0 and +0.0 are both 0. So a sparse bitmap is stored alongside, and the sign bit is put back at decode time. The design principle throughout: prefer the slow path over a wrong answer. Results The tool's own output (gridpack info, real run): Item Result ERA5 coverage scan 834 / 834 files, zero exceptions Lattice hit rate 100% across all 834 files, no fallbacks Disk saving (834-file benchmark) 1.88× measured, 46.8% Production migration (measured) 151 files: 83.44 GB .nc → 44.70 GB .gpz, 1.87×, 46.4% saved vs raw float32 median 10.27× roundtrip / xarray tests 46 + 10 items, ALL PASS The production migration is the number I trust most, because it is the tool running unattended on real files, and it reproduces the benchmark independently: Sample Files .nc .gpz Ratio Saved migration, first run 49 26.23 GB 14.04 GB 1.87× 46.5% migration, after restart 94 50.37 GB 26.94 GB 1.87× 46.5% all completed so far 151 83.44 GB 44.70 GB 1.87× 46.4% Three separate samples, all 1.87×. The benchmark said 1.88×; production says 1.87×. Whole-file, end to end, including the decode round-trip: File source .nc .gpz vs .nc vs raw bytes 200003_10u.nc 589.1 MB 298.5 MB 1.97× 10.35× 200003_2t.nc 483.1 MB 282.9 MB 1.71× 10.92× 200003_tp.nc 189.0 MB 108.8 MB 1.74× 28.40× 200001_10u.nc 1443.5 MB 747.0 MB 1.93× 3.94× radar.nc 181.1 MB 6.3 MB 28.83× 28.60× Performance: v0.4.0 → v0.5.3 In the previous version of this write-up we listed three suspected optimizations and admitted we had not done them. We have now done them. The headline first: 1.56× faster end to end, with byte-identical output. Of three suspects, only one was real Instrumentation had pointed at the auto mode, detect_lattice, and the zstd settings. In practice detect_lattice was the only confirmed bottleneck. It computes the trailing-zero count tz(sig) once per non-zero value, and the original implementation was a per-value loop with an np.where branch. Replacing it with a 16-bit lookup table: Stage v0.4.0 v0.5.3 Speedup _block_exp (598 MB of values) 13.96 s 4.96 s 2.81× detect_lattice (whole file) 19.6 s 8.2 s 2.4× compress_file (end to end) 148.6 s 95.3 s 1.56× detect_lattice's share of end-to-end wall clock dropped from 13% to 8.6%. We changed the speed, not the format In a lossless compressor, the risk when optimizing is that you quietly change the output too. So every step was gated on three checks: Output bytes identical: 296.1 MB, ratio 1.631×, verify_file exact=True Trailing-zero function exhaustively verified: every value in 1..2²⁰ plus 2×10⁵ random 32-bit values, matching the original implementation exactly Real data unchanged: detected E still −4; all 4 modes round-trip bitexact; all 56 tests pass CODEC_VERSION stays at 4 — files written by v0.4.0 read in v0.5.4, and vice versa. The release process changed too Since 0.5.3 it goes through GitHub Actions: push a tag → version gate → build + twine check + test suite → OIDC tokenless publish to PyPI. 0.5.4 came out of that pipeline. pip install "gridpack[xarray]>=0.5.4" While we were there, 0.5.4 fixed a real bug: gridpack info and gridpack verify crashed on every single invocation in 0.5.3. main() read args.chunk_t unconditionally, but those two subcommands' parsers never defined it — an immediate AttributeError. The fix routes the read through getattr(..., None). A command documented in your README that cannot run is worse than a slow one. Limitations, stated plainly 1. It is not a general-purpose compressor. The lattice win depends entirely on the data actually being on a quantization lattice. On-lattice: large win. Off-lattice: it falls back to the order-preserving path, which degrades to "differencing + byte-plane shuffle + zstd" — roughly what any general method does. 2. Full-precision months drop to around 2×. 4 of the 834 files (200001/200002 10u/10v) have a lattice exponent of −16 rather than −4 — those months were not quantized; they are full-precision data. They reach only 3.84–3.86× vs raw bytes, about 1.9× vs the .nc on disk. Not a bug, a property of the data, but it drags the average down. 3. It is slow. Even 1.56× faster, it is still an hours-scale job. Measured single-process throughput on one ERA5 field file is roughly 110–190 s; the completed 26.23 GB took 2 h 16 min, which works out to about 3.2 MB/s. Extrapolating linearly, the remaining 319 GB needs about 28 hours. Going faster means multiple processes, and whether that works depends on the memory allowance your environment gives you, not on the compressor. We hit MemoryError: Unable to allocate 63.4 MiB in one restricted environment: the machine had 35 GB free, but the shared allowance for the whole subprocess batch was exhausted. Single-process in the same environment ran fine. 4. A lossy mode exists, and is off by default. quant_bits>0 discards low bits of the order-preserving integer. Default is 0: strictly lossless. 5. It has mainly been validated on ERA5. Plus one radar file. Whether other reanalysis products, model output, or GRIB2 land on a lattice, we have not tested. Usage pip install "gridpack[xarray]>=0.5.4" gridpack compress input.nc -o output.gpz gridpack decompress output.gpz -o restored.nc xarray plugin (registers engine="gpz"): import xarray as xr ds = xr.open_dataset("output.gpz", engine="gpz") ds.gpz.to_file("another.gpz") Source: https://github.com/showkeyjar/gridpack One caveat for non-Chinese users: the CLI output and the README are currently Chinese-only. That is a real barrier and it is on the list to fix. What I would like you to attack v0.4.0 is frozen pending criticism; performance work goes in 0.5.x (now at 0.5.4, format-compatible with 0.4.0 in both directions). The four questions I most want answered: 1. What other data sits on a lattice? We have only validated ERA5. GRIB2, other reanalysis products, model output, radar base data — if you have a sample, we would like to know the result. 2. What should happen to full-precision months? The E=−16 files get none of the quantization dividend. Should they be detected and routed elsewhere, or just not compressed at all? 3. How much is left in the entropy coder? Right now the byte planes go through zstd. zstd is a general-purpose compressor and has no idea it is compressing floating-point residuals. How much would a context-mixing model add? We have not measured it. 4. How else can performance be recovered? detect_lattice is handled (1.56×). Two suspects remain: auto mode tries three representations at 2.4× the cost of a single mode, and in practice it usually picks dt for zero benefit; and verification reads the whole file back. If you have a better idea, say so. Data samples, benchmark results, or direct criticism of the design — issues are welcome. If this approach is useful to you, or if you have data that disproves it, we would like to hear about it.