A netCDF memory spike, found and fixed¶

A worked story from the PyPSA ecosystem — [issue #1555][issue], fixed in [PR #1830][pr]. Reading a network that carried shape geometries back from netCDF made memory spike ~10×, big enough to OOM-kill loads that used to fit. The timings barely moved; the peak is what gave it away.

Everything on this page runs live when the docs build. Here is the whole example — a length-skewed set of WKT polygons (many tiny regions, a few detailed coastlines) read two ways. Both return the same strings; only the dtype differs:

In [1]:

import numpy as np


def make_wkt_geometries(n_small: int = 500, n_detailed: int = 5, n_points: int = 6_000) -> list[str]:
    """A length-skewed set of WKT polygons: many tiny regions and a few detailed
    ones — the distribution that makes fixed-width padding pathological (the PR's
    own repro uses 200k points / ~6 GB; here it's dialled down to stay laptop-sized)."""
    tiny = ["POLYGON ((0 0, 1 0, 1 1, 0 0))"] * n_small
    detailed = ["POLYGON ((" + "0 0," * n_points + "0 0))"] * n_detailed
    return tiny + detailed


def to_fixed_width_array(geometries: list[str]) -> np.ndarray:
    """Before #1830: a fixed-width Unicode array. NumPy sizes the dtype to the
    longest string and pads every entry to it — n × 4 × len(longest) bytes."""
    return np.array(geometries, dtype="U")


def to_object_array(geometries: list[str]) -> np.ndarray:
    """After #1830: an object array. Each entry stays its own variable-length
    string, so the array costs the sum of the actual lengths."""
    return np.array(geometries, dtype=object)

import numpy as np def make_wkt_geometries(n_small: int = 500, n_detailed: int = 5, n_points: int = 6_000) -> list[str]: """A length-skewed set of WKT polygons: many tiny regions and a few detailed ones — the distribution that makes fixed-width padding pathological (the PR's own repro uses 200k points / ~6 GB; here it's dialled down to stay laptop-sized).""" tiny = ["POLYGON ((0 0, 1 0, 1 1, 0 0))"] * n_small detailed = ["POLYGON ((" + "0 0," * n_points + "0 0))"] * n_detailed return tiny + detailed def to_fixed_width_array(geometries: list[str]) -> np.ndarray: """Before #1830: a fixed-width Unicode array. NumPy sizes the dtype to the longest string and pads every entry to it — n × 4 × len(longest) bytes.""" return np.array(geometries, dtype="U") def to_object_array(geometries: list[str]) -> np.ndarray: """After #1830: an object array. Each entry stays its own variable-length string, so the array costs the sum of the actual lengths.""" return np.array(geometries, dtype=object)

Discovery — the peak flags the culprit¶

Point a benchmark at both readers and add --benchmark-memory. It's an ordinary pytest-benchmark test — the only new thing is the flag:

READERS = {"fixed-width (before #1830)": to_fixed_width_array, "object (after #1830)": to_object_array}

@pytest.mark.parametrize("reader", READERS, ids=list(READERS))
def test_read_shapes(benchmark, reader):
    benchmark(READERS[reader], make_wkt_geometries())

The timing table is business as usual; the peak columns, right of the divider, are not:

benchmark: 2 tests                                                                                                 
                                                                                                                   
  Name (time in us)                                             Mean   │   peak·min (KiB)   peak·mean    peak·max  
 ───────────────────────────────────────────────────────────────────────────────────────────────────────────────── 
  test_read_shapes[object (after #1830)]                8.4878 (1.0)   │             3.95        3.95        3.95  
  test_read_shapes[fixed-width (before #1830)]   3,084.3090 (363.38)   │        47,373.34   47,373.34   47,373.34  
                                                                                                                   
memory (right of │): a separate, untimed pass, not the timed rounds  •  also available via                         
--benchmark-memory-columns: allocated, allocations                                                                 
2 passed in 2.54s

Same 505 strings in, same strings out — and a 46 MiB peak against 4 KiB. In the real network that ratio was the difference between loading and dying. A timing-only suite would have shrugged: the mean is slower, but "slower" isn't the symptom users hit — the OOM is, and only the peak column names it.

Examination — why the array is 300× its payload¶

The peak says which reader is heavy; look at what it built to see why. Nothing memray-specific here — just the array's own dtype:

In [4]:

geometries = make_wkt_geometries()
arr = to_fixed_width_array(geometries)

print("dtype           :", arr.dtype)                          # sized to the longest string
print("array           :", round(arr.nbytes / 1e6, 1), "MB")
print("actual payload  :", round(sum(len(g) for g in geometries) / 1e6, 2), "MB")

geometries = make_wkt_geometries() arr = to_fixed_width_array(geometries) print("dtype :", arr.dtype) # sized to the longest string print("array :", round(arr.nbytes / 1e6, 1), "MB") print("actual payload :", round(sum(len(g) for g in geometries) / 1e6, 2), "MB")

dtype           : <U24015
array           : 48.5 MB
actual payload  : 0.14 MB

There it is. A fixed-width Unicode array sizes its dtype to the longest string — here one 24,015-character coastline — and pads every entry to it, at 4 bytes per character (UTF-32). So 505 strings whose real payload is 0.14 MB inflate to 505 × 4 × 24015 ≈ 48.5 MB. The tiny rectangles are each paying the coastline's weight. On the full network, with a 200k-point shape, that same padding is gigabytes.

Fix — object array, confirmed with compare¶

[PR #1830][pr]'s fix is one dtype: keep the strings as variable-length Python objects, so the array costs the sum of the real lengths, not the longest repeated. Run the reader before and after into two JSON files and diff them:

benchmem compare before.json after.json --columns peak --diff

bench_read.py::test_read_shapes
 name               peak base     after 
────────────────────────────────────────
 test_read_shapes    46.3 MiB   -100.0% 

The peak is gone. Wire that same diff into CI with --fail-on peak:20% (Catch regressions in CI) and the spike can't come back unnoticed — the next PR that reintroduces fixed-width padding turns the build red instead of shipping an OOM.

The loop, on your own code¶

Nothing here is shape-specific. Any allocation whose size dwarfs its payload — a padded array, a densified sparse matrix, a duplicated frame, a solver's working set — shows up the same way: a peak that timing can't see, a heavy line you can point at, and a compare that proves the fix. Start at the Quickstart with a benchmark you already have.