From 087386db51904fac046616aa0b47163bec231296 Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Sat, 15 Aug 2026 17:48:10 +0200 Subject: [PATCH 1/5] Let a geometry state named numeric columns, not only position An optical path could state a value varying along distance only if that value was a CRS position, and a CRS holds at most three axes. Borehole depth where the CRS is spent on easting, northing, and elevation; pipeline chainage; burial depth; fiber azimuth -- all of them fell in the gap, and the only place left for them was annotations, which are a set of intervals rather than a curve and hold no unit. Geometry.coordinates becomes a mapping of column name to values, with a units mapping beside it for the columns which are not axes. Which ones those are is the CRS's to say: a column it declares, or the canonical x/y/z alias of one, is that axis and takes the CRS's units; anything else is a quantity in its own right. A segment states every axis or none of them, so half a position never reaches a reader, and a segment which states no axis at all is now a perfectly ordinary thing to write. Each column is its own function track. Two segments may cover the same distance as long as they do not state the same column over it, which is what lets a borehole depth and a fiber azimuth be surveyed independently. Interpolation, half-open coverage, and the run-end rule are unchanged, and a column never bridges two segments. In a CSV, a header the CRS does not name is such a column and may carry its units in parentheses -- `depth (m)`. Text is refused there, with an error pointing at annotations.csv, which is where a value that varies along the fiber without being a number belongs. --- dascore/core/_spool_inventory.py | 18 +- dascore/core/inventory.py | 340 ++++++++++++++++++----- dascore/core/inventory_loader.py | 136 ++++++--- dascore/examples.py | 8 +- dascore/proc/inventory.py | 16 +- docs/recipes/tunnel_inventory.qmd | 2 +- docs/tutorial/inventory.qmd | 28 +- tests/test_core/test_inventory.py | 257 ++++++++++++++--- tests/test_core/test_inventory_loader.py | 86 +++++- tests/test_proc/test_proc_inventory.py | 96 ++++++- 10 files changed, 817 insertions(+), 170 deletions(-) diff --git a/dascore/core/_spool_inventory.py b/dascore/core/_spool_inventory.py index 27be9dc65..629eac569 100644 --- a/dascore/core/_spool_inventory.py +++ b/dascore/core/_spool_inventory.py @@ -668,12 +668,20 @@ def _get_geometry_coord(inventory, path, label, distances): return None if not path.geometry: return None - coords = path.coordinates_at(distances) + coords = path.coordinates_at(distances, crs) if index >= coords.shape[1]: return None return get_coord(data=coords[:, index], units=crs.units[index]) +def _get_geometry_column_coord(path, name, distances): + """Return one geometry column which is not a position, with its units.""" + values = path.column_at(name, distances) + if values is None: + return None + return get_coord(data=values, units=path.column_units(name) or None) + + def get_coord_values(inventory, path, name, distances): """Return the values of one requested coordinate, or None if undefined.""" if name == "distance": @@ -688,7 +696,13 @@ def get_coord_values(inventory, path, name, distances): path, track, field or TRACK_IDENTITY_FIELDS[track], distances ) if name in VALID_COORDINATE_LABELS: - return _get_geometry_coord(inventory, path, name, distances) + # A coordinate label this CRS does not declare is not an axis here, + # and is free to be a geometry column of its own -- depth, where the + # CRS is spent on easting, northing, and elevation. + if (coord := _get_geometry_coord(inventory, path, name, distances)) is not None: + return coord + if (coord := _get_geometry_column_coord(path, name, distances)) is not None: + return coord return _get_annotation_coord(path, name, distances) diff --git a/dascore/core/inventory.py b/dascore/core/inventory.py index ca53d7cda..e2ef9dc86 100644 --- a/dascore/core/inventory.py +++ b/dascore/core/inventory.py @@ -16,6 +16,7 @@ import itertools from collections.abc import Mapping, Sized +from contextlib import suppress from functools import cache from types import MappingProxyType, UnionType from typing import ( @@ -509,43 +510,93 @@ class Terminator(_OpticalComponentBase): class Geometry(InventoryModel): """ - Geometry for an interval of an optical path. + Measured curves along an interval of an optical path. A geometry is a piecewise segment placed by its ``distance`` array: at - least two strictly increasing optical distances, each paired with the - coordinate at that point (interpreted using the inventory CRS). Coverage - is the half-open span of the array; there is no separate length field. A - coil, or other "clump", is a segment whose coordinates repeat while - distance advances. - Interpolation between points is piecewise linear in the CRS and never - crosses segments; uncovered distance has undefined coordinates. + least two strictly increasing optical distances, each paired with one + value of every column the segment states. Coverage is the half-open span + of the array; there is no separate length field. A coil, or other + "clump", is a segment whose columns repeat while distance advances. + + A column whose name the inventory CRS declares -- or the canonical + ``x``, ``y``, ``z`` alias of one -- is that position axis, and its units + are the CRS's. Every other column is a numeric quantity along the fiber + in its own right: borehole depth where the CRS is spent on + easting/northing/elevation, pipeline chainage, burial depth, fiber + azimuth. Those carry their own entry in ``units``. + + Interpolation between points is piecewise linear and never crosses + segments; uncovered distance has undefined values. + + Examples + -------- + >>> from dascore.core.inventory import Geometry + >>> + >>> # A surveyed run, placed by the CRS's axes. + >>> trench = Geometry( + ... name="trench", + ... distance=(0.0, 100.0), + ... coordinates={"x": (0.0, 86.6), "y": (0.0, 0.0), "z": (-0.5, -0.5)}, + ... ) + >>> + >>> # A curve which is not a position at all. + >>> chainage = Geometry( + ... name="chainage", + ... distance=(0.0, 100.0), + ... coordinates={"chainage": (1200.0, 1290.0)}, + ... units={"chainage": "m"}, + ... ) """ _identity_field: ClassVar[str] = "name" name: str = Field(default="", description="Human-readable geometry name.") distance: tuple[float, ...] = Field( description=( - "Optical distances paired to coordinates; at least two strictly " + "Optical distances paired to values; at least two strictly " "increasing values whose span is the segment's coverage." ), ) - coordinates: tuple[tuple[float, ...], ...] = Field( - description="Coordinate points; same length as distance.", + coordinates: FrozenDictType[str, tuple[float, ...]] = Field( + description=( + "Columns measured along this segment, keyed by name; each holds " + "one value per distance. A name the CRS declares is that " + "position axis, any other is a numeric column of its own." + ), + ) + units: FrozenDictType[str, str] = Field( + default_factory=dict, + description=( + "Units of the columns which are not position axes; the CRS " + "states the units of the axes." + ), ) @model_validator(mode="after") def _validate_geometry(self) -> Self: """Enforce paired, strictly increasing control points.""" _check_control_points(self.distance, "Geometry distance", minimum=2) - if len(self.coordinates) != len(self.distance): - msg = "Geometry coordinates and distance must have the same length." + if not self.coordinates: + msg = "Geometry states no columns, so it describes nothing." raise InvalidInventoryError(msg) - dims = {len(coord) for coord in self.coordinates} - if len(dims) > 1 or 0 in dims: - msg = "Geometry coordinate points must share one nonzero dimensionality." + wrong = sorted( + name + for name, values in self.coordinates.items() + if len(values) != len(self.distance) + ) + if wrong: + msg = ( + f"Geometry column(s) {wrong} do not have one value per " + f"distance; distance states {len(self.distance)}." + ) raise InvalidInventoryError(msg) - if not np.all(np.isfinite(np.asarray(self.coordinates, dtype=float))): - msg = "Geometry coordinate values must be finite." + for name, values in self.coordinates.items(): + if not np.all(np.isfinite(np.asarray(values, dtype=float))): + msg = f"Geometry column {name!r} must hold finite values." + raise InvalidInventoryError(msg) + # Units name a column or nothing; the alternative is a unit sitting + # on a column the segment never states, which no reader would find. + if orphaned := sorted(set(self.units) - set(self.coordinates)): + msg = f"Geometry states units for {orphaned}, which it has no column for." raise InvalidInventoryError(msg) return self @@ -554,23 +605,25 @@ def interval(self) -> tuple[float, float]: """The (start, end) optical distance covered by this segment.""" return (self.distance[0], self.distance[-1]) - def interpolate(self, distances) -> np.ndarray: + def interpolate(self, distances) -> dict[str, np.ndarray]: """ - Return coordinates at the requested optical distances. + Return each column's values at the requested optical distances. - Distances outside this segment's coverage return NaN rows. Coverage is + Distances outside this segment's coverage return NaN. Coverage is half-open ``[first, last)``; inclusion of the outermost track endpoint is handled by the caller (`OpticalPath.coordinates_at`). """ dist = np.atleast_1d(np.asarray(distances, dtype=float)) - coords = np.asarray(self.coordinates, dtype=float) - out = np.full((len(dist), coords.shape[1]), np.nan) start, end = self.interval inside = (dist >= start) & (dist < end) - for dim in range(coords.shape[1]): - out[inside, dim] = np.interp( - dist[inside], np.asarray(self.distance), coords[:, dim] + knots = np.asarray(self.distance, dtype=float) + out = {} + for name, values in self.coordinates.items(): + column = np.full(len(dist), np.nan) + column[inside] = np.interp( + dist[inside], knots, np.asarray(values, dtype=float) ) + out[name] = column return out @@ -985,7 +1038,19 @@ def _overlapping_epochs(items, key) -> list[tuple]: return out -_MIXED_DIMS_MSG = "Geometry segments mix coordinate dimensionalities {dims}." +def axis_columns(segment, crs) -> dict[str, int]: + """ + Return which of a segment's columns name which canonical axis. + + A column resolves to an axis when the CRS declares its name, or when it + is the canonical ``x``/``y``/``z`` alias of an axis the CRS has. Every + other column is a quantity along the fiber rather than a position. + """ + out = {} + for name in segment.coordinates: + with suppress(InvalidInventoryError): + out[name] = crs.axis_index(name) + return out def _track_identity_fields() -> Mapping[str, str]: @@ -1023,6 +1088,12 @@ def _track_identity_fields() -> Mapping[str, str]: ) +# The reserved names a geometry column may not take. The coordinate labels +# are left out of it: a column named for one is how a segment states that +# axis, and one the CRS does not declare is free to be a column of its own. +_RESERVED_COLUMN_NAMES = RESERVED_GROUP_NAMES - set(VALID_COORDINATE_LABELS) + + def _times_equal(time1, time2) -> bool: """Compare two epoch times, treating unset (NaT) times as equal.""" null1, null2 = np.isnat(time1), np.isnat(time2) @@ -1116,34 +1187,88 @@ def component_intervals(self) -> tuple[tuple[float, float], ...]: position = nxt return tuple(out) - def coordinates_at(self, distances) -> np.ndarray: + def coordinates_at(self, distances, crs) -> np.ndarray: """ Return CRS coordinates at the requested optical distances. - Uncovered distance returns NaN rows. Segment coverage is half-open, - with the end of each coverage run included: a distance on a - segment's last control point belongs to that segment unless another - segment claims it. + Only the columns which name a position axis are assembled; a segment + stating none of them contributes no position, however much else it + measures. Uncovered distance returns NaN rows. Segment coverage is + half-open, with the end of each coverage run included: a distance on + a segment's last control point belongs to that segment unless + another segment claims it. + + Parameters + ---------- + distances + The optical distances to place. + crs + The inventory's coordinate reference system, which is what + decides that a column is an axis and which axis it is. """ dist = np.atleast_1d(np.asarray(distances, dtype=float)) + out = np.full((len(dist), len(crs.coordinate_labels)), np.nan) if not self.geometry: - return np.full((len(dist), 1), np.nan) - dims = {len(seg.coordinates[0]) for seg in self.geometry} - if len(dims) > 1: - raise InvalidInventoryError(_MIXED_DIMS_MSG.format(dims=sorted(dims))) - out = np.full((len(dist), dims.pop()), np.nan) + return out masks = interval_masks(dist, [x.interval for x in self.geometry]) for segment, mask in zip(self.geometry, masks, strict=True): + axes = axis_columns(segment, crs) + if axes and len(set(axes.values())) != len(crs.coordinate_labels): + # Half a position is not one. A checked inventory cannot get + # here; an unchecked one says so rather than handing back a + # row whose missing axis reads as uncovered distance. + msg = ( + f"Geometry {segment.name!r} states the axes {sorted(axes)} " + f"but the CRS declares {list(crs.coordinate_labels)}; a " + "segment states every axis or none of them." + ) + raise InvalidInventoryError(msg) + if not axes or not np.any(mask): + continue + values = segment.interpolate(dist[mask]) + rows = np.flatnonzero(mask) + for name, index in axes.items(): + column = values[name] + # interpolate() reports its own coverage, which excludes the + # run end the mask includes, so fill that from the last point. + column[np.isnan(column)] = segment.coordinates[name][-1] + out[rows, index] = column + return out + + def column_at(self, name: str, distances) -> np.ndarray | None: + """ + Return one geometry column's values at the requested distances. + + None when no segment states the column. Coverage follows + `OpticalPath.coordinates_at`, and values never bridge two segments: + distance between them is uncovered, whatever either side holds. + """ + stating = [x for x in self.geometry if name in x.coordinates] + if not stating: + return None + dist = np.atleast_1d(np.asarray(distances, dtype=float)) + out = np.full(len(dist), np.nan) + masks = interval_masks(dist, [x.interval for x in stating]) + for segment, mask in zip(stating, masks, strict=True): if not np.any(mask): continue - # interpolate() reports its own coverage, which excludes the - # run end the mask includes, so fill that from the last point. - seg_coords = segment.interpolate(dist[mask]) - last = np.asarray(segment.coordinates, dtype=float)[-1] - seg_coords[np.isnan(seg_coords[:, 0])] = last - out[mask] = seg_coords + column = segment.interpolate(dist[mask])[name] + column[np.isnan(column)] = segment.coordinates[name][-1] + out[np.flatnonzero(mask)] = column return out + def column_units(self, name: str) -> str: + """Return the units the segments stating a column agree on.""" + stated = {x.units[name] for x in self.geometry if name in x.units} + return stated.pop() if len(stated) == 1 else "" + + def geometry_columns(self) -> tuple[str, ...]: + """Return every column name this path's geometry segments state.""" + seen: dict[str, None] = {} + for segment in self.geometry: + seen.update(dict.fromkeys(segment.coordinates)) + return tuple(seen) + def check(self, tolerance: float = 1e-9) -> Self: """ Check track rules for this path. @@ -1168,22 +1293,63 @@ def check(self, tolerance: float = 1e-9) -> Self: f"{name} interval ({lo}, {hi}) extends past path " f"span ({start}, {end})." ) - dims = {len(seg.coordinates[0]) for seg in self.geometry} - if len(dims) > 1: - errors.append(_MIXED_DIMS_MSG.format(dims=sorted(dims))) - for name, spans in (("geometry", geo_spans), ("coupling", coup_spans)): - overlap = _intervals_overlap(spans) - if overlap is not None: - errors.append( - f"Overlapping {name} intervals {overlap[0]} and " - f"{overlap[1]}; {name} is a function track." - ) + overlap = _intervals_overlap(coup_spans) + if overlap is not None: + errors.append( + f"Overlapping coupling intervals {overlap[0]} and " + f"{overlap[1]}; coupling is a function track." + ) + errors.extend(self._check_geometry_columns()) errors.extend(self._check_annotation_groups()) if errors: msg = "Optical path validation failed:\n" + "\n".join(errors) raise InvalidInventoryError(msg) return self + def _check_geometry_columns(self) -> list[str]: + """ + Check the geometry columns of this path against each other. + + Each column is its own function track, so two segments may overlap + as long as they do not state the same column over the same distance. + The rules which need the CRS -- which columns are axes -- are the + inventory's, since only it knows what the axes are. + """ + errors = [] + spans: dict[str, list[tuple[float, float]]] = {} + units: dict[str, set[str]] = {} + for segment in self.geometry: + for name in segment.coordinates: + spans.setdefault(name, []).append(segment.interval) + if name in segment.units: + units.setdefault(name, set()).add(segment.units[name]) + for name in sorted(set(spans) & _RESERVED_COLUMN_NAMES): + errors.append( + f"Geometry column {name!r} is a reserved name; a column " + "becomes a coordinate and cannot shadow a structural " + "coordinate or a typed track." + ) + groups = {x.group for x in self.annotations if x.group} + for name in sorted(set(spans) & groups): + errors.append( + f"{name!r} is both a geometry column and an annotation " + "group; one name is one coordinate." + ) + for name in sorted(spans): + overlap = _intervals_overlap(spans[name]) + if overlap is not None: + errors.append( + f"Overlapping geometry intervals {overlap[0]} and " + f"{overlap[1]} for column {name!r}; a column is a " + "function track." + ) + if len(stated := units.get(name, set())) > 1: + errors.append( + f"Geometry column {name!r} is stated in " + f"{sorted(stated)}; a column has one unit." + ) + return errors + def _check_annotation_groups(self) -> list[str]: """Check that each annotation group holds one kind of value.""" groups: dict[str, list] = {} @@ -1247,19 +1413,15 @@ def select(self, *, distance: tuple[float | None, float | None]) -> Self: dist = np.asarray(seg.distance, dtype=float) inside = (dist > new_lo) & (dist < new_hi) new_dist = np.concatenate([[new_lo], dist[inside], [new_hi]]) - coords = np.asarray(seg.coordinates, dtype=float) - new_coords = np.stack( - [ - np.interp(new_dist, dist, coords[:, dim]) - for dim in range(coords.shape[1]) - ], - axis=1, - ) + new_coords = { + name: tuple(np.interp(new_dist, dist, np.asarray(values, dtype=float))) + for name, values in seg.coordinates.items() + } geometry.append( seg.model_copy( update={ "distance": tuple(new_dist), - "coordinates": tuple(map(tuple, new_coords)), + "coordinates": new_coords, } ) ) @@ -1301,12 +1463,14 @@ def flip(d): geometry = [] for seg in self.geometry: dist = np.asarray(seg.distance, dtype=float) - coords = np.asarray(seg.coordinates, dtype=float) geometry.append( seg.model_copy( update={ "distance": tuple(flip(dist)[::-1]), - "coordinates": tuple(map(tuple, coords[::-1])), + "coordinates": { + name: tuple(values[::-1]) + for name, values in seg.coordinates.items() + }, } ) ) @@ -1998,8 +2162,17 @@ def _coord_names(self) -> tuple[str, ...]: groups: dict[str, None] = {} tracks: dict[str, dict[str, None]] = {} shapes: dict[str, set[str]] = {} + crs = self.coordinate_reference_system + columns: dict[str, None] = {} for path in self._optical_paths(): groups.update(dict.fromkeys(x.group for x in path.annotations if x.group)) + # The axes are listed above whatever any path states; what a + # path adds is the columns which are not positions. + for segment in path.geometry: + axes = axis_columns(segment, crs) + columns.update( + dict.fromkeys(x for x in segment.coordinates if x not in axes) + ) for track in TRACK_IDENTITY_FIELDS: for item in getattr(path, track): fields = tracks.setdefault(track, {}) @@ -2016,6 +2189,7 @@ def _coord_names(self) -> tuple[str, ...]: out[track] = None names = (f"{track}.{x}" for x in fields) out.update(dict.fromkeys(x for x in names if x not in unusable)) + out.update(columns) out.update(groups) return tuple(out) @@ -2050,8 +2224,7 @@ def check_width(width, what): for array in net.fiber_arrays: for path in array.optical_paths: for segment in path.geometry: - what = f"Geometry {segment.name!r}" - check_width(len(segment.coordinates[0]), what) + errors.extend(self._check_segment_axes(segment, crs)) for station in net.stations: if station.coordinates is not None: check_width(len(station.coordinates), f"Station {station.code!r}") @@ -2062,6 +2235,41 @@ def check_width(width, what): check_width(len(channel.coordinates), what) return errors + @staticmethod + def _check_segment_axes(segment, crs) -> list[str]: + """ + Check one geometry segment's columns against the CRS. + + A segment states every position axis or none of them: a partial + position is not one, and deciding what the missing axis meant is not + the reader's job. Which columns those are is the CRS's to say, which + is why this lives here rather than on the path. + """ + what = f"Geometry {segment.name!r}" if segment.name else "A geometry" + axes = axis_columns(segment, crs) + errors = [] + spellings: dict[int, list[str]] = {} + for name, index in axes.items(): + spellings.setdefault(index, []).append(name) + for index, names in sorted(spellings.items()): + if len(names) > 1: + errors.append( + f"{what} states axis {crs.coordinate_labels[index]!r} " + f"twice, as {sorted(names)}." + ) + if axes and len(set(axes.values())) != len(crs.coordinate_labels): + errors.append( + f"{what} states the axes {sorted(axes)} but the inventory " + f"CRS declares {list(crs.coordinate_labels)}; a segment " + "states every axis or none of them." + ) + if on_axes := sorted(set(segment.units) & set(axes)): + errors.append( + f"{what} states units for the axis column(s) {on_axes}; the " + "CRS states the units of its own axes." + ) + return errors + def resolve(self, acquisition_key: str, time=None) -> ResolvedContext: """ Resolve an acquisition_key (and time) to its inventory context. diff --git a/dascore/core/inventory_loader.py b/dascore/core/inventory_loader.py index 4c7a4e132..7c8a7a1bc 100644 --- a/dascore/core/inventory_loader.py +++ b/dascore/core/inventory_loader.py @@ -478,6 +478,9 @@ class _Table(NamedTuple): # field, so nothing else records where a row sits and it must place # each row unambiguously. places: bool = False + # The field every column but the order gathers into, keyed by header. + # None where each column names a field of the object directly. + columns: str | None = None # Keyed by CSV stem, which is the attribute the table fills. Held as a @@ -488,7 +491,9 @@ class _Table(NamedTuple): "optical_components": _Table(order="sequence", places=True), "coupling": _Table(), "annotations": _Table(), - "geometry": _Table(points=True, group="segment", order="distance"), + "geometry": _Table( + points=True, group="segment", order="distance", columns="coordinates" + ), "distance_map": _Table(points=True, order="distance"), } @@ -684,14 +689,18 @@ def _object_rows(frame: pd.DataFrame, table: _Table, path: Path) -> list[dict]: return out -def _point_rows(frame: pd.DataFrame, table: _Table, path: Path, axes) -> list[dict]: +def _point_rows( + frame: pd.DataFrame, table: _Table, path: Path, units: Mapping[str, str] +) -> list[dict]: """ Read a table whose every row is one control point. Points gather into objects holding parallel arrays: one object per value of the grouping column, or a single object when the attribute - is one. Coordinate columns are named by the CRS and are stored on the - canonical axes, so the frame decides which column is which. + is one. Where the table gathers its columns (``geometry``), every + column but the ordering one becomes an entry of that mapping keyed by + its header, so a name the file never states is a column the segment + does not have. """ _require_columns(frame, [table.order, table.group], path) _require_stated(frame, [table.order, table.group], path) @@ -707,8 +716,9 @@ def _point_rows(frame: pd.DataFrame, table: _Table, path: Path, axes) -> list[di out = [] for name, rows in groups: point: dict[str, Any] = {} if name is None else {"name": str(name)} + gathered: dict[str, tuple] = {} for column in rows.columns: - if column == table.group or column in axes: + if column == table.group: continue stated = rows[column].notna() if not stated.any(): @@ -733,35 +743,18 @@ def _point_rows(frame: pd.DataFrame, table: _Table, path: Path, axes) -> list[di "is stated by every point or by none." ) raise InvalidInventoryError(msg) - point[column] = tuple(rows[column]) - if axes: - point["coordinates"] = _coordinates(rows, axes, path) + if table.columns is not None and column != table.order: + gathered[column] = tuple(rows[column]) + else: + point[column] = tuple(rows[column]) + if table.columns is not None: + point[table.columns] = gathered + if stated := {x: units[x] for x in gathered if x in units}: + point["units"] = stated out.append(point) return out -def _coordinates(rows: pd.DataFrame, axes: Mapping[str, int], path: Path): - """ - Gather a geometry table's labelled columns onto the canonical axes. - - The CRS names the axes and states their order, so a header is read by - which axis it names rather than by where it sits in the file. - """ - ordered = sorted(axes, key=lambda label: axes[label]) - out = [] - for _, row in rows.iterrows(): - stated = [row[label] for label in ordered] - if any(pd.isnull(x) for x in stated): - missing = [x for x, v in zip(ordered, stated, strict=True) if pd.isnull(v)] - msg = ( - f"{_quote(path)} leaves {', '.join(missing)} empty for a point; " - "a coordinate states every axis its frame declares." - ) - raise InvalidInventoryError(msg) - out.append(tuple(stated)) - return tuple(out) - - def _is_path_dir(child: Path) -> bool: """Return True if a directory name claims to be an optical path epoch.""" # is_dir() follows a link, and the stray walk steps over one, so a @@ -970,35 +963,94 @@ def _load_table(path: Path, table: _Table, stem: str, crs): if frame.empty: msg = f"{_quote(path)} states no rows, so it describes no {stem}." raise InvalidInventoryError(msg) - axes = _geometry_axes(frame, crs, path) if stem == "geometry" else {} + units: Mapping[str, str] = {} + if stem == "geometry": + frame, units = _geometry_columns(frame, crs, path) if not table.points: rows = _object_rows(frame, table, path) if stem == "annotations": _parse_annotations(rows, path) return rows - built = _point_rows(frame, table, path, axes) + built = _point_rows(frame, table, path, units) # A single object rather than a collection: the table has no grouping # column because every point belongs to the one map it describes. return built if table.group is not None else built[0] -def _geometry_axes(frame: pd.DataFrame, crs, path: Path) -> dict[str, int]: +# A scalar column may state its units in its header, `depth (m)`. The axes +# take theirs from the CRS, so a unit on one of those is refused below. +_UNIT_SUFFIX = re.compile(r"^(?P.*?)\s*\((?P[^()]*)\)$") + + +def _geometry_columns(frame: pd.DataFrame, crs, path: Path): """ - Return which column names which canonical axis. + Read a geometry table's headers, and refuse what cannot be a column. - Coordinates are stored on the canonical axes while a geometry table - names them the way its frame does, so the CRS decides both which - headers are legal and what each one means. + A header naming an axis the CRS declares is that axis; every other one + is a numeric column in its own right, which may carry its units in + parentheses. The axes are all stated or none are: a partial position is + not a position, and guessing the missing axis is not a reader's job. """ labels = tuple(crs.coordinate_labels) - stated = {x for x in frame.columns} - {"segment", "distance"} - if stated != set(labels): + renamed, units = {}, {} + for header in frame.columns: + if header in {"segment", "distance"}: + continue + name, unit = header, "" + if (match := _UNIT_SUFFIX.match(header)) is not None: + name, unit = match.group("name"), match.group("units").strip() + renamed[header] = name + if not unit: + continue + if name in labels or name in {"x", "y", "z"}: + msg = ( + f"{_quote(path)} states units for {name!r}, which is a " + "position axis; the CRS states the units of its own axes." + ) + raise InvalidInventoryError(msg) + units[name] = unit + stated = set(renamed.values()) + if len(stated) != len(renamed): + repeated = sorted( + {x for x in renamed.values() if list(renamed.values()).count(x) > 1} + ) msg = ( - f"{_quote(path)} states the coordinate columns {sorted(stated)}, " - f"but its frame declares {list(labels)}." + f"{_quote(path)} names the column(s) {repeated} more than once; " + "one column states one thing." ) raise InvalidInventoryError(msg) - return {label: index for index, label in enumerate(labels)} + axes = stated & set(labels) + if axes and axes != set(labels): + msg = ( + f"{_quote(path)} states the axis column(s) {sorted(axes)}, but " + f"its frame declares {list(labels)}; a segment states every axis " + "or none of them." + ) + raise InvalidInventoryError(msg) + frame = frame.rename(columns=renamed) + return _numeric_columns(frame, sorted(stated), path), units + + +def _numeric_columns(frame: pd.DataFrame, columns, path: Path) -> pd.DataFrame: + """ + Read a geometry table's columns as numbers, refusing text. + + Text along distance is what annotations are for, and a column of it + here would otherwise reach the model as a string it cannot place. + """ + frame = frame.copy() + for column in columns: + values = pd.to_numeric(frame[column], errors="coerce") + if (bad := frame[column].notna() & values.isna()).any(): + first = frame.loc[bad, column].iloc[0] + msg = ( + f"{_quote(path)} states {first!r} in column {column!r}, " + "which is not a number. A geometry column is numeric; text " + "which varies along the fiber belongs in annotations.csv." + ) + raise InvalidInventoryError(msg) + frame[column] = values + return frame def _parse_annotations(rows: list[dict], path: Path) -> None: diff --git a/dascore/examples.py b/dascore/examples.py index 84be8f7bb..3b113c112 100644 --- a/dascore/examples.py +++ b/dascore/examples.py @@ -805,7 +805,13 @@ def inventory_patch_pair(): Geometry( name="trench", distance=(100.0, 400.0), - coordinates=((-117.0, 40.0, 1500.0), (-117.0, 40.1, 1500.0)), + # The canonical axis names, so the segment states the CRS's + # axes whatever this inventory's CRS happens to call them. + coordinates={ + "x": (-117.0, -117.0), + "y": (40.0, 40.1), + "z": (1500.0, 1500.0), + }, ), ), coupling=( diff --git a/dascore/proc/inventory.py b/dascore/proc/inventory.py index a122554b0..a6d38c0b6 100644 --- a/dascore/proc/inventory.py +++ b/dascore/proc/inventory.py @@ -25,7 +25,12 @@ validate_enrich_selection, ) from dascore.core.coords import BaseCoord, get_coord -from dascore.core.inventory import Interrogator, Inventory, ResolvedContext +from dascore.core.inventory import ( + Interrogator, + Inventory, + ResolvedContext, + axis_columns, +) from dascore.exceptions import ( InvalidInventoryError, ParameterError, @@ -321,13 +326,18 @@ def _get_blanket_coord_names(inventory, path) -> list[str]: """ Return the coordinate names a blanket request copies. - The geometry axes and the annotation groups: what the path says about + The geometry columns and the annotation groups: what the path says about each channel. Optical distance and the typed-track fields are asked for by name, since they restate what the patch's own axis and the inventory already record. """ - labels = inventory.coordinate_reference_system.coordinate_labels + crs = inventory.coordinate_reference_system + labels = crs.coordinate_labels out = ["x", "y", "z"][: len(labels)] if path.geometry else [] + # The axes are copied under their canonical names, so a column which is + # one of them is already covered; the rest come under their own. + axes = {x for segment in path.geometry for x in axis_columns(segment, crs)} + out += [x for x in path.geometry_columns() if x not in axes] seen = dict.fromkeys(x.group for x in path.annotations) return out + [x for x in seen if x] diff --git a/docs/recipes/tunnel_inventory.qmd b/docs/recipes/tunnel_inventory.qmd index 7d9e8c972..41155c4e3 100644 --- a/docs/recipes/tunnel_inventory.qmd +++ b/docs/recipes/tunnel_inventory.qmd @@ -217,7 +217,7 @@ print("borehole 3 sensing fiber:", at["borehole 3 down"][0], "to", at["borehole # Where the fiber is -Geometry places optical distance in space, one row per control point, grouped into segments. The columns after `segment` and `distance` are the axes the coordinate reference system declared, and naming any others is an error rather than an extra. +Geometry holds the curves measured along the fiber, one row per control point, grouped into segments. A column whose name the coordinate reference system declares is that position axis, and a segment stating one states them all; any other column is a number in its own right, which this deployment has no use for but a borehole survey would. The survey points are the ones lettered in @fig-tunnel. Each straight run of fiber goes from one to the next, so the table is those points paired with the component that runs between them. diff --git a/docs/tutorial/inventory.qmd b/docs/tutorial/inventory.qmd index f455823a3..4da6fae61 100644 --- a/docs/tutorial/inventory.qmd +++ b/docs/tutorial/inventory.qmd @@ -36,7 +36,7 @@ flowchart LR OpticalPath -->|annotations| OpticalPathAnnotation ``` -The components are the ordered physical pieces the light passes through, and they are what gives the path its length: each one tiles the interval after the last, so the path ends where the final component does. The other three describe intervals of that length — the geometry says where each distance is in space, the coupling says how the fiber is attached to the ground there, and the annotations name anything else worth recording per interval. Each of those covers what it covers, and none of them has to cover the whole path. +The components are the ordered physical pieces the light passes through, and they are what gives the path its length: each one tiles the interval after the last, so the path ends where the final component does. The other three describe intervals of that length — the geometry holds the measured curves along the fiber, of which the columns the coordinate reference system names are position; the coupling says how the fiber is attached to the ground there; and the annotations name anything else worth recording per interval. Each of those covers what it covers, and none of them has to cover the whole path. Objects reused in several places — interrogators, cables, enclosures, measurements — are written once under the inventory's `resources` and referred to elsewhere by their `resource_id`. Such a field accepts either the object itself or that string, which is what the dashed edges below mean. @@ -100,6 +100,18 @@ files = { "start_distance,end_distance,coupling_type,medium\n" "100,400,trench,soil\n" ), + # A geometry column is any number measured along the fiber. This one + # is the trench's own chainage, and the middle segment is a coil: the + # chainage stands still while ten meters of fiber goes by. + "fiber_arrays/DAS.R2D1/path/geometry.csv": ( + "segment,distance,chainage (m)\n" + "to the coil,100,1200\n" + "to the coil,240,1340\n" + "slack coil,240,1340\n" + "slack coil,250,1340\n" + "past the coil,250,1340\n" + "past the coil,400,1490\n" + ), "fiber_arrays/DAS.R2D1/path/annotations.csv": ( "start_distance,end_distance,group,value\n" "100,250,zone,north\n" @@ -124,6 +136,10 @@ A few rules make that directory readable without a schema in front of you: - **`path` is the one reserved *container* name.** A directory called `path` addresses an optical path rather than serializing an attribute named "path". (`attrs` and `inventory` are reserved as file stems — an entity's own object file, and the envelope — and the top-level directory names above are fixed.) - **Files that participate in no convention are ignored.** Photos, field notes, and deployment logs can live in the inventory directory undisturbed. +Two tracks can hold something that varies along the fiber, and which one to use is decided by what the values are. **Numbers go in `geometry.csv`**, one column per quantity, interpolated between the control points which state them; a column may name its units in its header, `chainage (m)`. **Text and vocabulary go in `annotations.csv`**, which is a set of intervals rather than a curve — a zone name has no meaning between two of them. A column of text in a geometry table is refused, and says so. + +Position is not a separate kind of thing: a column whose name the coordinate reference system declares — or its canonical `x`, `y`, `z` alias — *is* that axis, takes its units from the CRS, and must be stated with all its siblings, since half a position is not a position. A segment which names none of them, like the chainage above, states no position at all and is perfectly valid. The `segment` names themselves are for humans; what a channel gets is the columns. + Loading checks the document, so a directory that is not a valid inventory says so at the point it is read: ```{python} @@ -147,7 +163,7 @@ patch, example_inventory = inventory_patch_pair() # What an inventory can contribute -[`Inventory.get_names`](`dascore.core.inventory.Inventory.get_names`) lists the names an inventory could put on a patch, split by where each one lands. `attrs` are the observing-system facts, which are one value per patch; `coords` take a value per channel, because they describe a position along the fiber. +[`Inventory.get_names`](`dascore.core.inventory.Inventory.get_names`) lists the names an inventory could put on a patch, split by where each one lands. `attrs` are the observing-system facts, which are one value per patch; `coords` take a value per channel, because they describe somewhere along the fiber. ```{python} names = inventory.get_names() @@ -155,14 +171,14 @@ names = inventory.get_names() print("attrs:", names.attrs[:4], "...") print("coords:", [x for x in names.coords if "." not in x]) -# The annotation groups named in the CSV are among them. -assert {"zone", "noisy"} <= set(names.coords) +# The annotation groups and the geometry column named in the CSVs. +assert {"zone", "noisy", "chainage"} <= set(names.coords) assert "gauge_length" in names.attrs ``` -`zone` and `noisy` are there because the annotations CSV named them: an annotation group becomes a coordinate under its own name. `coupling` and `optical_components` are there for the same reason — a track contributes its name only where a path actually describes it. `distance` is where each channel sits on the optical path, and `x`/`y`/`z` and `longitude`/`latitude`/`elevation` are the two spellings of the axes the inventory's coordinate reference system declares, which a geometry track resolves to. +`zone` and `noisy` are there because the annotations CSV named them, and `chainage` because the geometry CSV did: a group and a column each become a coordinate under their own name. `coupling`, `geometry`, and `optical_components` are there for the same reason — a track contributes its name only where a path actually describes it. `distance` is where each channel sits on the optical path, and `x`/`y`/`z` and `longitude`/`latitude`/`elevation` are the two spellings of the axes the inventory's coordinate reference system declares. -Listing a name is not promising a value for it. This example states no geometry, so the spatial names resolve to nothing until it gets one. +Listing a name is not promising a value for it. This example's geometry states chainage and no position, so the spatial names resolve to nothing until some segment states them. # Attaching an inventory to a spool diff --git a/tests/test_core/test_inventory.py b/tests/test_core/test_inventory.py index bdc6378a5..ae75e9e78 100644 --- a/tests/test_core/test_inventory.py +++ b/tests/test_core/test_inventory.py @@ -36,7 +36,7 @@ def build_inventory() -> inv.Inventory: geometry = inv.Geometry( name="survey", distance=(0.0, 100.0, 200.0), - coordinates=((0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (2.0, 0.0, 1.0)), + coordinates={"x": (0.0, 1.0, 2.0), "y": (0.0, 0.0, 0.0), "z": (0.0, 0.0, 1.0)}, ) path = inv.OpticalPath( name="main", @@ -96,7 +96,11 @@ def build_full_inventory() -> inv.Inventory: inv.Geometry( name="trench", distance=(100.0, 400.0), - coordinates=((-117.0, 40.0, 1500.0), (-117.0, 40.1, 1500.0)), + coordinates={ + "x": (-117.0, -117.0), + "y": (40.0, 40.1), + "z": (1500.0, 1500.0), + }, ), ), coupling=( @@ -178,38 +182,205 @@ def build_full_inventory() -> inv.Inventory: ).check() +class TestGeometryColumns: + """A geometry states named numeric columns, of which some are axes.""" + + @staticmethod + def _inventory(*geometry, crs=None, annotations=()): + """Wrap geometry segments in the smallest inventory holding them.""" + path = inv.OpticalPath( + optical_components=(inv.FiberSegment(optical_length=1000.0),), + geometry=geometry, + annotations=annotations, + ) + array = inv.FiberArray(code="L001", optical_paths=(path,)) + return inv.Inventory( + networks=(inv.Network(code="XX", fiber_arrays=(array,)),), + **({"coordinate_reference_system": crs} if crs else {}), + ) + + def test_a_column_which_is_not_a_position(self): + """Chainage is a curve along the fiber and no part of a position.""" + chainage = inv.Geometry( + name="chainage", + distance=(0.0, 100.0), + coordinates={"chainage": (1200.0, 1300.0)}, + units={"chainage": "m"}, + ) + assert self._inventory(chainage).check() is not None + out = chainage.interpolate([50.0]) + assert out["chainage"][0] == 1250.0 + + def test_a_segment_with_no_axes_is_legal(self): + """A path may be described without ever being placed in space.""" + depth = inv.Geometry( + distance=(0.0, 40.0), coordinates={"borehole_depth": (0.0, 40.0)} + ) + inventory = self._inventory(depth) + assert inventory.check() is inventory + crs = inventory.coordinate_reference_system + path = inventory.networks[0].fiber_arrays[0].optical_paths[0] + assert np.all(np.isnan(path.coordinates_at([20.0], crs))) + + def test_axes_are_all_or_none(self): + """Half a position is not a position.""" + partial = inv.Geometry( + name="partial", distance=(0.0, 10.0), coordinates={"x": (0.0, 1.0)} + ) + with pytest.raises(InvalidInventoryError, match="every axis or none"): + self._inventory(partial).check() + + def test_an_axis_stated_twice(self): + """`x` and the label the CRS gives it are one axis, not two.""" + crs = inv.CoordinateReferenceSystem( + coordinate_labels=("easting", "northing"), units=("meter", "meter") + ) + doubled = inv.Geometry( + name="doubled", + distance=(0.0, 10.0), + coordinates={ + "x": (0.0, 1.0), + "easting": (0.0, 1.0), + "northing": (0.0, 1.0), + }, + ) + with pytest.raises(InvalidInventoryError, match="twice"): + self._inventory(doubled, crs=crs).check() + + def test_overlap_is_refused_per_column(self): + """Two segments may overlap unless they state the same column.""" + first = inv.Geometry(distance=(0.0, 60.0), coordinates={"depth": (0.0, 6.0)}) + second = inv.Geometry(distance=(50.0, 80.0), coordinates={"depth": (5.0, 8.0)}) + with pytest.raises(InvalidInventoryError, match="for column 'depth'"): + self._inventory(first, second).check() + + def test_different_columns_may_overlap(self): + """Each column is its own function track, so they are independent.""" + depth = inv.Geometry(distance=(0.0, 60.0), coordinates={"depth": (0.0, 6.0)}) + azimuth = inv.Geometry( + distance=(50.0, 80.0), coordinates={"azimuth": (5.0, 8.0)} + ) + inventory = self._inventory(depth, azimuth) + assert inventory.check() is inventory + + def test_a_reserved_column_name(self): + """A column becomes a coordinate, so it cannot shadow one.""" + clash = inv.Geometry(distance=(0.0, 10.0), coordinates={"time": (0.0, 1.0)}) + with pytest.raises(InvalidInventoryError, match="reserved name"): + self._inventory(clash).check() + + def test_a_column_which_is_also_an_annotation_group(self): + """One name is one coordinate, whichever track would define it.""" + column = inv.Geometry(distance=(0.0, 10.0), coordinates={"zone": (0.0, 1.0)}) + annotation = inv.OpticalPathAnnotation( + start_distance=0.0, end_distance=10.0, group="zone", value=1.0 + ) + with pytest.raises(InvalidInventoryError, match="one name is one coordinate"): + self._inventory(column, annotations=(annotation,)).check() + + def test_units_on_an_axis_are_refused(self): + """The CRS states the units of its own axes.""" + segment = inv.Geometry( + name="axed", + distance=(0.0, 10.0), + coordinates={"x": (0.0, 1.0), "y": (0.0, 1.0), "z": (0.0, 1.0)}, + units={"x": "furlong"}, + ) + with pytest.raises(InvalidInventoryError, match="states the units of its own"): + self._inventory(segment).check() + + def test_units_for_a_column_which_is_not_there(self): + """A unit sitting on nothing is a typo no reader would find.""" + with pytest.raises(ValidationError, match="no column for"): + inv.Geometry( + distance=(0.0, 10.0), + coordinates={"depth": (0.0, 1.0)}, + units={"dpeth": "m"}, + ) + + def test_one_unit_per_column(self): + """Two segments cannot measure one column in two units.""" + meters = inv.Geometry( + distance=(0.0, 10.0), + coordinates={"depth": (0.0, 1.0)}, + units={"depth": "m"}, + ) + feet = inv.Geometry( + distance=(20.0, 30.0), + coordinates={"depth": (0.0, 1.0)}, + units={"depth": "ft"}, + ) + with pytest.raises(InvalidInventoryError, match="has one unit"): + self._inventory(meters, feet).check() + + def test_a_crs_label_the_crs_does_not_declare(self): + """`depth` is a column of its own where the CRS spends z on elevation.""" + crs = inv.CoordinateReferenceSystem( + coordinate_labels=("easting", "northing", "elevation"), + units=("meter", "meter", "meter"), + ) + segment = inv.Geometry( + distance=(0.0, 40.0), + coordinates={"depth": (0.0, 40.0)}, + units={"depth": "m"}, + ) + inventory = self._inventory(segment, crs=crs) + assert inventory.check() is inventory + assert "depth" in inventory.get_names().coords + + def test_a_column_never_bridges_two_segments(self): + """Distance between two segments is uncovered, whatever they hold.""" + first = inv.Geometry(distance=(0.0, 10.0), coordinates={"depth": (0.0, 1.0)}) + second = inv.Geometry(distance=(20.0, 30.0), coordinates={"depth": (2.0, 3.0)}) + path = self._inventory(first, second).networks[0].fiber_arrays[0] + values = path.optical_paths[0].column_at("depth", [5.0, 15.0, 25.0]) + assert not np.isnan(values[0]) and not np.isnan(values[2]) + assert np.isnan(values[1]) + + def test_a_column_no_segment_states(self): + """None, so the caller's on_missing policy rules rather than a nan.""" + segment = inv.Geometry(distance=(0.0, 10.0), coordinates={"depth": (0.0, 1.0)}) + path = self._inventory(segment).networks[0].fiber_arrays[0].optical_paths[0] + assert path.column_at("azimuth", [5.0]) is None + + class TestGeometry: """Geometry segment rules.""" def test_requires_two_points(self): """Requires two points.""" with pytest.raises(ValidationError, match="at least 2 control points"): - inv.Geometry(distance=(1.0,), coordinates=((0.0, 0.0),)) + inv.Geometry(distance=(1.0,), coordinates={"x": (0.0,), "y": (0.0,)}) def test_strictly_increasing(self): """Strictly increasing.""" with pytest.raises(ValidationError, match="strictly increasing"): - inv.Geometry(distance=(1.0, 1.0), coordinates=((0.0, 0.0), (1.0, 1.0))) + inv.Geometry( + distance=(1.0, 1.0), coordinates={"x": (0.0, 1.0), "y": (0.0, 1.0)} + ) def test_paired_lengths(self): """Paired lengths.""" - with pytest.raises(ValidationError, match="same length"): - inv.Geometry(distance=(0.0, 1.0), coordinates=((0.0, 0.0),)) + with pytest.raises(ValidationError, match="one value per distance"): + inv.Geometry(distance=(0.0, 1.0), coordinates={"x": (0.0,), "y": (0.0,)}) def test_coil_repeated_coordinates(self): """A coil interpolates to a constant coordinate.""" coil = inv.Geometry( distance=(1200.0, 1300.0), - coordinates=((500.0, 120.0), (500.0, 120.0)), + coordinates={"x": (500.0, 500.0), "y": (120.0, 120.0)}, ) out = coil.interpolate([1200.0, 1250.0, 1299.0]) - assert np.allclose(out, [[500.0, 120.0]] * 3) + assert np.allclose(out["x"], [500.0] * 3) + assert np.allclose(out["y"], [120.0] * 3) def test_uncovered_is_nan(self): """Uncovered is nan.""" - geo = inv.Geometry(distance=(10.0, 20.0), coordinates=((0.0, 0.0), (1.0, 1.0))) + geo = inv.Geometry( + distance=(10.0, 20.0), coordinates={"x": (0.0, 1.0), "y": (0.0, 1.0)} + ) out = geo.interpolate([5.0, 25.0]) - assert np.all(np.isnan(out)) + assert all(np.all(np.isnan(x)) for x in out.values()) class TestPathTracks: @@ -238,7 +409,7 @@ def test_coupling_overlap_raises(self): def test_geometry_overlap_raises(self): """Geometry overlap raises.""" - seg = dict(coordinates=((0.0, 0.0), (1.0, 1.0))) + seg = dict(coordinates={"x": (0.0, 1.0), "y": (0.0, 1.0)}) path = inv.OpticalPath( optical_components=(inv.FiberSegment(optical_length=100.0),), geometry=( @@ -341,14 +512,18 @@ def test_out_of_bounds_raises(self): def test_outer_endpoint_included(self): """The outermost covered endpoint of the geometry track resolves.""" - path = build_inventory().networks[0].fiber_arrays[0].optical_paths[0] - coords = path.coordinates_at([200.0]) + inventory = build_inventory() + path = inventory.networks[0].fiber_arrays[0].optical_paths[0] + crs = inventory.coordinate_reference_system + coords = path.coordinates_at([200.0], crs) assert np.allclose(coords, [[2.0, 0.0, 1.0]]) def test_uncovered_distance_is_nan(self): """Uncovered distance is nan.""" - path = build_inventory().networks[0].fiber_arrays[0].optical_paths[0] - assert np.all(np.isnan(path.coordinates_at([225.0]))) + inventory = build_inventory() + path = inventory.networks[0].fiber_arrays[0].optical_paths[0] + crs = inventory.coordinate_reference_system + assert np.all(np.isnan(path.coordinates_at([225.0], crs))) class TestDistanceMap: @@ -754,10 +929,10 @@ def test_duplicate_resource_ids_raise(self): resources=[inv.Cable(resource_id="x"), inv.Cable(resource_id="x")] ) - def test_zero_dim_geometry_raises(self): - """Zero dim geometry raises.""" - with pytest.raises(ValidationError, match="nonzero"): - inv.Geometry(distance=(0.0, 1.0), coordinates=((), ())) + def test_columnless_geometry_raises(self): + """A segment which measures nothing describes nothing.""" + with pytest.raises(ValidationError, match="states no columns"): + inv.Geometry(distance=(0.0, 1.0), coordinates={}) def test_station_extra_fields_forbidden(self): """Coordinates are canonical (x, y, z); label fields are not stored.""" @@ -1007,22 +1182,26 @@ def test_duplicate_array_codes_raise(self): with pytest.raises(InvalidInventoryError, match="Duplicate fiber array"): net.check() - def test_mixed_dimensionality_geometry_raises(self): - """Segments with different coordinate dims fail the path check.""" + def test_partial_axes_fail_the_inventory_check(self): + """A segment stating some axes and not others fails the check.""" path = inv.OpticalPath( optical_components=(inv.FiberSegment(optical_length=100.0),), geometry=( inv.Geometry( - distance=(0.0, 10.0), coordinates=((0.0, 0.0), (1.0, 1.0)) + distance=(0.0, 10.0), coordinates={"x": (0.0, 1.0), "y": (0.0, 1.0)} ), inv.Geometry( distance=(20.0, 30.0), - coordinates=((0.0, 0.0, 0.0), (1.0, 1.0, 1.0)), + coordinates={"x": (0.0, 1.0), "y": (0.0, 1.0), "z": (0.0, 1.0)}, ), ), ) - with pytest.raises(InvalidInventoryError, match="dimensionalities"): - path.check() + array = inv.FiberArray(code="L001", optical_paths=(path,)) + inventory = inv.Inventory( + networks=(inv.Network(code="XX", fiber_arrays=(array,)),) + ) + with pytest.raises(InvalidInventoryError, match="every axis or none"): + inventory.check() def test_inventory_function_dispatch(self): """dc.inventory handles bad input with clear errors.""" @@ -1076,7 +1255,9 @@ def test_nonfinite_interval_values_raise(self): start_distance=np.nan, end_distance=10.0, coupling_type="trench" ) with pytest.raises(ValidationError, match="finite"): - inv.Geometry(distance=(0.0, np.inf), coordinates=((0.0, 0.0), (1.0, 1.0))) + inv.Geometry( + distance=(0.0, np.inf), coordinates={"x": (0.0, 1.0), "y": (0.0, 1.0)} + ) with pytest.raises(ValidationError, match="finite"): inv.DistanceMap(channel=(0.0, np.inf), distance=(0.0, 1.0)) @@ -1146,11 +1327,12 @@ def test_coordinates_at_without_geometry(self): path = inv.OpticalPath( optical_components=(inv.FiberSegment(optical_length=10.0),) ) - assert np.all(np.isnan(path.coordinates_at([5.0]))) + crs = inv.CoordinateReferenceSystem() + assert np.all(np.isnan(path.coordinates_at([5.0], crs))) def test_select_drops_out_of_range_geometry(self): """Selection drops segments entirely outside the clip.""" - seg = dict(coordinates=((0.0, 0.0), (1.0, 1.0))) + seg = dict(coordinates={"x": (0.0, 1.0), "y": (0.0, 1.0)}) path = inv.OpticalPath( optical_components=(inv.FiberSegment(optical_length=100.0),), geometry=( @@ -1776,16 +1958,17 @@ def test_coordinates_at_rejects_mixed_dimensions(self): optical_components=(inv.FiberSegment(optical_length=100.0),), geometry=( inv.Geometry( - distance=(0.0, 10.0), coordinates=((0.0, 0.0), (1.0, 1.0)) + distance=(0.0, 10.0), coordinates={"x": (0.0, 1.0), "y": (0.0, 1.0)} ), inv.Geometry( distance=(20.0, 30.0), - coordinates=((0.0, 0.0, 0.0), (1.0, 1.0, 1.0)), + coordinates={"x": (0.0, 1.0), "y": (0.0, 1.0), "z": (0.0, 1.0)}, ), ), ) - with pytest.raises(InvalidInventoryError, match="mix coordinate"): - path.coordinates_at([5.0]) + crs = inv.CoordinateReferenceSystem() + with pytest.raises(InvalidInventoryError, match="every axis or none"): + path.coordinates_at([5.0], crs) def test_coordinate_width_must_match_crs(self): """Coordinates are read through the CRS, so they must fit its axes.""" @@ -1795,7 +1978,7 @@ def test_coordinate_width_must_match_crs(self): inv.Geometry( name="flat", distance=(0.0, 10.0), - coordinates=((0.0, 0.0), (1.0, 1.0)), + coordinates={"x": (0.0, 1.0), "y": (0.0, 1.0)}, ), ), ) @@ -1803,7 +1986,7 @@ def test_coordinate_width_must_match_crs(self): inventory = inv.Inventory( networks=(inv.Network(code="XX", fiber_arrays=(array,)),) ) - with pytest.raises(InvalidInventoryError, match="CRS declares 3 axes"): + with pytest.raises(InvalidInventoryError, match="every axis or none"): inventory.check() def test_station_coordinate_width_must_match_crs(self): @@ -1845,8 +2028,10 @@ def test_dip_range(self): def test_geometry_coordinates_must_be_finite(self): """A nan control point would read as uncovered distance.""" - with pytest.raises(ValidationError, match="must be finite"): - inv.Geometry(distance=(0.0, 1.0), coordinates=((np.nan, 0.0), (1.0, 1.0))) + with pytest.raises(ValidationError, match="must hold finite"): + inv.Geometry( + distance=(0.0, 1.0), coordinates={"x": (np.nan, 1.0), "y": (0.0, 1.0)} + ) def test_point_coordinates_must_be_finite(self): """Stations and channels name real positions.""" diff --git a/tests/test_core/test_inventory_loader.py b/tests/test_core/test_inventory_loader.py index 087228a87..030d6d58e 100644 --- a/tests/test_core/test_inventory_loader.py +++ b/tests/test_core/test_inventory_loader.py @@ -1166,13 +1166,81 @@ def test_an_empty_cell_is_unset(self, make_inventory): assert path.coupling[0].description == "" assert path.coupling[1].description == "backfilled" + def test_a_column_which_is_not_an_axis(self, make_inventory): + """A header the CRS does not name is a curve along the fiber.""" + files = { + **MINIMAL, + **TRACKS, + "fiber_arrays/DAS.L001/path/geometry.csv": ( + "segment,distance,longitude,latitude,elevation,borehole_depth\n" + "S100,100.0,-117.0,40.0,687.0,0.0\n" + "S100,102.0,-117.1,40.1,685.0,2.0\n" + ), + } + geometry = one_path(make_inventory(files)).geometry[0] + assert geometry.coordinates["borehole_depth"] == (0.0, 2.0) + assert geometry.units == {} + + def test_a_column_states_its_units_in_its_header(self, make_inventory): + """`depth (m)` is the column `depth`, measured in meters.""" + files = { + **MINIMAL, + **TRACKS, + "fiber_arrays/DAS.L001/path/geometry.csv": ( + "segment,distance,depth (m)\nS100,100.0,0.0\nS100,102.0,2.0\n" + ), + } + geometry = one_path(make_inventory(files)).geometry[0] + assert geometry.coordinates["depth"] == (0.0, 2.0) + assert geometry.units["depth"] == "m" + + def test_units_on_an_axis_header(self, make_inventory): + """The CRS states the units of its axes, so a header may not.""" + files = { + **MINIMAL, + **TRACKS, + "fiber_arrays/DAS.L001/path/geometry.csv": ( + "segment,distance,longitude,latitude,elevation (ft)\n" + "S100,100.0,-117.0,40.0,687.0\n" + "S100,102.0,-117.1,40.1,685.0\n" + ), + } + with pytest.raises(InvalidInventoryError, match="units of its own axes"): + make_inventory(files) + + def test_a_column_of_text_is_refused(self, make_inventory): + """Text along distance is what annotations are for.""" + files = { + **MINIMAL, + **TRACKS, + "fiber_arrays/DAS.L001/path/geometry.csv": ( + "segment,distance,zone\nS100,100.0,north\nS100,102.0,south\n" + ), + } + with pytest.raises(InvalidInventoryError, match=r"annotations\.csv"): + make_inventory(files) + + def test_one_column_stated_twice(self, make_inventory): + """A unit suffix is not a second column, however it is spelled.""" + files = { + **MINIMAL, + **TRACKS, + "fiber_arrays/DAS.L001/path/geometry.csv": ( + "segment,distance,depth,depth (m)\nS100,100.0,0.0,0.0\n" + "S100,102.0,2.0,2.0\n" + ), + } + with pytest.raises(InvalidInventoryError, match="more than once"): + make_inventory(files) + def test_a_geometry_names_its_axes_the_way_its_frame_does(self, make_inventory): """Coordinates are stored canonically, whatever the CRS calls them.""" path = one_path(make_inventory({**MINIMAL, **TRACKS})) geometry = path.geometry[0] assert geometry.distance == (100.0, 102.0) - # longitude, latitude, elevation are axes 0, 1, 2 of the default CRS. - assert geometry.coordinates[0] == (-117.0, 40.0, 687.0) + # A column is keyed by the header which states it. + assert geometry.coordinates["longitude"] == (-117.0, -117.1) + assert geometry.coordinates["elevation"] == (687.0, 685.0) def test_a_declared_frame_renames_the_axes(self, make_inventory): """The envelope decides which headers a geometry table may state.""" @@ -1194,18 +1262,20 @@ def test_a_declared_frame_renames_the_axes(self, make_inventory): ), } geometry = one_path(make_inventory(files)).geometry[0] - assert geometry.coordinates[0] == (2562048.25, 1137365.53, 687.0) + assert geometry.coordinates["x"] == (2562048.25, 2562048.17) + assert geometry.coordinates["elevation"] == (687.0, 685.0) - def test_a_header_the_frame_does_not_declare(self, make_inventory): - """A geometry header disagreeing with the frame raises.""" + def test_a_partial_set_of_axes(self, make_inventory): + """A segment states every axis the frame declares, or none of them.""" files = { **MINIMAL, **TRACKS, "fiber_arrays/DAS.L001/path/geometry.csv": ( - "segment,distance,x,y,z\nS100,100.0,1.0,2.0,3.0\n" + "segment,distance,longitude,latitude\nS100,100.0,1.0,2.0\n" + "S100,102.0,1.1,2.1\n" ), } - with pytest.raises(InvalidInventoryError, match="its frame declares"): + with pytest.raises(InvalidInventoryError, match="every axis or none"): make_inventory(files) def test_a_point_missing_an_axis(self, make_inventory): @@ -1807,7 +1877,7 @@ def test_points_gather_into_the_segment_which_names_them(self, make_inventory): # Interleaved in the file; gathered by name and ordered by distance. assert geometry["S100"].distance == (100.0, 102.0) assert geometry["S120"].distance == (200.0, 202.0) - assert geometry["S100"].coordinates[1] == (-117.1, 40.1, 685.0) + assert geometry["S100"].coordinates["latitude"] == (40.0, 40.1) def test_a_wrong_cell_names_the_field_it_could_not_take(self, make_inventory): """Matching only 'Could not read' would pass for any path failure.""" diff --git a/tests/test_proc/test_proc_inventory.py b/tests/test_proc/test_proc_inventory.py index 24dbd9b3d..fdfbfe44c 100644 --- a/tests/test_proc/test_proc_inventory.py +++ b/tests/test_proc/test_proc_inventory.py @@ -322,6 +322,81 @@ def test_distance_collision_raises(self, patch, inventory): patch.enrich(inventory, attrs=False, coords=("distance",)) +class TestGeometryColumns: + """A geometry column which is not a position still reaches the patch.""" + + @staticmethod + def _with_depth(inventory, **kwargs): + """Add a borehole-depth column over part of the path.""" + segment = Geometry( + name="hole", + distance=(100.0, 200.0), + coordinates={"borehole_depth": (0.0, 100.0)}, + units={"borehole_depth": "m"}, + **kwargs, + ) + path = inventory.networks[0].fiber_arrays[0].optical_paths[0] + return _replace_path(inventory, geometry=(*path.geometry, segment)) + + def test_a_column_becomes_a_coordinate(self, patch, inventory): + """With the units the segment states for it.""" + inv = self._with_depth(inventory) + out = patch.enrich(inv, attrs=False, coords=("borehole_depth",)) + coord = out.get_coord("borehole_depth") + # Channel 0 is path distance 100, the top of the hole. + assert coord.values[0] == 0.0 + assert coord.units is not None and "m" in str(coord.units) + + def test_uncovered_channels_are_nan(self, patch, inventory): + """The column covers 100 to 200; the path runs to 400.""" + inv = self._with_depth(inventory) + out = patch.enrich(inv, attrs=False, coords=("borehole_depth",)) + values = out.get_coord("borehole_depth").values + assert not np.isnan(values[0]) + assert np.isnan(values[-1]) + + def test_a_blanket_request_includes_it(self, patch, inventory): + """It is one of the things the path says about a channel.""" + inv = self._with_depth(inventory) + out = patch.enrich(inv, attrs=False) + assert "borehole_depth" in set(out.coords.coord_map) + + def test_get_names_lists_it(self, inventory): + """So a reader can find it without opening the CSV.""" + inv = self._with_depth(inventory) + assert "borehole_depth" in inv.get_names().coords + + def test_selection_trims_channels(self, patch, inventory): + """Selecting on a column is selecting on the fiber it describes.""" + inv = self._with_depth(inventory) + spool = dc.spool(patch).attach_inventory(inv) + out = spool.select(borehole_depth=(0, 50))[0] + depth = out.get_coord("distance") + # 100 to 150 m of path is the upper half of the hole, and the + # acquisition maps path distance 100 onto channel 0. + assert depth.max() < patch.get_coord("distance").max() + + def test_a_column_does_not_bridge_segments(self, patch, inventory): + """Two holes are two holes, and the fiber between them is neither.""" + first = Geometry( + name="hole 1", + distance=(100.0, 150.0), + coordinates={"borehole_depth": (0.0, 50.0)}, + ) + second = Geometry( + name="hole 2", + distance=(300.0, 350.0), + coordinates={"borehole_depth": (0.0, 50.0)}, + ) + path = inventory.networks[0].fiber_arrays[0].optical_paths[0] + inv = _replace_path(inventory, geometry=(*path.geometry, first, second)) + out = patch.enrich(inv, attrs=False, coords=("borehole_depth",)) + values = out.get_coord("borehole_depth").values + assert not np.isnan(values[10]) # in the first hole + assert not np.isnan(values[200]) # channel 200 is the top of the second + assert np.isnan(values[100]) # between them, and in neither + + class TestCoords: """What the optical path projects onto the patch.""" @@ -536,11 +611,18 @@ def test_empty_track_is_missing(self, patch, inventory): patch.enrich(inv, attrs=False, coords=("coupling.medium",)) def test_axis_missing_from_geometry(self, patch, inventory): - """A two-dimensional geometry has no third axis to return.""" + """A two-axis CRS has no third axis to return.""" flat = Geometry( - name="flat", distance=(100.0, 400.0), coordinates=((0.0, 0.0), (1.0, 1.0)) + name="flat", + distance=(100.0, 400.0), + coordinates={"x": (0.0, 1.0), "y": (0.0, 1.0)}, + ) + crs = CoordinateReferenceSystem( + coordinate_labels=("easting", "northing"), units=("meter", "meter") + ) + inv = _replace_path( + inventory.new(coordinate_reference_system=crs), geometry=(flat,) ) - inv = _replace_path(inventory, geometry=(flat,)) with pytest.raises(PatchError, match="defines no 'z'"): patch.enrich(inv, attrs=False, coords=("z",)) @@ -680,8 +762,12 @@ def test_endpoint_belongs_to_its_own_run(self, patch, inventory): def test_geometry_endpoint_is_local(self, patch, inventory): """The same rule holds for the geometry track.""" - first = Geometry(distance=(100.0, 200.0), coordinates=((0.0, 0.0), (1.0, 1.0))) - second = Geometry(distance=(300.0, 400.0), coordinates=((3.0, 3.0), (4.0, 4.0))) + columns = {"x": (0.0, 1.0), "y": (0.0, 1.0), "z": (0.0, 1.0)} + first = Geometry(distance=(100.0, 200.0), coordinates=columns) + second = Geometry( + distance=(300.0, 400.0), + coordinates={"x": (3.0, 4.0), "y": (3.0, 4.0), "z": (3.0, 4.0)}, + ) inv = _replace_path(inventory, geometry=(first, second)) out = patch.enrich(inv, attrs=False, coords=("x",)) # channel 100 is path distance 200, the last point of the first segment From c2a728c7a54bef7de23cbef2fd9c7eca857e3857 Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Sat, 15 Aug 2026 18:04:01 +0200 Subject: [PATCH 2/5] Take a review of the column form Eight findings, of which two placed the fiber wrongly rather than loudly. A segment stating `x`, `y`, `z` and another stating `longitude`, `latitude`, `elevation` over the same distance are two spellings of one axis, and checking columns by name alone let them both through; the axes are now checked again against what the CRS says they are. A segment stating `x` *and* `longitude` has three distinct axes and passed the partial-position guard, after which whichever spelling the mapping held last won -- both are now refused where the position is assembled, not only where it is checked. The rest: a segment which measures without placing no longer claims the position track's run end, or offers an axis of nan where nothing places the fiber at all; a dotted column name is refused, since that is how a field of a typed track is asked for; select and reverse revalidate rather than copying past the validators, so the columns of a selected segment stay frozen and checked; `distance (m)` beside `distance` is a duplicate rather than a pandas TypeError; and `z` is a column of its own under a CRS which declares only two axes, in a CSV as it already was in the model. --- dascore/core/_spool_inventory.py | 5 +- dascore/core/inventory.py | 132 +++++++++++++++-------- dascore/core/inventory_loader.py | 26 +++-- dascore/proc/inventory.py | 6 +- tests/test_core/test_inventory.py | 102 ++++++++++++++++++ tests/test_core/test_inventory_loader.py | 31 ++++++ tests/test_proc/test_proc_inventory.py | 15 +++ 7 files changed, 259 insertions(+), 58 deletions(-) diff --git a/dascore/core/_spool_inventory.py b/dascore/core/_spool_inventory.py index 629eac569..55cbc2e16 100644 --- a/dascore/core/_spool_inventory.py +++ b/dascore/core/_spool_inventory.py @@ -34,6 +34,7 @@ Inventory, ResolvedContext, _annotation_kind, + axis_columns, interval_masks, ) from dascore.core.inventory_loader import BLESSED_NAME, find_inventory @@ -666,7 +667,9 @@ def _get_geometry_coord(inventory, path, label, distances): index = crs.axis_index(label) except InvalidInventoryError: return None - if not path.geometry: + # Geometry which places nothing defines no axis, so the caller's + # on_missing policy rules rather than a column of nan. + if not any(axis_columns(x, crs) for x in path.geometry): return None coords = path.coordinates_at(distances, crs) if index >= coords.shape[1]: diff --git a/dascore/core/inventory.py b/dascore/core/inventory.py index e2ef9dc86..b82a6b873 100644 --- a/dascore/core/inventory.py +++ b/dascore/core/inventory.py @@ -1053,6 +1053,39 @@ def axis_columns(segment, crs) -> dict[str, int]: return out +def _axis_set_errors(segment, axes: Mapping[str, int], crs) -> list[str]: + """Return what is wrong with the set of axes one segment states.""" + what = f"Geometry {segment.name!r}" if segment.name else "A geometry" + spellings: dict[int, list[str]] = {} + for name, index in axes.items(): + spellings.setdefault(index, []).append(name) + errors = [ + f"{what} states axis {crs.coordinate_labels[index]!r} twice, as " + f"{sorted(names)}." + for index, names in sorted(spellings.items()) + if len(names) > 1 + ] + if axes and len(spellings) != len(crs.coordinate_labels): + errors.append( + f"{what} states the axes {sorted(axes)} but the inventory CRS " + f"declares {list(crs.coordinate_labels)}; a segment states every " + "axis or none of them." + ) + return errors + + +def _check_axis_set(segment, axes: Mapping[str, int], crs) -> None: + """ + Refuse a segment whose axes are partial or doubly spelled. + + A checked inventory cannot reach this; an unchecked one says so rather + than filling a position from whichever spelling the mapping happened to + hold last, or leaving an axis reading as uncovered distance. + """ + if errors := _axis_set_errors(segment, axes, crs): + raise InvalidInventoryError(" ".join(errors)) + + def _track_identity_fields() -> Mapping[str, str]: """ Map each typed track of an optical path to the field its name means. @@ -1208,22 +1241,18 @@ def coordinates_at(self, distances, crs) -> np.ndarray: """ dist = np.atleast_1d(np.asarray(distances, dtype=float)) out = np.full((len(dist), len(crs.coordinate_labels)), np.nan) - if not self.geometry: + placing = [x for x in self.geometry if axis_columns(x, crs)] + if not placing: return out - masks = interval_masks(dist, [x.interval for x in self.geometry]) - for segment, mask in zip(self.geometry, masks, strict=True): + # Only the segments which place the fiber decide the position track's + # coverage. A depth-only segment starting where one of them ends is + # not a rival for that distance, and letting it claim the run end + # would take the last surveyed point off the position it belongs to. + masks = interval_masks(dist, [x.interval for x in placing]) + for segment, mask in zip(placing, masks, strict=True): axes = axis_columns(segment, crs) - if axes and len(set(axes.values())) != len(crs.coordinate_labels): - # Half a position is not one. A checked inventory cannot get - # here; an unchecked one says so rather than handing back a - # row whose missing axis reads as uncovered distance. - msg = ( - f"Geometry {segment.name!r} states the axes {sorted(axes)} " - f"but the CRS declares {list(crs.coordinate_labels)}; a " - "segment states every axis or none of them." - ) - raise InvalidInventoryError(msg) - if not axes or not np.any(mask): + _check_axis_set(segment, axes, crs) + if not np.any(mask): continue values = segment.interpolate(dist[mask]) rows = np.flatnonzero(mask) @@ -1329,6 +1358,12 @@ def _check_geometry_columns(self) -> list[str]: "becomes a coordinate and cannot shadow a structural " "coordinate or a typed track." ) + for name in sorted(x for x in spans if "." in x): + errors.append( + f"Geometry column {name!r} states a dotted name, which is " + "how a field of a typed track is asked for; a column is " + "asked for by a name of its own." + ) groups = {x.group for x in self.annotations if x.group} for name in sorted(set(spans) & groups): errors.append( @@ -1417,14 +1452,10 @@ def select(self, *, distance: tuple[float | None, float | None]) -> Self: name: tuple(np.interp(new_dist, dist, np.asarray(values, dtype=float))) for name, values in seg.coordinates.items() } - geometry.append( - seg.model_copy( - update={ - "distance": tuple(new_dist), - "coordinates": new_coords, - } - ) - ) + # new(), not model_copy(): a copy skips the validators, and + # would leave `coordinates` a plain mutable dict whose columns + # nothing had checked against the new distance array. + geometry.append(seg.new(distance=tuple(new_dist), coordinates=new_coords)) outer = self.end_distance coupling = _clip_intervals(self.coupling, lo, hi, outer) annotations = _clip_intervals(self.annotations, lo, hi, outer) @@ -1464,14 +1495,12 @@ def flip(d): for seg in self.geometry: dist = np.asarray(seg.distance, dtype=float) geometry.append( - seg.model_copy( - update={ - "distance": tuple(flip(dist)[::-1]), - "coordinates": { - name: tuple(values[::-1]) - for name, values in seg.coordinates.items() - }, - } + seg.new( + distance=tuple(flip(dist)[::-1]), + coordinates={ + name: tuple(values[::-1]) + for name, values in seg.coordinates.items() + }, ) ) geometry.sort(key=lambda s: s.distance[0]) @@ -2225,6 +2254,7 @@ def check_width(width, what): for path in array.optical_paths: for segment in path.geometry: errors.extend(self._check_segment_axes(segment, crs)) + errors.extend(self._check_axis_overlap(path, crs)) for station in net.stations: if station.coordinates is not None: check_width(len(station.coordinates), f"Station {station.code!r}") @@ -2247,22 +2277,7 @@ def _check_segment_axes(segment, crs) -> list[str]: """ what = f"Geometry {segment.name!r}" if segment.name else "A geometry" axes = axis_columns(segment, crs) - errors = [] - spellings: dict[int, list[str]] = {} - for name, index in axes.items(): - spellings.setdefault(index, []).append(name) - for index, names in sorted(spellings.items()): - if len(names) > 1: - errors.append( - f"{what} states axis {crs.coordinate_labels[index]!r} " - f"twice, as {sorted(names)}." - ) - if axes and len(set(axes.values())) != len(crs.coordinate_labels): - errors.append( - f"{what} states the axes {sorted(axes)} but the inventory " - f"CRS declares {list(crs.coordinate_labels)}; a segment " - "states every axis or none of them." - ) + errors = _axis_set_errors(segment, axes, crs) if on_axes := sorted(set(segment.units) & set(axes)): errors.append( f"{what} states units for the axis column(s) {on_axes}; the " @@ -2270,6 +2285,31 @@ def _check_segment_axes(segment, crs) -> list[str]: ) return errors + @staticmethod + def _check_axis_overlap(path, crs) -> list[str]: + """ + Check that two segments do not place the same axis twice. + + The path checks its columns by name, which is all it can do; two + spellings of one axis are two names there and one axis here, so the + overlap has to be looked for again against what the CRS says. + """ + spans: dict[int, list[tuple[float, float]]] = {} + for segment in path.geometry: + for index in set(axis_columns(segment, crs).values()): + spans.setdefault(index, []).append(segment.interval) + errors = [] + for index in sorted(spans): + overlap = _intervals_overlap(spans[index]) + if overlap is not None: + errors.append( + f"Overlapping geometry intervals {overlap[0]} and " + f"{overlap[1]} for axis " + f"{crs.coordinate_labels[index]!r}; an axis is a " + "function track." + ) + return errors + def resolve(self, acquisition_key: str, time=None) -> ResolvedContext: """ Resolve an acquisition_key (and time) to its inventory context. diff --git a/dascore/core/inventory_loader.py b/dascore/core/inventory_loader.py index 7c8a7a1bc..970155190 100644 --- a/dascore/core/inventory_loader.py +++ b/dascore/core/inventory_loader.py @@ -992,6 +992,15 @@ def _geometry_columns(frame: pd.DataFrame, crs, path: Path): not a position, and guessing the missing axis is not a reader's job. """ labels = tuple(crs.coordinate_labels) + + def is_axis(name: str) -> bool: + """Whether the CRS reads this header as one of its own axes.""" + try: + crs.axis_index(name) + except InvalidInventoryError: + return False + return True + renamed, units = {}, {} for header in frame.columns: if header in {"segment", "distance"}: @@ -1002,25 +1011,26 @@ def _geometry_columns(frame: pd.DataFrame, crs, path: Path): renamed[header] = name if not unit: continue - if name in labels or name in {"x", "y", "z"}: + if is_axis(name): msg = ( f"{_quote(path)} states units for {name!r}, which is a " "position axis; the CRS states the units of its own axes." ) raise InvalidInventoryError(msg) units[name] = unit - stated = set(renamed.values()) - if len(stated) != len(renamed): - repeated = sorted( - {x for x in renamed.values() if list(renamed.values()).count(x) > 1} - ) + # Counted against the structural columns as well: `distance (m)` renames + # to a column the table already has, and two of them would reach pandas + # rather than this message. + written = [*renamed.values(), "segment", "distance"] + if repeated := sorted({x for x in written if written.count(x) > 1}): msg = ( f"{_quote(path)} names the column(s) {repeated} more than once; " "one column states one thing." ) raise InvalidInventoryError(msg) - axes = stated & set(labels) - if axes and axes != set(labels): + stated = set(renamed.values()) + axes = {x for x in stated if is_axis(x)} + if axes and len(axes) != len(labels): msg = ( f"{_quote(path)} states the axis column(s) {sorted(axes)}, but " f"its frame declares {list(labels)}; a segment states every axis " diff --git a/dascore/proc/inventory.py b/dascore/proc/inventory.py index a6d38c0b6..3a3e2ab05 100644 --- a/dascore/proc/inventory.py +++ b/dascore/proc/inventory.py @@ -333,10 +333,10 @@ def _get_blanket_coord_names(inventory, path) -> list[str]: """ crs = inventory.coordinate_reference_system labels = crs.coordinate_labels - out = ["x", "y", "z"][: len(labels)] if path.geometry else [] - # The axes are copied under their canonical names, so a column which is - # one of them is already covered; the rest come under their own. + # The axes are copied under their canonical names, and only where some + # segment actually places the fiber; the rest come under their own. axes = {x for segment in path.geometry for x in axis_columns(segment, crs)} + out = ["x", "y", "z"][: len(labels)] if axes else [] out += [x for x in path.geometry_columns() if x not in axes] seen = dict.fromkeys(x.group for x in path.annotations) return out + [x for x in seen if x] diff --git a/tests/test_core/test_inventory.py b/tests/test_core/test_inventory.py index ae75e9e78..23173d4f9 100644 --- a/tests/test_core/test_inventory.py +++ b/tests/test_core/test_inventory.py @@ -344,6 +344,108 @@ def test_a_column_no_segment_states(self): assert path.column_at("azimuth", [5.0]) is None +class TestGeometryColumnReviewFindings: + """What a review of the column form found, kept as regressions.""" + + @staticmethod + def _path(*geometry): + return inv.OpticalPath( + optical_components=(inv.FiberSegment(optical_length=100.0),), + geometry=geometry, + ) + + def _inventory(self, *geometry, crs=None): + array = inv.FiberArray(code="L001", optical_paths=(self._path(*geometry),)) + return inv.Inventory( + networks=(inv.Network(code="XX", fiber_arrays=(array,)),), + **({"coordinate_reference_system": crs} if crs else {}), + ) + + def test_two_spellings_of_one_axis_overlap(self): + """Checking columns by name alone would let these two through.""" + canonical = inv.Geometry( + name="a", + distance=(0.0, 10.0), + coordinates={"x": (0.0, 10.0), "y": (0.0, 10.0), "z": (0.0, 10.0)}, + ) + labelled = inv.Geometry( + name="b", + distance=(5.0, 15.0), + coordinates={ + "longitude": (100.0, 110.0), + "latitude": (100.0, 110.0), + "elevation": (100.0, 110.0), + }, + ) + with pytest.raises(InvalidInventoryError, match="for axis"): + self._inventory(canonical, labelled).check() + + def test_an_axis_stated_twice_is_refused_when_placing(self): + """Otherwise whichever spelling came last would win, silently.""" + doubled = inv.Geometry( + name="doubled", + distance=(0.0, 10.0), + coordinates={ + "x": (0.0, 10.0), + "longitude": (100.0, 110.0), + "y": (0.0, 1.0), + "z": (0.0, 1.0), + }, + ) + crs = inv.CoordinateReferenceSystem() + with pytest.raises(InvalidInventoryError, match="twice"): + self._path(doubled).coordinates_at([5.0], crs) + + def test_a_column_segment_does_not_claim_a_position_run_end(self): + """The position track's coverage is the segments which place it.""" + placed = inv.Geometry( + name="placed", + distance=(0.0, 10.0), + coordinates={"x": (0.0, 1.0), "y": (0.0, 1.0), "z": (0.0, 1.0)}, + ) + measured = inv.Geometry( + name="measured", + distance=(10.0, 20.0), + coordinates={"borehole_depth": (0.0, 10.0)}, + ) + crs = inv.CoordinateReferenceSystem() + out = self._path(placed, measured).coordinates_at([10.0], crs) + assert np.allclose(out, [[1.0, 1.0, 1.0]]) + + def test_a_dotted_column_name(self): + """A dotted name is how a field of a typed track is asked for.""" + dotted = inv.Geometry( + distance=(0.0, 10.0), coordinates={"coupling.depth": (0.0, 1.0)} + ) + with pytest.raises(InvalidInventoryError, match="dotted name"): + self._inventory(dotted).check() + + def test_select_revalidates_the_columns_it_writes(self): + """model_copy would skip the validators and leave a mutable dict.""" + segment = inv.Geometry( + name="run", + distance=(0.0, 10.0), + coordinates={"x": (0.0, 1.0), "y": (0.0, 1.0), "z": (0.0, 1.0)}, + ) + out = self._path(segment).select(distance=(2.0, 8.0)) + coordinates = out.geometry[0].coordinates + assert isinstance(coordinates, Mapping) + with pytest.raises(TypeError): + coordinates["x"] = (0.0,) + + def test_a_canonical_name_the_crs_has_no_axis_for(self): + """`z` is a column of its own where the CRS declares two axes.""" + crs = inv.CoordinateReferenceSystem( + coordinate_labels=("easting", "northing"), units=("meter", "meter") + ) + segment = inv.Geometry( + distance=(0.0, 10.0), coordinates={"z": (0.0, 1.0)}, units={"z": "m"} + ) + inventory = self._inventory(segment, crs=crs) + assert inventory.check() is inventory + assert "z" in inventory.get_names().coords + + class TestGeometry: """Geometry segment rules.""" diff --git a/tests/test_core/test_inventory_loader.py b/tests/test_core/test_inventory_loader.py index 030d6d58e..9a6c575e2 100644 --- a/tests/test_core/test_inventory_loader.py +++ b/tests/test_core/test_inventory_loader.py @@ -1220,6 +1220,37 @@ def test_a_column_of_text_is_refused(self, make_inventory): with pytest.raises(InvalidInventoryError, match=r"annotations\.csv"): make_inventory(files) + def test_a_structural_column_restated_with_units(self, make_inventory): + """`distance (m)` renames onto a column the table already has.""" + files = { + **MINIMAL, + **TRACKS, + "fiber_arrays/DAS.L001/path/geometry.csv": ( + "segment,distance,distance (m)\nS100,100.0,100.0\nS100,102.0,102.0\n" + ), + } + with pytest.raises(InvalidInventoryError, match="more than once"): + make_inventory(files) + + def test_a_canonical_name_the_crs_has_no_axis_for(self, make_inventory): + """`z` is a column of its own where the CRS declares two axes.""" + files = { + **MINIMAL, + **TRACKS, + "inventory.yaml": ( + "object_type: Inventory\n" + "coordinate_reference_system:\n" + " coordinate_labels: [easting, northing]\n" + " units: [meter, meter]\n" + ), + "fiber_arrays/DAS.L001/path/geometry.csv": ( + "segment,distance,z (m)\nS100,100.0,0.0\nS100,102.0,2.0\n" + ), + } + geometry = one_path(make_inventory(files)).geometry[0] + assert geometry.coordinates["z"] == (0.0, 2.0) + assert geometry.units["z"] == "m" + def test_one_column_stated_twice(self, make_inventory): """A unit suffix is not a second column, however it is spelled.""" files = { diff --git a/tests/test_proc/test_proc_inventory.py b/tests/test_proc/test_proc_inventory.py index fdfbfe44c..de5349f6f 100644 --- a/tests/test_proc/test_proc_inventory.py +++ b/tests/test_proc/test_proc_inventory.py @@ -376,6 +376,21 @@ def test_selection_trims_channels(self, patch, inventory): # acquisition maps path distance 100 onto channel 0. assert depth.max() < patch.get_coord("distance").max() + def test_an_axis_is_missing_where_nothing_places_the_fiber(self, patch, inventory): + """Geometry which measures but does not place defines no axis.""" + segment = Geometry( + name="hole", + distance=(100.0, 200.0), + coordinates={"borehole_depth": (0.0, 100.0)}, + ) + inv = _replace_path(inventory, geometry=(segment,)) + with pytest.raises(PatchError, match="defines no 'x'"): + patch.enrich(inv, attrs=False, coords=("x",)) + # And a blanket request does not offer one either. + out = patch.enrich(inv, attrs=False) + assert "x" not in set(out.coords.coord_map) + assert "borehole_depth" in set(out.coords.coord_map) + def test_a_column_does_not_bridge_segments(self, patch, inventory): """Two holes are two holes, and the fiber between them is neither.""" first = Geometry( From 674d94dc7da26e3c2a19134d0c05c7306ec8797b Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Sat, 15 Aug 2026 19:28:34 +0200 Subject: [PATCH 3/5] Make one stretch of fiber have one geometry name Letting two segments overlap where they state different columns left the bare `geometry` coordinate ambiguous: it is each segment's name over its interval, so a channel covered by a depth survey and an azimuth survey took whichever name the tuple held last, and reversing the tuple changed what `select(geometry=...)` matched. Segments which overlap are two measurements of one stretch of fiber, so they state its name -- both "hole 1" rather than "depth survey" and "azimuth survey". Overlapping segments whose names differ are refused, which leaves the identity single-valued without taking the overlap away. --- dascore/core/inventory.py | 32 ++++++++++++++++++++++++++++++- tests/test_core/test_inventory.py | 21 ++++++++++++++++++-- 2 files changed, 50 insertions(+), 3 deletions(-) diff --git a/dascore/core/inventory.py b/dascore/core/inventory.py index b82a6b873..41412bef6 100644 --- a/dascore/core/inventory.py +++ b/dascore/core/inventory.py @@ -526,7 +526,10 @@ class Geometry(InventoryModel): azimuth. Those carry their own entry in ``units``. Interpolation between points is piecewise linear and never crosses - segments; uncovered distance has undefined values. + segments; uncovered distance has undefined values. Two segments may + cover the same distance as long as they state different columns -- a + depth survey and an azimuth survey of one borehole -- and then they + share a name, being two measurements of one stretch of fiber. Examples -------- @@ -1370,6 +1373,7 @@ def _check_geometry_columns(self) -> list[str]: f"{name!r} is both a geometry column and an annotation " "group; one name is one coordinate." ) + errors.extend(self._check_overlapping_names()) for name in sorted(spans): overlap = _intervals_overlap(spans[name]) if overlap is not None: @@ -1385,6 +1389,32 @@ def _check_geometry_columns(self) -> list[str]: ) return errors + def _check_overlapping_names(self) -> list[str]: + """ + Check that segments covering one stretch of fiber share a name. + + Two segments may cover the same distance now, as long as they state + different columns -- a depth survey and an azimuth survey of one + borehole. They are two measurements of one stretch of fiber, so they + are one segment by name: the bare ``geometry`` coordinate is that + name, and a channel with two of them would take whichever the tuple + happened to hold last. + """ + errors = [] + for first, second in itertools.combinations(self.geometry, 2): + if first.name == second.name: + continue + lo = max(first.interval[0], second.interval[0]) + hi = min(first.interval[1], second.interval[1]) + if lo < hi: + errors.append( + f"Geometry segments {first.name!r} and {second.name!r} " + f"both cover ({lo}, {hi}); segments which overlap state " + "different columns of one stretch of fiber, so they " + "share its name." + ) + return errors + def _check_annotation_groups(self) -> list[str]: """Check that each annotation group holds one kind of value.""" groups: dict[str, list] = {} diff --git a/tests/test_core/test_inventory.py b/tests/test_core/test_inventory.py index 23173d4f9..2f994cf56 100644 --- a/tests/test_core/test_inventory.py +++ b/tests/test_core/test_inventory.py @@ -256,13 +256,30 @@ def test_overlap_is_refused_per_column(self): def test_different_columns_may_overlap(self): """Each column is its own function track, so they are independent.""" - depth = inv.Geometry(distance=(0.0, 60.0), coordinates={"depth": (0.0, 6.0)}) + depth = inv.Geometry( + name="hole 1", distance=(0.0, 60.0), coordinates={"depth": (0.0, 6.0)} + ) azimuth = inv.Geometry( - distance=(50.0, 80.0), coordinates={"azimuth": (5.0, 8.0)} + name="hole 1", distance=(50.0, 80.0), coordinates={"azimuth": (5.0, 8.0)} ) inventory = self._inventory(depth, azimuth) assert inventory.check() is inventory + def test_overlapping_segments_share_a_name(self): + """The bare `geometry` coordinate is that name, so there is one.""" + depth = inv.Geometry( + name="depth survey", + distance=(0.0, 60.0), + coordinates={"depth": (0.0, 6.0)}, + ) + azimuth = inv.Geometry( + name="azimuth survey", + distance=(50.0, 80.0), + coordinates={"azimuth": (5.0, 8.0)}, + ) + with pytest.raises(InvalidInventoryError, match="share its name"): + self._inventory(depth, azimuth).check() + def test_a_reserved_column_name(self): """A column becomes a coordinate, so it cannot shadow one.""" clash = inv.Geometry(distance=(0.0, 10.0), coordinates={"time": (0.0, 1.0)}) From 1c1f22602d70927f02715f7efdb53bd8eac73703 Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Sat, 15 Aug 2026 19:45:17 +0200 Subject: [PATCH 4/5] Say the geometry rules once each The validation had grown a helper per finding rather than a function per idea. The path's rules are one function now and the inventory's are another, split where they have to be: the path checks its columns by name, and the CRS-dependent ones can only be checked where the CRS is. Two other duplications went with them. `coordinates_at` and `column_at` each carried their own copy of "interpolate this column, then fill the run end the mask includes from the last control point", which is now one helper. The loader read its headers and then coerced them in two passes over the same columns, which is now one. `column_units` and the separate raise-wrapper around the axis-set rule each had a single caller and are inlined. The prose came down too, mostly by not saying in a field description what the class docstring above it already says. --- dascore/core/_spool_inventory.py | 5 +- dascore/core/inventory.py | 239 ++++++++++++------------------ dascore/core/inventory_loader.py | 47 +++--- tests/test_core/test_inventory.py | 2 +- 4 files changed, 115 insertions(+), 178 deletions(-) diff --git a/dascore/core/_spool_inventory.py b/dascore/core/_spool_inventory.py index 55cbc2e16..7756bbb44 100644 --- a/dascore/core/_spool_inventory.py +++ b/dascore/core/_spool_inventory.py @@ -682,7 +682,10 @@ def _get_geometry_column_coord(path, name, distances): values = path.column_at(name, distances) if values is None: return None - return get_coord(data=values, units=path.column_units(name) or None) + # The path has already refused two segments measuring one column in two + # units, so whichever states it states the one they agree on. + units = next((x.units[name] for x in path.geometry if name in x.units), None) + return get_coord(data=values, units=units) def get_coord_values(inventory, path, name, distances): diff --git a/dascore/core/inventory.py b/dascore/core/inventory.py index 7874af967..bda48adff 100644 --- a/dascore/core/inventory.py +++ b/dascore/core/inventory.py @@ -519,11 +519,11 @@ class Geometry(InventoryModel): "clump", is a segment whose columns repeat while distance advances. A column whose name the inventory CRS declares -- or the canonical - ``x``, ``y``, ``z`` alias of one -- is that position axis, and its units - are the CRS's. Every other column is a numeric quantity along the fiber - in its own right: borehole depth where the CRS is spent on - easting/northing/elevation, pipeline chainage, burial depth, fiber - azimuth. Those carry their own entry in ``units``. + ``x``, ``y``, ``z`` alias of one -- is that position axis, and takes the + CRS's units. Every other column is a numeric quantity along the fiber in + its own right, carrying its own entry in ``units``: borehole depth where + the CRS is spent on easting/northing/elevation, pipeline chainage, fiber + azimuth. Interpolation between points is piecewise linear and never crosses segments; uncovered distance has undefined values. Two segments may @@ -562,16 +562,12 @@ class Geometry(InventoryModel): coordinates: FrozenDictType[str, tuple[float, ...]] = Field( description=( "Columns measured along this segment, keyed by name; each holds " - "one value per distance. A name the CRS declares is that " - "position axis, any other is a numeric column of its own." + "one value per distance." ), ) units: FrozenDictType[str, str] = Field( default_factory=dict, - description=( - "Units of the columns which are not position axes; the CRS " - "states the units of the axes." - ), + description="Units of the columns which are not position axes.", ) @model_validator(mode="after") @@ -1087,6 +1083,19 @@ def axis_columns(segment, crs) -> dict[str, int]: return out +def _placed(segment, name: str, distances) -> np.ndarray: + """ + One column of a segment, over distances already claimed for it. + + ``interpolate`` stops at its own half-open coverage, which excludes the + run end a mask includes, so that point is filled from the last control + point rather than left NaN. + """ + column = segment.interpolate(distances)[name] + column[np.isnan(column)] = segment.coordinates[name][-1] + return column + + def _axis_set_errors(segment, axes: Mapping[str, int], crs) -> list[str]: """Return what is wrong with the set of axes one segment states.""" what = f"Geometry {segment.name!r}" if segment.name else "A geometry" @@ -1108,18 +1117,6 @@ def _axis_set_errors(segment, axes: Mapping[str, int], crs) -> list[str]: return errors -def _check_axis_set(segment, axes: Mapping[str, int], crs) -> None: - """ - Refuse a segment whose axes are partial or doubly spelled. - - A checked inventory cannot reach this; an unchecked one says so rather - than filling a position from whichever spelling the mapping happened to - hold last, or leaving an axis reading as uncovered distance. - """ - if errors := _axis_set_errors(segment, axes, crs): - raise InvalidInventoryError(" ".join(errors)) - - def _track_identity_fields() -> Mapping[str, str]: """ Map each typed track of an optical path to the field its name means. @@ -1258,44 +1255,30 @@ def coordinates_at(self, distances, crs) -> np.ndarray: """ Return CRS coordinates at the requested optical distances. - Only the columns which name a position axis are assembled; a segment - stating none of them contributes no position, however much else it - measures. Uncovered distance returns NaN rows. Segment coverage is - half-open, with the end of each coverage run included: a distance on - a segment's last control point belongs to that segment unless - another segment claims it. - - Parameters - ---------- - distances - The optical distances to place. - crs - The inventory's coordinate reference system, which is what - decides that a column is an axis and which axis it is. + The CRS is what decides a column is an axis, and only those are + assembled: a segment naming none of them contributes no position, + however much else it measures, and does not decide the position + track's coverage either. Uncovered distance returns NaN rows. + Coverage is half-open, with the end of each run included: a distance + on a segment's last control point belongs to it unless another + segment claims it. """ dist = np.atleast_1d(np.asarray(distances, dtype=float)) out = np.full((len(dist), len(crs.coordinate_labels)), np.nan) placing = [x for x in self.geometry if axis_columns(x, crs)] - if not placing: - return out - # Only the segments which place the fiber decide the position track's - # coverage. A depth-only segment starting where one of them ends is - # not a rival for that distance, and letting it claim the run end - # would take the last surveyed point off the position it belongs to. masks = interval_masks(dist, [x.interval for x in placing]) for segment, mask in zip(placing, masks, strict=True): axes = axis_columns(segment, crs) - _check_axis_set(segment, axes, crs) + if errors := _axis_set_errors(segment, axes, crs): + # A checked inventory cannot get here; an unchecked one says + # so rather than filling a position from whichever spelling + # the mapping held last. + raise InvalidInventoryError(" ".join(errors)) if not np.any(mask): continue - values = segment.interpolate(dist[mask]) rows = np.flatnonzero(mask) for name, index in axes.items(): - column = values[name] - # interpolate() reports its own coverage, which excludes the - # run end the mask includes, so fill that from the last point. - column[np.isnan(column)] = segment.coordinates[name][-1] - out[rows, index] = column + out[rows, index] = _placed(segment, name, dist[mask]) return out def column_at(self, name: str, distances) -> np.ndarray | None: @@ -1313,18 +1296,10 @@ def column_at(self, name: str, distances) -> np.ndarray | None: out = np.full(len(dist), np.nan) masks = interval_masks(dist, [x.interval for x in stating]) for segment, mask in zip(stating, masks, strict=True): - if not np.any(mask): - continue - column = segment.interpolate(dist[mask])[name] - column[np.isnan(column)] = segment.coordinates[name][-1] - out[np.flatnonzero(mask)] = column + if np.any(mask): + out[np.flatnonzero(mask)] = _placed(segment, name, dist[mask]) return out - def column_units(self, name: str) -> str: - """Return the units the segments stating a column agree on.""" - stated = {x.units[name] for x in self.geometry if name in x.units} - return stated.pop() if len(stated) == 1 else "" - def geometry_columns(self) -> tuple[str, ...]: """Return every column name this path's geometry segments state.""" seen: dict[str, None] = {} @@ -1375,10 +1350,11 @@ def _check_geometry_columns(self) -> list[str]: Each column is its own function track, so two segments may overlap as long as they do not state the same column over the same distance. - The rules which need the CRS -- which columns are axes -- are the - inventory's, since only it knows what the axes are. + Overlapping segments do share a name, being two measurements of one + stretch of fiber, or the bare ``geometry`` coordinate would have two + of them for a channel. The rules needing the CRS -- which columns + are axes -- are the inventory's, since only it knows the axes. """ - errors = [] spans: dict[str, list[tuple[float, float]]] = {} units: dict[str, set[str]] = {} for segment in self.geometry: @@ -1386,64 +1362,46 @@ def _check_geometry_columns(self) -> list[str]: spans.setdefault(name, []).append(segment.interval) if name in segment.units: units.setdefault(name, set()).add(segment.units[name]) - for name in sorted(set(spans) & _RESERVED_COLUMN_NAMES): - errors.append( - f"Geometry column {name!r} is a reserved name; a column " - "becomes a coordinate and cannot shadow a structural " - "coordinate or a typed track." - ) - for name in sorted(x for x in spans if "." in x): - errors.append( - f"Geometry column {name!r} states a dotted name, which is " - "how a field of a typed track is asked for; a column is " - "asked for by a name of its own." - ) groups = {x.group for x in self.annotations if x.group} - for name in sorted(set(spans) & groups): - errors.append( - f"{name!r} is both a geometry column and an annotation " - "group; one name is one coordinate." - ) - errors.extend(self._check_overlapping_names()) - for name in sorted(spans): - overlap = _intervals_overlap(spans[name]) - if overlap is not None: - errors.append( - f"Overlapping geometry intervals {overlap[0]} and " - f"{overlap[1]} for column {name!r}; a column is a " - "function track." - ) - if len(stated := units.get(name, set())) > 1: - errors.append( - f"Geometry column {name!r} is stated in " - f"{sorted(stated)}; a column has one unit." - ) - return errors - - def _check_overlapping_names(self) -> list[str]: - """ - Check that segments covering one stretch of fiber share a name. - - Two segments may cover the same distance now, as long as they state - different columns -- a depth survey and an azimuth survey of one - borehole. They are two measurements of one stretch of fiber, so they - are one segment by name: the bare ``geometry`` coordinate is that - name, and a channel with two of them would take whichever the tuple - happened to hold last. - """ - errors = [] + errors = [ + f"Geometry column {name!r} is a reserved name; a column becomes " + "a coordinate and cannot shadow a structural coordinate or a " + "typed track." + for name in sorted(set(spans) & _RESERVED_COLUMN_NAMES) + ] + errors += [ + f"Geometry column {name!r} states a dotted name, which is how a " + "field of a typed track is asked for; a column is asked for by a " + "name of its own." + for name in sorted(x for x in spans if "." in x) + ] + errors += [ + f"{name!r} is both a geometry column and an annotation group; " + "one name is one coordinate." + for name in sorted(set(spans) & groups) + ] for first, second in itertools.combinations(self.geometry, 2): - if first.name == second.name: - continue lo = max(first.interval[0], second.interval[0]) hi = min(first.interval[1], second.interval[1]) - if lo < hi: + if first.name != second.name and lo < hi: errors.append( f"Geometry segments {first.name!r} and {second.name!r} " f"both cover ({lo}, {hi}); segments which overlap state " "different columns of one stretch of fiber, so they " "share its name." ) + for name in sorted(spans): + if (overlap := _intervals_overlap(spans[name])) is not None: + errors.append( + f"Overlapping geometry intervals {overlap[0]} and " + f"{overlap[1]} for column {name!r}; a column is a " + "function track." + ) + if len(stated := units.get(name, set())) > 1: + errors.append( + f"Geometry column {name!r} is stated in {sorted(stated)}; " + "a column has one unit." + ) return errors def _check_annotation_groups(self) -> list[str]: @@ -2313,9 +2271,7 @@ def check_width(width, what): for net in self.networks: for array in net.fiber_arrays: for path in array.optical_paths: - for segment in path.geometry: - errors.extend(self._check_segment_axes(segment, crs)) - errors.extend(self._check_axis_overlap(path, crs)) + errors.extend(self._check_geometry_axes(path, crs)) for station in net.stations: if station.coordinates is not None: check_width(len(station.coordinates), f"Station {station.code!r}") @@ -2327,47 +2283,36 @@ def check_width(width, what): return errors @staticmethod - def _check_segment_axes(segment, crs) -> list[str]: - """ - Check one geometry segment's columns against the CRS. - - A segment states every position axis or none of them: a partial - position is not one, and deciding what the missing axis meant is not - the reader's job. Which columns those are is the CRS's to say, which - is why this lives here rather than on the path. - """ - what = f"Geometry {segment.name!r}" if segment.name else "A geometry" - axes = axis_columns(segment, crs) - errors = _axis_set_errors(segment, axes, crs) - if on_axes := sorted(set(segment.units) & set(axes)): - errors.append( - f"{what} states units for the axis column(s) {on_axes}; the " - "CRS states the units of its own axes." - ) - return errors - - @staticmethod - def _check_axis_overlap(path, crs) -> list[str]: + def _check_geometry_axes(path, crs) -> list[str]: """ - Check that two segments do not place the same axis twice. - - The path checks its columns by name, which is all it can do; two - spellings of one axis are two names there and one axis here, so the - overlap has to be looked for again against what the CRS says. + Check one path's geometry against the CRS. + + Which columns are axes is the CRS's to say, so the rules needing it + live here rather than on the path: that a segment states every axis + or none, that it does not spell one twice, that it leaves the axes' + units to the CRS, and that two segments do not place the same axis + over one distance -- which the path cannot see, two spellings of an + axis being two names to it. """ + errors = [] spans: dict[int, list[tuple[float, float]]] = {} for segment in path.geometry: - for index in set(axis_columns(segment, crs).values()): + axes = axis_columns(segment, crs) + errors += _axis_set_errors(segment, axes, crs) + if on_axes := sorted(set(segment.units) & set(axes)): + what = f"Geometry {segment.name!r}" if segment.name else "A geometry" + errors.append( + f"{what} states units for the axis column(s) {on_axes}; " + "the CRS states the units of its own axes." + ) + for index in set(axes.values()): spans.setdefault(index, []).append(segment.interval) - errors = [] for index in sorted(spans): - overlap = _intervals_overlap(spans[index]) - if overlap is not None: + if (overlap := _intervals_overlap(spans[index])) is not None: errors.append( f"Overlapping geometry intervals {overlap[0]} and " - f"{overlap[1]} for axis " - f"{crs.coordinate_labels[index]!r}; an axis is a " - "function track." + f"{overlap[1]} for axis {crs.coordinate_labels[index]!r}; " + "an axis is a function track." ) return errors diff --git a/dascore/core/inventory_loader.py b/dascore/core/inventory_loader.py index 970155190..7b4bd99ec 100644 --- a/dascore/core/inventory_loader.py +++ b/dascore/core/inventory_loader.py @@ -986,12 +986,13 @@ def _geometry_columns(frame: pd.DataFrame, crs, path: Path): """ Read a geometry table's headers, and refuse what cannot be a column. - A header naming an axis the CRS declares is that axis; every other one - is a numeric column in its own right, which may carry its units in - parentheses. The axes are all stated or none are: a partial position is - not a position, and guessing the missing axis is not a reader's job. + A header naming an axis the CRS declares is that axis and takes the + CRS's units; every other one is a numeric column in its own right, and + may carry its units in parentheses. The axes are all stated or none + are, since guessing the missing one is not a reader's job. Text is + refused: a value which varies along the fiber without being a number + is what annotations are for. """ - labels = tuple(crs.coordinate_labels) def is_axis(name: str) -> bool: """Whether the CRS reads this header as one of its own axes.""" @@ -1001,6 +1002,7 @@ def is_axis(name: str) -> bool: return False return True + labels = tuple(crs.coordinate_labels) renamed, units = {}, {} for header in frame.columns: if header in {"segment", "distance"}: @@ -1009,15 +1011,14 @@ def is_axis(name: str) -> bool: if (match := _UNIT_SUFFIX.match(header)) is not None: name, unit = match.group("name"), match.group("units").strip() renamed[header] = name - if not unit: - continue - if is_axis(name): + if unit and is_axis(name): msg = ( f"{_quote(path)} states units for {name!r}, which is a " "position axis; the CRS states the units of its own axes." ) raise InvalidInventoryError(msg) - units[name] = unit + if unit: + units[name] = unit # Counted against the structural columns as well: `distance (m)` renames # to a column the table already has, and two of them would reach pandas # rather than this message. @@ -1028,8 +1029,7 @@ def is_axis(name: str) -> bool: "one column states one thing." ) raise InvalidInventoryError(msg) - stated = set(renamed.values()) - axes = {x for x in stated if is_axis(x)} + axes = {x for x in renamed.values() if is_axis(x)} if axes and len(axes) != len(labels): msg = ( f"{_quote(path)} states the axis column(s) {sorted(axes)}, but " @@ -1037,30 +1037,19 @@ def is_axis(name: str) -> bool: "or none of them." ) raise InvalidInventoryError(msg) - frame = frame.rename(columns=renamed) - return _numeric_columns(frame, sorted(stated), path), units - - -def _numeric_columns(frame: pd.DataFrame, columns, path: Path) -> pd.DataFrame: - """ - Read a geometry table's columns as numbers, refusing text. - - Text along distance is what annotations are for, and a column of it - here would otherwise reach the model as a string it cannot place. - """ - frame = frame.copy() - for column in columns: + frame = frame.rename(columns=renamed).copy() + for column in sorted(set(renamed.values())): values = pd.to_numeric(frame[column], errors="coerce") if (bad := frame[column].notna() & values.isna()).any(): - first = frame.loc[bad, column].iloc[0] msg = ( - f"{_quote(path)} states {first!r} in column {column!r}, " - "which is not a number. A geometry column is numeric; text " - "which varies along the fiber belongs in annotations.csv." + f"{_quote(path)} states {frame.loc[bad, column].iloc[0]!r} in " + f"column {column!r}, which is not a number. A geometry column " + "is numeric; text which varies along the fiber belongs in " + "annotations.csv." ) raise InvalidInventoryError(msg) frame[column] = values - return frame + return frame, units def _parse_annotations(rows: list[dict], path: Path) -> None: diff --git a/tests/test_core/test_inventory.py b/tests/test_core/test_inventory.py index ca85eaaf2..64cfb9ee7 100644 --- a/tests/test_core/test_inventory.py +++ b/tests/test_core/test_inventory.py @@ -308,7 +308,7 @@ def test_units_on_an_axis_are_refused(self): def test_units_for_a_column_which_is_not_there(self): """A unit sitting on nothing is a typo no reader would find.""" - with pytest.raises(ValidationError, match="no column for"): + with pytest.raises(ValidationError, match="has no column"): inv.Geometry( distance=(0.0, 10.0), coordinates={"depth": (0.0, 1.0)}, From 732b8c1998a7cbb6af832c83fe47061b32a076e4 Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Sat, 15 Aug 2026 20:32:13 +0200 Subject: [PATCH 5/5] Drop a guard which can no longer fire `coordinates_at` used to return as many columns as the geometry happened to state, so an axis the segments did not reach had to be checked for. It returns one column per axis the CRS declares now, and `axis_index` refuses a label the CRS has no axis for, so the index is always one of those columns and the guard was unreachable -- which is what took project coverage off 100%. --- dascore/core/_spool_inventory.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dascore/core/_spool_inventory.py b/dascore/core/_spool_inventory.py index 7756bbb44..2831f30b0 100644 --- a/dascore/core/_spool_inventory.py +++ b/dascore/core/_spool_inventory.py @@ -671,9 +671,9 @@ def _get_geometry_coord(inventory, path, label, distances): # on_missing policy rules rather than a column of nan. if not any(axis_columns(x, crs) for x in path.geometry): return None + # axis_index refuses a label this CRS has no axis for, and the array is + # one column per axis, so the index is always one of its columns. coords = path.coordinates_at(distances, crs) - if index >= coords.shape[1]: - return None return get_coord(data=coords[:, index], units=crs.units[index])