Tighten the stored format, open each file once, and read what ObsPy reads - #82
Merged
Conversation
open_dataarray passed phony_dims="sort" to silence xarray's transitional phony_dims warning (0f43721). But "sort" makes h5netcdf run _determine_phony_dimensions(), an eager recursive walk of every group in the file, on each open -- and open_datacollection calls open_dataarray once per array. On a 42-array collection that is 42 full-tree walks of 95 groups and 787 datasets: 11.6 s to open, 8.2 s of it in that scan. "access" is xarray's own default, names phony dims lazily per group, skips the walk, and silences the warning just the same. No xdas-written file has unlabeled dimensions, so no phony dim is ever created either way.
A geometry column holding a single value — the absent `starts_k` and `steps_k`, and the `sizes_k` of an acquisition of equal-length files — is detected at construction and kept as a zero-stride broadcast view instead of a full-length array. Their tile boundaries follow: `_edges` becomes a small per-axis object holding either the eager `cumsum` of a varying column, exactly as before, or the constant tile width, from which `total`, `at`, `searchsorted` and `span` are closed forms. Detected, never assumed: `combine` starts a new acquisition whenever the geometry changes, so the size column cannot vary inside a group, but a truncated last file would give a second value and must stay on the eager path. The manifest dataset is untouched, so a round trip rewrites exactly what it read. On a 23-million-tile archive, opening drops from 1.67 GB to 1.11 GB resident and from 3.70 s to 3.31 s; locating a tile along a constant axis is a division instead of a binary search (2.4 us to 0.3 us).
Reading walks a single xr.open_datatree of the whole file instead of one targeted xr.open_dataset per group; writing lands every group's metadata through a single DataTree.to_netcdf and every data variable through one writable h5netcdf handle. One writable handle is the fix for #81: h5netcdf walks the whole file on every writable open to find the next free dimension id, so per-array opens made saving quadratic in the collection size — now flat at ~23 ms per event at N=800. The manifest of a tile-backed array becomes just another node of the same tree, so VirtualBackend.finalize_save (a reopen per array) is replaced by sibling_datasets, which returns the stored form instead of writing it. _get_depth and its visit() walk retire with the reader: a node holding variables is a data array, one holding only groups is a nesting level. Deviation from the plan sketch: eager variables are written in the h5netcdf pass rather than through to_netcdf, which would have changed the accepted encoding keys (h5netcdf's "chunks" and hdf5plugin filter dicts vs xarray's "chunksizes").
A geometry column holding the same value for every tile was folded to a zero-stride view in memory but still written, read and held per tile: 185 MB of `sizes_0` on a 23 M-tile archive, all of it the number 6. Fold it in the stored manifest too. A column holding its kind's default is left out, one holding a single value everywhere is stored 0-d, and only a varying one stays 1-D. Folding a column stops it sizing its axis, and nothing else is obliged to: `concat([a, a])` on a single-file array folds every column, and would reopen as one tile instead of two. So a folded `sizes_k` carries an `ntiles` attribute — the manifest states the shape of an axis it no longer measures, instead of leaving it to be inferred from whichever variable happens to still carry the dimension. It rides on `sizes_k` rather than on the dataset, whose every attribute is a user attribute, and being an attribute it costs no HDF5 object, unlike a variable or a dimension of its own. Folding happens in the constructor, beside the recoding of strings to fixed-width bytes: every entry point funnels through it, so scans, slices, flips, concatenations and stored manifests all converge on the same form, and a manifest written before this change folds on open instead of needing a rescan. On a 42-array, 23.07 M-tile archive: 765.0 -> 579.7 MB on disk, 1.112 -> 0.918 GB resident after open, 3.32 -> 2.94 s to open.
The stored manifest spread what a tiling needs over four spellings: an `ntiles` attribute on each folded `sizes_k`, a `root` variable with its own string-width dimension, an `axes`/`source_ndim` pair of variables with their own dimension, and the engine on the placeholder in a `__tile_array__` attribute -- while the dtype was nowhere, inferred from the placeholder's HDF5 type. Collect all of it into one JSON document, the `header` attribute of the `__tiles__` group: `ntiles`, `engine`, `dtype`, and `root` and `axes` absent at their defaults. One `json.loads`, natively typed inside, and the manifest's variables become exactly the per-tile columns. `root` rides through JSON rather than as a bytes attribute, which h5netcdf rejects outright: `json.dumps(os.fsdecode(root))` is pure ASCII even for a non-UTF-8 directory, and `os.fsencode` recovers the bytes exactly. `axes` becomes the full virtual arrangement in numpy indexing-key notation, `null` where an axis is inserted -- `[1, null, 0]` reads as `transpose(a)[:, np.newaxis]` -- so the source rank stops being a separate statement: it is the grid rank less the number of nulls. That costs two invariants, both enforced in `_canonical`: synthetic axes stay trailing and numbered in their virtual order, so a transpose that reorders them renumbers the grid instead of the map; and an integer key on a synthetic axis deletes it from the grid rather than hiding it, a hidden synthetic axis being contentless once the pin has selected into `paths`. Manifests are canonical either way, and `TileArray(dataset)` now takes one argument -- `dtype` and `engine` stay as optional arguments that must agree with the header when there is one. The array itself stops carrying user attributes. It is a duck array; numpy arrays have no `.attrs`, no scanner ever set them, and the metadata belongs to the enclosing DataArray, which is where it already was in practice. The placeholder now points at its describing group with a `__tiling__ = "__tiles__"` attribute, in the fashion of CF's `grid_mapping`: dispatch is attribute-driven, the group name demotes to a conventional default, and several tiled variables could share a file. It is dunder-fenced because the placeholder's attributes are the DataArray's user namespace -- the file has exactly two reservation registries, CF's plain words and xdas's `__*__`, and the reader strips the latter by pattern. While the CF face is being rewritten, make it true. The old spelling was CF-shaped but did not validate: `coordinate_interpolation` groups must end with the interpolation variable, the mapping attribute is the singular `tie_point_mapping` naming the interpolated dimension, the index variable and the subsampled dimension, and `computational_precision` is required. Files declare CF-1.13. The reader accepts both grammars, telling them apart by whether a group's last word names a variable carrying `interpolation_name`. One deliberate break, so files change once: everything written before still opens and upgrades in place, but files written now cannot be read by earlier versions. On the 42-array, 23.07 M-tile archive the rewrite is byte-for-byte equivalent -- same arrays, same geometry, same reads -- and opening goes from 3.22 s / 1.72 GB to 2.57 s / 0.94 GB, the drop being the pre-fold columns the reader no longer has to fold on the way in.
The sampled format mirrored CF-1.13 coordinate subsampling only loosely: it carried the spacing as the sampling variable's own value, with private `dtype` and `units` attributes to rebuild a timedelta, and its `tie_point_mapping` listed the tie point coordinate variable where the grammar expects the index variable. Word it like the interpolated case instead. The sampling variable is now a container, as CF's interpolation variable is, and its mapping puts the segment length variable in the tie point index variable's slot — the one deliberate departure, since a sampled axis describes its segments by length rather than by end index. The spacing travels as attributes through the shared `encode_delta` / `decode_delta`, the way an interpolated coordinate already carries its regular metadata. No file on disk uses the previous spelling, so nothing reads it. Also stop the interpolated reader raising a `KeyError` on the oldest files, which spelled the mapping without writing any interpolation variable at all.
The header made the break one-way in intent, but the reader still carried the whole earlier spelling: the `__tile_array__` placeholder attribute, the `root` / `axes` / `source_ndim` variables, the per-column `ntiles` attributes, and the canonicalization that turned them into a header. That is a second format to keep correct for files that a single rewrite with 0.2.8 converts. The upgrade path goes; what remains is what a dataset with no header says by itself: the tile counts its declared dimensions give, the identity arrangement, the paths stored whole. That is a hand-assembled manifest, which the scan-time assembly and the tests build, not a stored one — every file this release writes carries its header.
No engine has emitted a dask graph since tiles landed: the two formats that needed one, Silixa and MiniSEED, produce tile manifests, which describe the same mapping as plain array data instead of as a computation graph, and build in a fraction of the time. What was left was a serializer, a deserializer and a deprecated write path, all of them dead weight on the format and on the writer. The `xdas.dask` module goes, with the `msgpack` dependency it alone needed, and `__dask_array__` is no longer read. A dask array is still valid data for a `DataArray` — it is now written like any other lazy value, computed on the spot, and `virtual=True` rejects it as it does any other non-virtual array.
`trim_overlaps` cuts a data array at its overlaps and keeps one copy of each duplicated sample, lazily and on sample boundaries. Claims accumulate as spans rather than a watermark, so a segment enveloped in a lower-precedence one keeps a run on each side of it instead of losing everything past the overlap. `concat` opening a new dimension now checks that the inputs agree on their other coordinates, and promotes the scalar ones that vary to a coordinate along that dimension. `DataCollection.query` reports fields it does not know instead of silently returning everything, `fields` walks the whole subtree, and `select` names the operation the way obspy does.
The miniseed engine described a whole file as one tile and classified it at scan time as synchronized or unsynchronized, which rejected files obspy reads fine, left no pointer to an individual trace, and turned overlap handling into a destructive global flag baked into every stored manifest. The obspy engine emits one lazy data array per obspy Trace instead, addressed by the data's own address — the four SEED identifiers plus both time bounds — so a re-segmenting obspy version cannot silently designate different samples. The collection nests on the SEED hierarchy, giving Stream.select semantics through DataCollection.select. Merging contiguous traces, moving gaps into the coordinate and separating acquisition epochs are combine_by_coords' job, which `open` now runs whether it opened one file or many. Named for the library, not the format: engine="miniseed" stays as an alias and pre-rename manifests keep decoding, but everything obspy reads now goes through the same path. The engine is registered last so it does not shadow the format-specific ones during auto-detection. Along the way: a mapping keyed by a zero-padded code — a SEED location — no longer reads back from netCDF as a sequence, and `query` applies an indexer wherever its level sits rather than only at the root.
The miniSEED guide becomes the obspy guide, rewritten around what the engine now returns: read a glob without describing the directory layout, select by SEED level, recover the original segments from a gap, resolve overlaps, then stack channels and stations into an N-dimensional array — all of it lazy.
Restored verbatim under its own name, so views it wrote keep decoding and code written against it keeps running unchanged — `ignore_last_sample` and the `dim="station"` idiom included. It is no longer an alias for the obspy engine: the two describe the same files with different shapes, so one name cannot mean both. Both take part in auto-detection, `obspy` registered first. `xd.open` asks for a collection, which only the obspy engine provides, so it wins; `xd.open_dataarray` asks for a single array, which the obspy engine cannot give for a multi-trace file, and the legacy engine stacks it. Each function keeps returning what it is named for. The legacy decode moves out of `ObsPyEngine.load_tile` and back onto the engine that wrote it, which loses its `method`/`ignore_last_sample` branch. The stream converters and the band-code table were never miniSEED-specific: they stay with the obspy engine and are re-exported.
Each entry is one or two sentences, as 0.2.8 reads: the rationale belongs in the docstrings and the user guide, which already carry it. Three substantive corrections along the way: the element type read from the miniSEED encoding is a property of the new obspy engine, not a fix to the legacy one, so it moves into that entry; sampled coordinates stored by 0.2.8 are not read back, which the notes did not say; and the collection I/O that is linear in its size again (#81) was missing altogether.
Dask virtualization is gone; the only thing left was one isinstance in the data array repr, holding a pinned `dask<2025.4.0` in the runtime dependencies. The repr duck-types the module name instead, and dask moves to the test group, where the two tests that build one live.
The sequence branch of `query` accepts an integer or a slice and told the caller it wanted a string, which is the mapping branch's rule.
`headonly=True` leaves `tr.data` an empty float64 array whatever the file holds, so a STEIM-compressed int32 file was scanned as float64. The dask path never checked, and delivered int32 under a float64 label; the tile array does check, so the read failed outright — STEIM being the common case, `engine="miniseed"` was unusable on it.
Interpolated coordinates kept a compatibility branch across the format change and sampled ones did not, so a 0.2.8 file with `ctype="sampled"` failed with a bare KeyError naming a variable the reader had invented. The two spellings are told apart the same way as for interpolation: by whether the container carries the attribute the new grammar puts there.
The rule had three copies, and the fix for zero-padded keys reached two of them: `xd.open_datacollection(path, engine="xdas")` still parsed the keys as integers and raised KeyError on a collection keyed by a SEED location. The nested reader was a single-use wrapper around the same rule, so it folds into its caller.
The sweep discarded every exception, so a file an engine recognised and then rejected for a stateable reason was reported as merely unreadable. No single refusal is the reason — engines are tried in a cache-warmed order, not in order of likelihood — so the message lists them all, and leaves out the engines that only said they do not offer that shape.
`from_stream` built its start time from `.datetime`, truncating to microseconds where the rest of the engine works in nanoseconds; `_bare_header` was a single-use helper whose caller already explains itself; and a doc code block was written so that `ruff format` moved a comment onto its own line.
It was dropped when the obspy engine replaced the miniseed one, but ffb95f6 restored the legacy engine with the parameter intact, and the breaking-change entry promises it still works.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## dev #82 +/- ##
==========================================
Coverage 100.00% 100.00%
==========================================
Files 47 45 -2
Lines 6017 6252 +235
Branches 1016 1067 +51
==========================================
+ Hits 6017 6252 +235 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
The stored format
One header per tiling. Everything about a tiling that is not a per-tile
column — the tile counts, the engine specification, the element type, the
common source directory, the axis arrangement — now travels in a single JSON
headerattribute on the__tiles__group, replacing the__tile_array__placeholder attribute, the
root/axes/source_ndimvariables and theper-column
ntilesattributes. The manifest variables are now exactly theper-tile columns, and the placeholder points at its group through a
CF-
grid_mapping-style__tiling__attribute instead of being tied to thegroup name. The reader does not accept the earlier spelling — see the
migration note below.
Constant tile geometry costs one element, not one per tile. A
sizes_k,starts_korsteps_kcolumn holding a single value — what an acquisition ofequal-length files gives, and always the case for the absent origin and stride
columns — is kept as a broadcast view in memory and a 0-d variable on disk, and
its tile boundaries as a closed form instead of a full
cumsum. On a23-million-tile archive: 1.67 GB → 1.11 GB resident, and locating a tile along
such an axis becomes a division instead of a binary search.
Coordinates are stored as CF-1.13 actually defines them. Interpolated
coordinates get the singular
tie_point_mappingand the mandatorycomputational_precision; sampled ones are stored as a deliberate variation onthe same grammar. Both still read the earlier spelling.
Dask virtualization is gone, reader and writer alike, along with the
xdas.daskmodule. No engine has emitted it since tiles landed. A Dask arrayremains valid
DataArraydata — it is computed on write like any other eagerarray — and this removes the
FutureWarningdeprecation that was indev.Opening and saving
One open per direction — fixes #81. Reading walks a single
xr.open_datatree; writing lands every group's metadata through oneDataTree.to_netcdfand the data variables through one writableh5netcdfhandle. That matters on write because every writable
h5netcdfopen walks thewhole file to find the next free dimension id, so one open per data array made
saving a collection quadratic in its size.
DataCollection.to_netcdf(virtual=True), the case #81 reports:devdevdoubles its per-event cost every time the collection doubles; this branchis flat, so the gap keeps widening with size.
The open path benefits from the same shared handle (595 s → 44 s at 400
groups) but is not linear yet — 35 → 110 ms/event over that range. The
residual has a different, independent cause:
h5netcdfresolvesDIMENSION_LISTreferences throughH5Iget_name, which is O(objects in file)per group and sits outside xdas. Worth its own issue; this PR does not claim it.
Reading seismological files
The
obspyengine, named for the library rather than for a format: decodingis
obspy.read, so miniSEED, SAC, GSE2, SEG-2 and everything else ObsPysupports goes through it. Each contiguous
Tracebecomes one lazyDataArrayand the collection mirrors the
Stream, nestednetwork / station / location / channel. Files the old engine rejected (twosampling rates, duplicated ids, interleaved acquisitions) are readable, and each
tile points at an individual trace instead of the whole file. The
"miniseed"engine is kept unchanged beside it — views it wrote keep decoding, code written
against it keeps running,
ignore_last_sampleincluded — and auto-detectionreaches
"obspy"first.Supporting changes:
xdas.trim_overlapsresolves overlaps at the manifestlevel (keeping the later copy by default, ObsPy's
merge(method=1));xdas.concatopening a new dimension now checks the other coordinates agreeand promotes the scalar ones that vary, so
xd.concat(traces, "channel")givesback the stacked shape lazily;
DataCollection.select(alias ofquery) anda
fieldsthat reports the whole subtree giveobspy.Stream.selectsemantics.Breaking changes
write them back with this branch. Stored coordinates are unaffected — both
spellings still read.
xd.openon a seismological file returns a nested collection, not astacked array.
xd.concat(traces, "channel")is the one-liner back;engine="miniseed"still gives the old shape.xd.opennow combines whether it opened one file or many, so the returnedshape no longer depends on the file count; leaf sequences are named
acquisition.DataCollection.queryraisesKeyErroron an indexer naming no level of thecollection, and applies an indexer wherever its level sits in the tree.
TileArraycarries no user attributes."--".