diff --git a/dascore/core/_spool_inventory.py b/dascore/core/_spool_inventory.py index bb8190467..242749dc1 100644 --- a/dascore/core/_spool_inventory.py +++ b/dascore/core/_spool_inventory.py @@ -33,6 +33,7 @@ VALID_COORDINATE_LABELS, Inventory, ResolvedContext, + axis_columns, ) from dascore.core.inventory_loader import BLESSED_NAME, find_inventory from dascore.exceptions import ( @@ -665,14 +666,27 @@ def _get_geometry_coord(inventory, path, label, distances): index = crs.axis_index(label) except InvalidInventoryError: return None - if not path.geometry: - return None - coords = path.coordinates_at(distances) - if index >= coords.shape[1]: + # 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 + # 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) 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 + # 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): """Return the values of one requested coordinate, or None if undefined.""" if name == "distance": @@ -687,7 +701,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 45f87269a..babce305a 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 ( @@ -506,43 +507,92 @@ 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 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 + 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 + -------- + >>> 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." + ), + ) + units: FrozenDictType[str, str] = Field( + default_factory=dict, + description="Units of the columns which are not position 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 @@ -551,23 +601,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 @@ -988,7 +1040,53 @@ 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 _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" + 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 _track_identity_fields() -> Mapping[str, str]: @@ -1026,6 +1124,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) @@ -1098,34 +1202,62 @@ 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 + 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)) - 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) - masks = interval_masks(dist, [x.interval for x in self.geometry]) - for segment, mask in zip(self.geometry, masks, strict=True): + out = np.full((len(dist), len(crs.coordinate_labels)), np.nan) + placing = [x for x in self.geometry if axis_columns(x, crs)] + 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 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 - # 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 + rows = np.flatnonzero(mask) + for name, index in axes.items(): + out[rows, index] = _placed(segment, name, dist[mask]) + 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 np.any(mask): + out[np.flatnonzero(mask)] = _placed(segment, name, dist[mask]) return out + 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. @@ -1150,22 +1282,79 @@ 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. + 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. + """ + 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]) + groups = {x.group for x in self.annotations if x.group} + 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): + lo = max(first.interval[0], second.interval[0]) + hi = min(first.interval[1], second.interval[1]) + 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]: """Check that each annotation group holds one kind of value.""" groups: dict[str, list] = {} @@ -1229,22 +1418,14 @@ 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, - ) - geometry.append( - seg.model_copy( - update={ - "distance": tuple(new_dist), - "coordinates": tuple(map(tuple, new_coords)), - } - ) - ) + new_coords = { + name: tuple(np.interp(new_dist, dist, np.asarray(values, dtype=float))) + for name, values in seg.coordinates.items() + } + # 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) @@ -1283,13 +1464,13 @@ 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])), - } + 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]) @@ -1957,8 +2138,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, {}) @@ -1975,6 +2165,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) @@ -2008,9 +2199,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: - what = f"Geometry {segment.name!r}" - check_width(len(segment.coordinates[0]), what) + 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}") @@ -2021,6 +2210,40 @@ def check_width(width, what): check_width(len(channel.coordinates), what) return errors + @staticmethod + def _check_geometry_axes(path, crs) -> list[str]: + """ + 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: + 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) + for index in sorted(spans): + if (overlap := intervals_overlap(spans[index])) is not None: + errors.append( + f"Overlapping geometry intervals {overlap[0]} and " + f"{overlap[1]} for axis {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 5377f2360..15af9fd35 100644 --- a/dascore/core/inventory_loader.py +++ b/dascore/core/inventory_loader.py @@ -477,6 +477,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 @@ -487,7 +490,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"), } @@ -535,14 +540,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) @@ -558,8 +567,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(): @@ -584,35 +594,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 @@ -835,35 +828,93 @@ def _read_track_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 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. """ + + 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 + 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 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) + 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. + written = [*renamed.values(), "segment", "distance"] + if repeated := sorted({x for x in written if written.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 = {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 " + f"its frame declares {list(labels)}; a segment states every axis " + "or none of them." + ) + raise InvalidInventoryError(msg) + 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(): + msg = ( + 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, units 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..3a3e2ab05 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 - out = ["x", "y", "z"][: len(labels)] if path.geometry else [] + crs = inventory.coordinate_reference_system + labels = crs.coordinate_labels + # 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/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 e2ed18b93..64cfb9ee7 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,324 @@ 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( + name="hole 1", distance=(0.0, 60.0), coordinates={"depth": (0.0, 6.0)} + ) + azimuth = inv.Geometry( + 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)}) + 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="has no column"): + 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 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.""" 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 +528,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 +631,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 +1048,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 +1301,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 +1374,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 +1446,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 +2077,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 +2097,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 +2105,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 +2147,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..9a6c575e2 100644 --- a/tests/test_core/test_inventory_loader.py +++ b/tests/test_core/test_inventory_loader.py @@ -1166,13 +1166,112 @@ 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_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 = { + **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 +1293,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 +1908,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..de5349f6f 100644 --- a/tests/test_proc/test_proc_inventory.py +++ b/tests/test_proc/test_proc_inventory.py @@ -322,6 +322,96 @@ 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_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( + 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 +626,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 +777,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