perf(collection): stream DatasetCollection.to_netcdf instead of materialising the full cube

[MAfarrag]

Context

DatasetCollection.to_netcdf writes the collection's (T, B, Y, X) cube to a single NetCDF via pyramids' own
GDAL multidim writer (no third-party engine). It is eager: it materialises the whole cube in RAM before writing.
The docstring and an ARC-46 inline warning already flag this and point users at to_zarr for large cubes. This
issue tracks making the writer streaming so to_netcdf no longer needs the full cube resident.

Related history:

• feat(collection): atomic zarr append + validate=True headers; track LazyDatasetCollection #874 (closed) — atomic/streamed to_zarr append + LazyDatasetCollection tracking; same lazy-datacube theme.
• ARC-46 (planning/architecture-review — arc-collection.md) — the eager-materialisation finding this addresses.
• to_zarr already streams chunk-by-chunk; this brings to_netcdf to parity.

Problem / Current Behaviour

Peak memory is roughly 2× the full cube, from two layers of eagerness:

1. In to_netcdf — the timesteps are stacked into one array:

   # src/pyramids/dataset/collection.py (~L1597)
   cube = np.stack([np.asarray(ds.read_array()) for ds in self.datasets], axis=0)

   self.datasets are already lazy per-timestep handles and read_array() reads only that one raster — so it is
   purely this np.stack that forces every timestep into RAM at once.

2. In the writer — interop._build_multidim writes the full arrays into an in-memory MEM multidim dataset
   (md_arr.Write(full_array)), then _create_copy_to_netcdf does netCDF.CreateCopy(path, mem_src) — a second
   full copy of the cube.

The existing ARC-46 guard (collection.py ~L1578) only warns above a size threshold; it does not avoid the
allocation.

Affected locations

File	Symbol	Notes
src/pyramids/dataset/collection.py	DatasetCollection.to_netcdf	np.stack of all timesteps (~L1597); ARC-46 warn (~L1578)
src/pyramids/netcdf/engines/interop.py	_build_multidim (~L353)	writes full arrays into a MEM multidim dataset
src/pyramids/netcdf/engines/interop.py	_create_copy_to_netcdf (~L463)	CreateCopy MEM → netCDF (second full copy)
src/pyramids/netcdf/engines/interop.py	write_multidim_netcdf (~L483)	public entrypoint; takes fully-materialised data_vars

Feasibility — empirically confirmed (GDAL 3.13.1, this env)

GDAL's netCDF driver reports DCAP_CREATE_MULTIDIMENSIONAL = YES, and a multidim array can be created at full shape
and written slab by slab (one timestep per write), never holding the whole cube:

ds  = gdal.GetDriverByName("netCDF").CreateMultiDimensional(path)
rg  = ds.GetRootGroup()
arr = rg.CreateMDArray("data", [dt, dy, dx], edt)   # full shape, no data yet
for t in range(T):
    arr.Write(one_timestep, array_start_idx=[t, 0, 0], count=[1, Y, X])   # returns rc=0
ds.Close()

Round-trip verified: 3 timesteps written as [1, Y, X] slabs, reopened → shape (3, 4, 5) with correct per-timestep
values. Gotchas found:

• Reopen needs the netCDF driver hint (gdal.OpenEx(path, gdal.OF_MULTIDIM_RASTER, allowed_drivers=["netCDF"]));
  pyramids' NetCDF.read_file already opens netCDF correctly, so consumers are unaffected — this was only a
  test-harness detail.
• FORMAT=NC4C is rejected ("Operation not allowed in define mode"); default and FORMAT=NC4 both work.

Proposed Solution

Add a streaming writer alongside write_multidim_netcdf, and rework to_netcdf to feed it timestep-by-timestep:

• New writer creates the netCDF multidim file directly via GetDriverByName("netCDF").CreateMultiDimensional(path),
  writes the small 1-D coordinate arrays (time/y/x) whole, sets root/variable attributes via SetAttribute
  (reuse _apply_md_array_attrs), and exposes the data MDArray(s) for the caller to fill by slab.
• to_netcdf replaces the np.stack with a loop over enumerate(self.datasets), writing each timestep's slab:
  [t, :, :] into each per-band variable (var_per_band=True), or [t, :, :, :] into the single 4-D data
  variable (var_per_band=False).

Result: peak RAM drops from ~2×(T·B·Y·X) to ~one timestep (B·Y·X) + the coordinate arrays. Output file is
byte-for-byte equivalent in structure (same dims/vars/attrs), only the write path changes.

Out of Scope

• The to_zarr path (already streaming) and the broader LazyDatasetCollection (see closed feat(collection): atomic zarr append + validate=True headers; track LazyDatasetCollection #874).
• Any change to the on-disk NetCDF layout, variable names, or CF attributes — behaviour-preserving output only.

Effort Estimate

Size: M — one new writer function + the to_netcdf loop + tests. No new dependency (pure GDAL multidim).

Definition of Done

• a streaming multidim writer creates the netCDF and fills the data array by per-timestep slab
• to_netcdf no longer builds the full (T, B, Y, X) cube (no np.stack of all timesteps; no MEM CreateCopy)
• output is unchanged: same dims/vars, and nodata / crs_wkt / epsg / GeoTransform / Conventions round-trip
• both var_per_band=True and var_per_band=False paths stream correctly
• tests assert the round-trip equals the previous eager output and cover the time-axis encodings (int index +
  datetime64 → CF nanoseconds since units)
• peak-memory behaviour documented (drop the ARC-46 "materialises the full cube" caveat from the docstring)