From f4f0544cc9aa7385ff21c13edce8abffe73b606d Mon Sep 17 00:00:00 2001 From: "William Mark Katzenmeyer, P.E., C.F.M." Date: Sun, 30 Aug 2026 23:16:04 -0400 Subject: [PATCH 1/5] Add exact 2D geometry rebuild and containment validation --- docs/api/geometry.md | 19 + docs/user-guide/geometry-operations.md | 126 +++ ras_commander/__init__.py | 6 +- ras_commander/geom/GeomBcLines.py | 37 +- ras_commander/geom/GeomMesh.py | 824 ++++++++++++++- ras_commander/geom/GeomMeshDataclasses.py | 64 +- ras_commander/geom/GeomParser.py | 13 +- ras_commander/geom/GeomReferenceFeatures.py | 191 +++- ras_commander/geom/GeomStorage.py | 229 ++++- ras_commander/geom/__init__.py | 11 +- .../gui/workflows/mesh_regeneration.py | 960 ++++++++++++++++-- ras_commander/hdf/HdfBndry.py | 36 +- ras_commander/hdf/HdfMesh.py | 42 +- ras_commander/hdf/HdfResultsPlan.py | 58 +- setup.py | 6 +- tests/test_geom_bc_lines.py | 19 +- tests/test_geom_mesh.py | 540 +++++++++- tests/test_geom_reference_features.py | 91 ++ .../test_geom_storage_2d_flow_area_writer.py | 109 +- tests/test_gui_mesh_regeneration_exact.py | 451 ++++++++ tests/test_hdf_bndry_logging.py | 42 + tests/test_hdf_results_plan_logging.py | 39 + 22 files changed, 3731 insertions(+), 182 deletions(-) create mode 100644 tests/test_gui_mesh_regeneration_exact.py diff --git a/docs/api/geometry.md b/docs/api/geometry.md index 5ded65a20..190a6ac8a 100644 --- a/docs/api/geometry.md +++ b/docs/api/geometry.md @@ -163,6 +163,7 @@ reference-line output. - `add_reference_lines(geom_file, lines, storage_area)` - Insert manually supplied reference lines into a `.g##` file +- `replace_reference_lines(geom_file, lines, storage_area, expected_existing_names=...)` - Atomically replace or remove one 2D area's complete reference-line collection while preserving other areas, with an optional ordered concurrency guard - `generate_reference_lines_from_longitudinal_line(...)` - Generate transverse reference-line dictionaries at regular station intervals along a named longitudinal line @@ -199,10 +200,17 @@ back to normal-to-line orientation unless `orientation_fallback="raise"` is set. Headless 2D mesh generation helpers and compiled geometry HDF refinement-region utilities. +### Domain and Mesh Methods + +- `audit_domain_containment(geom_number, mesh_name=..., cell_size=..., ras_object=...)` - Fail closed unless every breakline, refinement region, and structure associated with the selected 2D area is wholly covered by the exact compiled perimeter buffered **inward** by one base mesh-cell spacing. BC lines are intentionally excluded because they are authored on the perimeter and require a separate association/overlap audit. +- `generate(geom_number, mesh_name=..., ras_object=...)` - Regenerate the mesh and automatically run the same inward one-cell containment gate before loading native RAS Mapper dependencies. +- `compute_property_tables(geom_number, mesh_name=..., ras_object=...)` - Compute face profiles, Manning's n assignments, face hydraulic tables, and cell properties against the restored geometry associations. + ### Refinement Region Methods - `add_refinement_region(geom_number, polygon, spacing_dx, ...)` - Add one refinement polygon to an existing compiled geometry HDF. - `add_flowline_refinement_regions(geom_number, flowlines, buffer_width, ...)` - Buffer GeoDataFrame or LineString channel flowlines into refinement-region polygons, optionally simplify/trim them, write them through `add_refinement_region()`, and return FID/name/spacing mappings. +- `replace_refinement_regions(geom_number, regions, expected_existing_names=..., ...)` - Atomically replace or remove the complete HDF refinement-region collection, with an optional optimistic-concurrency guard. - `get_refinement_regions(geom_number)` - Read refinement-region FID, name, and spacing values from a compiled geometry HDF. - `set_refinement_region_spacing(geom_number, spacing_dx, ...)` - Update spacing for one or more existing refinement regions. - `set_refinement_region_name(geom_number, new_name, ...)` - Rename an existing refinement region. @@ -359,6 +367,17 @@ Storage area and 2D flow area geometry parsing and writing. - `get_2d_flow_area_settings(geom_file)` - Read 2D flow area computation settings - `set_2d_flow_area_settings(geom_file, area_name, **settings)` - Write 2D flow area settings (subgrid sampling, composite classification) - `write_2d_flow_area_perimeter(geom_file, area_name, coordinates, ...)` - Write 2D flow area perimeter +- `replace_breaklines(geom_file, flow_area_name, breaklines, expected_existing_names=..., ...)` - Atomically replace the geometry-global breakline collection while preserving supplied near/far spacing, near-repeat, and protection-radius values. + +## MeshRegenerationWorkflow + +Exact RAS Mapper geometry import and legacy mesh-regeneration GUI workflows. + +### Methods + +- `refresh_geometry_hdf_from_text(geom_number=..., geometry_name=..., flow_area_name=..., ras_object=..., ...)` - Transactionally displace one exact geometry HDF, let the explicitly initialized HEC-RAS version rebuild it from task-local `.g##` text, validate the exact 2D perimeter and sibling-HDF isolation, and roll back on failure. This imports geometry features but does not create computation cells. +- `regenerate_mesh(geom_number=..., geometry_name=..., flow_area_name=..., ras_object=..., ...)` - Open/save and validate an already-current exact geometry and compiled mesh. +- `regenerate_mesh_iterative(...)` - Legacy retry workflow; exact geometry selectors are supported and no first-registration fallback is used. ## GeomLevee diff --git a/docs/user-guide/geometry-operations.md b/docs/user-guide/geometry-operations.md index 375975434..79fca4fd5 100644 --- a/docs/user-guide/geometry-operations.md +++ b/docs/user-guide/geometry-operations.md @@ -194,6 +194,132 @@ regional = RasGeo.get_regional_mannings("01", "2D Flow Area") RasGeo.set_base_mannings_table("01", updated_table) ``` +## Rebuilding a copied 2D geometry from text + +For a task-local breakout model, initialize the copied project with the HEC-RAS +version that will perform the work. An existing geometry HDF remains +authoritative in RAS Mapper, so opening and saving it does **not** import an +externally edited `.g##` perimeter or breakline collection. Use the exact, +transactional import workflow instead: + +```python +from pathlib import Path +from ras_commander import ( + GeomMesh, + GeomReferenceFeatures, + GeomStorage, + init_ras_project, +) +from ras_commander.gui.workflows import MeshRegenerationWorkflow + +ras = init_ras_project( + Path(r"C:\tasks\breakout\Model.prj"), + "6.6", + ras_object="new", + load_results_summary=False, +) + +# The caller has already edited only the task-local cloned g03 text. +GeomStorage.set_2d_flow_area_perimeter( + ras.project_folder / "Model.g03", + "Breakout Area", + reduced_domain_polygon, +) +GeomStorage.replace_breaklines( + ras.project_folder / "Model.g03", + "Breakout Area", + retained_and_clipped_breaklines, + expected_existing_names=source_breakline_names, +) + +# Remove parent reference lines outside the reduced domain and replace the +# retained collection in one guarded text mutation. +GeomReferenceFeatures.replace_reference_lines( + ras.project_folder / "Model.g03", + retained_reference_lines, + storage_area="Breakout Area", + expected_existing_names=source_reference_line_names, +) + +refresh = MeshRegenerationWorkflow.refresh_geometry_hdf_from_text( + geom_number="03", + geometry_name="Breakout Geometry", + flow_area_name="Breakout Area", + ras_object=ras, +) +if not refresh.success: + raise refresh.error + +# HDF-only collections are replaced after RAS Mapper imports the text. +GeomMesh.replace_refinement_regions( + "03", + retained_refinement_regions, + expected_existing_names=source_refinement_names, + ras_object=ras, +) + +# The admissible feature envelope is smaller than the new domain: the exact +# compiled perimeter buffered inward by one full base-cell spacing. +containment = GeomMesh.audit_domain_containment( + "03", + mesh_name="Breakout Area", + ras_object=ras, +) +if not containment.ok: + raise ValueError(containment.violations) + +mesh = GeomMesh.generate( + "03", + mesh_name="Breakout Area", + ras_object=ras, + hecras_dir=Path(ras.ras_exe_path).parent, +) +if not mesh.ok: + raise RuntimeError(mesh.error_message) + +GeomMesh.compute_property_tables( + "03", + mesh_name="Breakout Area", + ras_object=ras, +) +``` + +The import defaults to the geometry referenced by the sole current plan when +`geom_number` is omitted. If both number and name are supplied, they must +identify the same unique RAS Mapper layer. Only that HDF is displaced; failure +restores it, non-target geometry HDFs must retain their size and modification +time, and the post-save perimeter must match the text geometrically. The GUI +process tree is supervised from the owned `Ras.exe` PID and no global process +cleanup is performed. Terrain, land-cover, infiltration, and sediment +associations are captured before the import and restored and validated on the +replacement HDF. A missing associated artifact or failed restoration rolls the +transaction back. + +The containment buffer is **inward**, not outward. With a 200-foot base mesh, +the eligible feature envelope is the new perimeter eroded by 200 feet. +Breaklines, refinement regions, and structures must be wholly inside that +envelope before meshing begins. Do not apply this gate to BC lines: inflow and +outflow lines belong on the perimeter and instead need exact 2D-area +association, external-face coverage, endpoint clearance, and mutual-overlap +checks. Reference lines are result-extraction features rather than meshing +constraints; remove any parent line that no longer intersects the reduced +domain so HEC-RAS cannot fail during results processing. + +`HdfBndry.get_breaklines()` exposes `cell_spacing_near`, +`cell_spacing_far`, `near_repeats`, and `protection_radius` so a clipped +collection can preserve the parent BLE meshing controls. Multipart intersections +must be emitted as uniquely named single-part breaklines because the text format +stores one polyline per breakline block. + +Treat the resulting compiled geometry and property tables as reusable run +inputs. A compute may still perform ordinary plan preparation, but do not clear +geometry-preprocessor artifacts for every hydrograph or rating-curve ordinate. +On a real 200,226-cell HEC-RAS 6.6 qualification model, a two-hour Diffusion +Wave solve used about 9 seconds while complete plan preparation and results +processing used about 44 minutes. The operational optimization is therefore to +prepare and certify each immutable reduced geometry once, then reuse it across +the flow series unless geometry-owned inputs change. + ## Hydraulic Tables (HTAB) Extract property tables from preprocessed geometry HDF: diff --git a/ras_commander/__init__.py b/ras_commander/__init__.py index 7bd443376..eb8353d32 100644 --- a/ras_commander/__init__.py +++ b/ras_commander/__init__.py @@ -207,7 +207,8 @@ def agent_guide_path(): 'GeomLateral', 'GeomInlineWeir', 'GeomBridge', 'GeomCulvert', 'GeomCulvertGIS', 'GeomReferenceFeatures', 'GeomBcLines', 'GeomMesh', 'GeomPipeNetwork', 'MeshResult', 'BCConflict', - 'BCFixResult', + 'BCFixResult', 'DomainContainmentResult', + 'DomainContainmentViolation', ) }, **{ @@ -370,7 +371,8 @@ def __getattribute__(self, name): 'GeomInlineWeir', 'GeomBridge', 'GeomCulvert', 'GeomCulvertGIS', 'GeomReferenceFeatures', 'GeomBcLines', 'GeomMesh', 'GeomPipeNetwork', - 'MeshResult', 'BCConflict', 'BCFixResult', + 'MeshResult', 'BCConflict', 'BCFixResult', 'DomainContainmentResult', + 'DomainContainmentViolation', # Deprecated geometry classes (will be removed before v1.0) 'RasGeo', 'RasGeometry', 'RasGeometryUtils', diff --git a/ras_commander/geom/GeomBcLines.py b/ras_commander/geom/GeomBcLines.py index 223a2b5dd..ffb9958c4 100644 --- a/ras_commander/geom/GeomBcLines.py +++ b/ras_commander/geom/GeomBcLines.py @@ -47,7 +47,6 @@ from __future__ import annotations -import shutil from pathlib import Path from typing import Any, Dict, List, Optional, Union @@ -55,6 +54,7 @@ from ..Decorators import log_call from ..LoggingConfig import get_logger +from .GeomParser import GeomParser from .GeomReferenceFeatures import _format_coord_line logger = get_logger(__name__) @@ -285,11 +285,6 @@ def add_bc_lines( if not lines: raise ValueError("lines must contain at least one BC line spec") - # Backup - backup_path = Path(str(geom_path) + ".bak") - shutil.copy2(geom_path, backup_path) - logger.debug(f"Created backup: {backup_path}") - with open(geom_path, "r", encoding="utf-8", errors="ignore", newline="") as f: file_lines = f.readlines() line_ending = _detect_line_ending(file_lines) @@ -372,8 +367,12 @@ def add_bc_lines( file_lines[insert_idx:insert_idx] = new_text_lines - with open(geom_path, "w", encoding="utf-8", newline="") as f: - f.writelines(file_lines) + backup_path = GeomParser.safe_write_geometry( + geom_path, + file_lines, + create_backup=True, + ) + logger.debug("Created backup: %s", backup_path) inserted = [name for name, _, _ in prepared if name not in replaced] logger.info( @@ -433,9 +432,6 @@ def delete_bc_line( if not clean_name: raise ValueError("name is required") - backup_path = Path(str(geom_path) + ".bak") - shutil.copy2(geom_path, backup_path) - with open(geom_path, "r", encoding="utf-8", errors="ignore", newline="") as f: file_lines = f.readlines() @@ -448,8 +444,12 @@ def delete_bc_line( lines_removed = end - start del file_lines[start:end] - with open(geom_path, "w", encoding="utf-8", newline="") as f: - f.writelines(file_lines) + backup_path = GeomParser.safe_write_geometry( + geom_path, + file_lines, + create_backup=True, + ) + logger.debug("Created backup: %s", backup_path) logger.info( "Deleted BC line %s from %s (%d lines)", @@ -508,9 +508,6 @@ def rename_bc_line( if clean_old == clean_new: raise ValueError("old_name and new_name are identical") - backup_path = Path(str(geom_path) + ".bak") - shutil.copy2(geom_path, backup_path) - with open(geom_path, "r", encoding="utf-8", errors="ignore", newline="") as f: file_lines = f.readlines() line_ending = _detect_line_ending(file_lines) @@ -531,8 +528,12 @@ def rename_bc_line( # ending. file_lines[start] = f"{_BC_NAME_KEY}{clean_new:<40s}{line_ending}" - with open(geom_path, "w", encoding="utf-8", newline="") as f: - f.writelines(file_lines) + backup_path = GeomParser.safe_write_geometry( + geom_path, + file_lines, + create_backup=True, + ) + logger.debug("Created backup: %s", backup_path) logger.info( "Renamed BC line %s -> %s in %s", diff --git a/ras_commander/geom/GeomMesh.py b/ras_commander/geom/GeomMesh.py index f67a2ff43..e2ae2b1b1 100644 --- a/ras_commander/geom/GeomMesh.py +++ b/ras_commander/geom/GeomMesh.py @@ -22,23 +22,28 @@ Only overwrite if the caller explicitly passes bl_spacing_near/bl_spacing_far. 2. **Existing HDF validation** — Require a current .g##.hdf compiled by HEC-RAS/Ras.exe before any mesh-generation work begins. -3. **Text → HDF sync** — Sync per-breakline spacing from text into the HDF so +3. **Containment gate** — Erode the exact compiled 2D perimeter inward by one + base-cell spacing. Require every associated breakline, refinement region, + and structure to be wholly covered by that admissible polygon. Boundary + condition lines are checked separately because they belong on the perimeter. +4. **Text → HDF sync** — Sync per-breakline spacing from text into the HDF so RegenerateMeshPoints (which reads HDF, not text) uses correct values. -4. **Load .NET geometry** — RASGeometry(hdf_path) → D2FlowArea → perimeter, +5. **Load .NET geometry** — RASGeometry(hdf_path) → D2FlowArea → perimeter, breaklines (merged BreakLines + Regions + Structures via _build_breaklines). -5. **Generate seeds** — Primary: RegenerateMeshPoints (private .NET method via +6. **Generate seeds** — Primary: RegenerateMeshPoints (private .NET method via reflection) produces breakline-aware seeds. Fallback: PointGenerator. GeneratePoints(perim, cell_size) for base-grid seeds. -6. **Fix loop** (matches TryAutoFix tier ordering): - - Tier 0: Pre-flight removal of short perimeter segments +7. **Fix loop** (matches TryAutoFix tier ordering): + - Tier 0: Diagnose short perimeter segments while preserving the exact + authored perimeter - Tier 1: DuplicatePoints → remove duplicate seed points - Tier 2 first: MaxFacesPerCellExceeded → add midpoint seeds - Tier 3: FacePerimeterConnectionError → remove bad perimeter vertices - Tier 4: Ratio escalation [0.05 → 0.10 → 0.15 → 0.25] - Tier 5: Douglas-Peucker perimeter simplification (last resort) -7. **Extract cell centers** — geom.Save() + h5py read (fast), or .NET Cell(i) +8. **Extract cell centers** — geom.Save() + h5py read (fast), or .NET Cell(i) iteration (slow fallback). -8. **Write .g01 text** — _patch_text_seeds() writes cell centers as the sole +9. **Write .g01 text** — _patch_text_seeds() writes cell centers as the sole persistent output. _set_point_generation_data() updates the seed count header. Requires: @@ -54,15 +59,14 @@ from __future__ import annotations -import logging import os import platform import shutil import sys -from dataclasses import dataclass, field from numbers import Number from pathlib import Path -from typing import Any, Callable, List, Optional, Tuple, Union +from typing import Any, Callable, List, Mapping, Optional, Sequence, Tuple, TYPE_CHECKING, Union +from uuid import uuid4 from ..Decorators import log_call from .._gdal_runtime import ( @@ -78,7 +82,18 @@ resolve_association_attr_path as _shared_resolve_association_attr_path, ) from ..LoggingConfig import get_logger -from .GeomMeshDataclasses import BCConflict, BCFixResult, MeshResult +from .GeomMeshDataclasses import ( + BCConflict, + BCFixResult, + DomainContainmentResult, + DomainContainmentViolation, + MeshResult, +) + +if TYPE_CHECKING: + import numpy + from RasMapperLib import PointMs + from ..RasPrj import RasPrj logger = get_logger(__name__) @@ -454,6 +469,26 @@ def _remove_seed_indexes(seeds_pm, indexes: set[int], ns: dict): +def _seed_indexes_outside_perimeter(seeds_pm, perimeter) -> set[int]: + """Return seed indexes outside an exact .NET perimeter polygon.""" + from shapely.geometry import Point, Polygon + + perimeter_coords = [ + (float(perimeter.PointM(index).X), float(perimeter.PointM(index).Y)) + for index in range(perimeter.Count) + ] + polygon = Polygon(perimeter_coords) + if polygon.is_empty or not polygon.is_valid: + return set() + return { + index + for index in range(seeds_pm.Count) + if not polygon.covers( + Point(float(seeds_pm[index].X), float(seeds_pm[index].Y)) + ) + } + + def _autofix_max_faces(mesh, seeds_as_list: list, ns: dict) -> Tuple[list, int, list]: """Add midpoints of longest 2 faces (with no internal points) per cell exceeding MAX_FACES. @@ -880,9 +915,9 @@ def _set_breakline_spacing_impl( if breakline_name is not None: matches = [ - i for i, l in enumerate(lines) - if l.startswith("BreakLine Name=") - and l.split("=", 1)[1].strip() == breakline_name + i for i, line in enumerate(lines) + if line.startswith("BreakLine Name=") + and line.split("=", 1)[1].strip() == breakline_name ] if len(matches) > 1: raise ValueError( @@ -1496,8 +1531,6 @@ def _patch_text_seeds( When mesh_name is given, only the block under the matching ``Storage Area=`` header is replaced (multi-area safe). """ - import numpy as _np - lines = geom_text_path.read_text(encoding="utf-8", errors="replace").splitlines( keepends=True ) @@ -1766,11 +1799,411 @@ def _mesh_metadata_is_meaningful( ) +def _text_flow_area_perimeter( + geom_text_path: Path, + mesh_name: str, +) -> Optional["numpy.ndarray"]: + """Read one exact 2D-area perimeter without requiring GeoPandas.""" + import numpy as np + from .GeomStorage import GeomStorage + + lines = geom_text_path.read_text(encoding="utf-8", errors="replace").splitlines( + keepends=True + ) + block = GeomStorage._find_storage_area_block(lines, mesh_name) + if block is None: + return None + _start, _end, block_lines = block + info = GeomStorage._inspect_storage_area_block(block_lines) + if not info["is_2d"] or info["surface_line_idx"] is None: + return None + return np.asarray( + GeomStorage._parse_surface_line_coords(block_lines, info), + dtype=float, + ) + + +def _hdf_flow_area_perimeter( + hdf_path: Path, + mesh_name: str, +) -> "numpy.ndarray": + """Read one collection-level HDF perimeter before or after meshing.""" + import h5py + import numpy as np + + with h5py.File(str(hdf_path), "r") as hdf: + group = hdf["Geometry/2D Flow Areas"] + attributes = group["Attributes"][()] + names = [] + for value in attributes["Name"]: + if isinstance(value, (bytes, np.bytes_)): + names.append( + bytes(value).decode("utf-8", errors="replace").rstrip("\x00").strip() + ) + else: + names.append(str(value).rstrip("\x00").strip()) + matches = [index for index, name in enumerate(names) if name == mesh_name] + if len(matches) != 1: + raise ValueError( + f"2D flow area {mesh_name!r} resolved to {len(matches)} HDF rows" + ) + point_start, point_count = ( + int(value) for value in group["Polygon Info"][matches[0], :2] + ) + return np.asarray( + group["Polygon Points"][point_start : point_start + point_count], + dtype=float, + ) + + +def _decode_hdf_text(value: Any) -> str: + """Decode one HDF compound-field value without guessing an encoding.""" + if isinstance(value, (bytes, bytearray)): + return bytes(value).decode("utf-8", errors="replace").rstrip("\x00").strip() + try: + import numpy as np + + if isinstance(value, np.bytes_): + return bytes(value).decode("utf-8", errors="replace").rstrip("\x00").strip() + except ImportError: # pragma: no cover - NumPy is a required dependency + pass + return str(value).rstrip("\x00").strip() + + +def _hdf_attribute_text( + attributes: Any, + index: int, + candidates: Sequence[str], +) -> str: + """Read the first available text field from one compound HDF row.""" + fields = attributes.dtype.names or () + for candidate in candidates: + if candidate in fields: + return _decode_hdf_text(attributes[candidate][index]) + return "" + + +def _hdf_polyline_features( + hdf: Any, + group_path: str, + *, + info_name: str = "Polyline Info", + parts_name: str = "Polyline Parts", + points_name: str = "Polyline Points", +) -> list[dict[str, Any]]: + """Read a native HEC-RAS polyline collection with strict schema checks.""" + from shapely.geometry import LineString, MultiLineString + + if group_path not in hdf: + return [] + group = hdf[group_path] + required = ("Attributes", info_name, points_name) + if not any(name in group for name in required): + if group_path == "Geometry/Structures": + count_attributes = ( + "Bridge/Culvert Count", + "Connection Count", + "Inline Structure Count", + "Lateral Structure Count", + ) + observed_counts = [ + int(group.attrs[name]) + for name in count_attributes + if name in group.attrs + ] + if observed_counts and all(count == 0 for count in observed_counts): + return [] + missing = [name for name in required if name not in group] + if missing: + raise RuntimeError( + f"{group_path} is missing required dataset(s): {', '.join(missing)}" + ) + attributes = group["Attributes"][()] + info = group[info_name][()] + points = group[points_name][()] + if len(attributes) != len(info): + raise RuntimeError( + f"{group_path} Attributes rows ({len(attributes)}) do not match " + f"{info_name} rows ({len(info)})" + ) + parts = group[parts_name][()] if parts_name in group else None + features: list[dict[str, Any]] = [] + for index, row in enumerate(info): + if len(row) < 2: + raise RuntimeError(f"{group_path}/{info_name} row {index} is malformed") + point_start, point_count = int(row[0]), int(row[1]) + part_start = int(row[2]) if len(row) > 2 else 0 + part_count = int(row[3]) if len(row) > 3 else 1 + feature_points = points[point_start : point_start + point_count] + geometry = None + parse_error = "" + if len(feature_points) < 2: + parse_error = "fewer_than_two_points" + elif part_count <= 1: + geometry = LineString(feature_points) + else: + if parts is None: + raise RuntimeError( + f"{group_path}/{parts_name} is required for multipart feature {index}" + ) + line_parts = [] + for part_row in parts[part_start : part_start + part_count]: + raw_start, relative_count = int(part_row[0]), int(part_row[1]) + relative_start = ( + raw_start - point_start + if point_start <= raw_start < point_start + point_count + else raw_start + ) + part_points = feature_points[ + relative_start : relative_start + relative_count + ] + if len(part_points) >= 2: + line_parts.append(LineString(part_points)) + if line_parts: + geometry = ( + line_parts[0] + if len(line_parts) == 1 + else MultiLineString(line_parts) + ) + else: + parse_error = "multipart_feature_has_no_valid_parts" + features.append( + { + "index": index, + "name": _hdf_attribute_text( + attributes, + index, + ("Name", "Structure Name", "Connection Name"), + ), + "structure_type": _hdf_attribute_text( + attributes, + index, + ("Type", "Structure Type", "Connection Type"), + ), + "associations": [ + _hdf_attribute_text(attributes, index, (field,)) + for field in ("SA-2D", "US SA/2D", "DS SA/2D", "2D Area Name") + if field in (attributes.dtype.names or ()) + ], + "geometry": geometry, + "parse_error": parse_error, + } + ) + return features + + +def _hdf_polygon_features( + hdf: Any, + group_path: str, +) -> list[dict[str, Any]]: + """Read a native HEC-RAS polygon collection with strict schema checks.""" + from shapely.geometry import MultiPolygon, Polygon + + if group_path not in hdf: + return [] + group = hdf[group_path] + required = ("Attributes", "Polygon Info", "Polygon Points") + missing = [name for name in required if name not in group] + if missing: + raise RuntimeError( + f"{group_path} is missing required dataset(s): {', '.join(missing)}" + ) + attributes = group["Attributes"][()] + info = group["Polygon Info"][()] + points = group["Polygon Points"][()] + if len(attributes) != len(info): + raise RuntimeError( + f"{group_path} Attributes rows ({len(attributes)}) do not match " + f"Polygon Info rows ({len(info)})" + ) + parts = group["Polygon Parts"][()] if "Polygon Parts" in group else None + features: list[dict[str, Any]] = [] + for index, row in enumerate(info): + if len(row) < 2: + raise RuntimeError(f"{group_path}/Polygon Info row {index} is malformed") + point_start, point_count = int(row[0]), int(row[1]) + part_start = int(row[2]) if len(row) > 2 else 0 + part_count = int(row[3]) if len(row) > 3 else 1 + feature_points = points[point_start : point_start + point_count] + geometry = None + parse_error = "" + if len(feature_points) < 4: + parse_error = "fewer_than_four_polygon_points" + elif part_count <= 1: + geometry = Polygon(feature_points) + else: + if parts is None: + raise RuntimeError( + f"{group_path}/Polygon Parts is required for multipart feature {index}" + ) + polygons = [] + for part_row in parts[part_start : part_start + part_count]: + raw_start, relative_count = int(part_row[0]), int(part_row[1]) + relative_start = ( + raw_start - point_start + if point_start <= raw_start < point_start + point_count + else raw_start + ) + part_points = feature_points[ + relative_start : relative_start + relative_count + ] + if len(part_points) >= 4: + polygons.append(Polygon(part_points)) + if polygons: + geometry = polygons[0] if len(polygons) == 1 else MultiPolygon(polygons) + else: + parse_error = "multipart_feature_has_no_valid_parts" + features.append( + { + "index": index, + "name": _hdf_attribute_text(attributes, index, ("Name",)), + "structure_type": "", + "associations": [], + "geometry": geometry, + "parse_error": parse_error, + } + ) + return features + + +def _audit_domain_containment_hdf( + hdf_path: Path, + mesh_name: str, + base_cell_spacing: float, +) -> DomainContainmentResult: + """Audit mesh-owned features against a one-cell inward domain buffer.""" + import h5py + from shapely.geometry import Polygon + + spacing = float(base_cell_spacing) + if not spacing > 0: + raise ValueError("base_cell_spacing must be positive") + perimeter = Polygon(_hdf_flow_area_perimeter(hdf_path, mesh_name)) + if perimeter.is_empty or not perimeter.is_valid or perimeter.area <= 0: + raise RuntimeError(f"2D flow area {mesh_name!r} has an invalid HDF perimeter") + admissible = perimeter.buffer(-spacing) + if admissible.is_empty or not admissible.is_valid or admissible.area <= 0: + raise RuntimeError( + f"Inward buffer of {spacing:g} project units empties or invalidates " + f"2D flow area {mesh_name!r}" + ) + + collections: list[tuple[str, list[dict[str, Any]]]] = [] + with h5py.File(str(hdf_path), "r") as hdf: + collections.append( + ( + "breakline", + _hdf_polyline_features(hdf, "Geometry/2D Flow Area Break Lines"), + ) + ) + collections.append( + ( + "refinement_region", + _hdf_polygon_features( + hdf, + "Geometry/2D Flow Area Refinement Regions", + ), + ) + ) + structures = _hdf_polyline_features( + hdf, + "Geometry/Structures", + info_name="Centerline Info", + parts_name="Centerline Parts", + points_name="Centerline Points", + ) + relevant_structures = [] + for feature in structures: + associations = [value for value in feature["associations"] if value] + if associations and mesh_name not in associations: + continue + relevant_structures.append(feature) + collections.append(("structure", relevant_structures)) + + checked_counts: dict[str, int] = {} + violations: list[DomainContainmentViolation] = [] + for feature_type, features in collections: + checked_counts[feature_type] = len(features) + for feature in features: + geometry = feature["geometry"] + reason = feature["parse_error"] + if geometry is not None and not reason: + if geometry.is_empty or not geometry.is_valid: + reason = "invalid_or_empty_geometry" + elif not admissible.covers(geometry): + reason = "outside_one_cell_inward_buffer" + if not reason: + continue + outside_length = 0.0 + outside_area = 0.0 + geometry_type = "" + if geometry is not None: + geometry_type = geometry.geom_type + if reason == "outside_one_cell_inward_buffer": + outside = geometry.difference(admissible) + outside_length = float(outside.length) + outside_area = float(outside.area) + violations.append( + DomainContainmentViolation( + feature_type=feature_type, + feature_name=feature["name"], + feature_index=int(feature["index"]), + geometry_type=geometry_type, + reason=reason, + outside_length=outside_length, + outside_area=outside_area, + structure_type=feature["structure_type"], + ) + ) + + return DomainContainmentResult( + mesh_name=mesh_name, + geom_hdf_path=str(hdf_path), + base_cell_spacing=spacing, + inward_buffer_distance=spacing, + admissible_geometry_type=admissible.geom_type, + checked_counts=checked_counts, + violations=violations, + ) + + +def _closed_rings_match(left, right) -> bool: + """Compare closed XY rings independent of start vertex and orientation.""" + import numpy as np + + left = np.asarray(left, dtype=float) + right = np.asarray(right, dtype=float) + if len(left) > 1 and np.allclose(left[0], left[-1], rtol=0.0, atol=1e-9): + left = left[:-1] + if len(right) > 1 and np.allclose(right[0], right[-1], rtol=0.0, atol=1e-9): + right = right[:-1] + if left.shape != right.shape or left.ndim != 2 or left.shape[1] != 2: + return False + if len(left) == 0: + return False + + scale = max(float(np.ptp(left[:, 0])), float(np.ptp(left[:, 1])), 1.0) + tolerance = max(1e-6, scale * 1e-9) + first = left[0] + candidates = np.flatnonzero( + np.max(np.abs(right - first), axis=1) <= tolerance + ) + for index in candidates: + forward = np.roll(right, -int(index), axis=0) + if np.allclose(left, forward, rtol=0.0, atol=tolerance): + return True + reverse = np.roll(right[::-1], int(index) + 1, axis=0) + if np.allclose(left, reverse, rtol=0.0, atol=tolerance): + return True + return False + + def _mesh_hdf_consistency_issues( geom_text_path: Path, hdf_path: Path, *, mesh_name: str | None = None, + ignore_seed_count: bool = False, ) -> list[str]: """Return content mismatches between geometry text seeds and HDF metadata.""" text_metadata = _read_mesh_metadata_from_text(geom_text_path) @@ -1799,7 +2232,8 @@ def _mesh_hdf_consistency_issues( hdf_cell_count = hdf_area.get("cell_count") hdf_count_source = hdf_area.get("cell_count_source") if ( - isinstance(text_seed_count, int) + not ignore_seed_count + and isinstance(text_seed_count, int) and isinstance(hdf_cell_count, int) and hdf_count_source == "Attributes/Cell Count" and text_seed_count != hdf_cell_count @@ -1826,6 +2260,18 @@ def _mesh_hdf_consistency_issues( f"but HDF Spacing {axis}={float(hdf_spacing):g}" ) + try: + text_perimeter = _text_flow_area_perimeter(geom_text_path, area_name) + if text_perimeter is not None: + hdf_perimeter = _hdf_flow_area_perimeter(hdf_path, area_name) + if not _closed_rings_match(text_perimeter, hdf_perimeter): + issues.append( + f"{area_name}: text perimeter ({len(text_perimeter)} points) " + f"does not match HDF perimeter ({len(hdf_perimeter)} points)" + ) + except (OSError, KeyError, TypeError, ValueError) as exc: + issues.append(f"{area_name}: perimeter comparison failed; {exc}") + # Breakline geometry. HEC-RAS stores breaklines as a flat, global list (not # per-area) in both the .g01 text and the compiled HDF. GeomMesh.generate's .NET # RegenerateMeshPoints seeds the refined corridor from the *HDF* breaklines, so a @@ -1888,6 +2334,7 @@ def _ensure_hdf( ras_object=None, mesh_name: str | None = None, recompile_via_rasexe: bool = False, + ignore_seed_count: bool = False, ) -> Path: """Return an existing compiled HDF for *geom_text_path*. @@ -1914,6 +2361,7 @@ def _ensure_hdf( geom_text_path, hdf_path, mesh_name=mesh_name, + ignore_seed_count=ignore_seed_count, ) except Exception as exc: consistency_issues = [f"unreadable; {exc}"] @@ -1930,6 +2378,7 @@ def _ensure_hdf( geom_text_path, hdf_path, mesh_name=mesh_name, + ignore_seed_count=ignore_seed_count, ) if not consistency_issues: return hdf_path @@ -2303,6 +2752,75 @@ def setup_gdal_bridge( return True + @staticmethod + @log_call + def audit_domain_containment( + geom_number: Union[str, Number, Path], + mesh_name: Optional[str] = None, + mesh_index: int = 0, + cell_size: Optional[float] = None, + ras_object=None, + recompile_via_rasexe: bool = False, + ) -> DomainContainmentResult: + """Verify mesh-owned features stay one base cell inside the 2D perimeter. + + The admissible geometry is the exact compiled 2D perimeter buffered + *inward* by one base mesh-cell spacing. Every breakline, refinement + region, and structure associated with the selected 2D area must be + wholly covered by that eroded polygon. Boundary-condition lines are + intentionally excluded because external BC lines are authored at the + perimeter and require a separate association audit. + + This read-only check is cross-platform and does not load RasMapperLib. + A missing, stale, malformed, or ambiguous geometry HDF fails closed. + """ + geom_path = _resolve_geom_text_path(geom_number, ras_object) + hdf_path = _ensure_hdf( + geom_path, + require_current=True, + ras_object=ras_object, + mesh_name=mesh_name, + recompile_via_rasexe=recompile_via_rasexe, + ignore_seed_count=True, + ) + hdf_metadata = _read_mesh_metadata_from_hdf(hdf_path) + area_names = list(hdf_metadata) + if mesh_name is None: + if not isinstance(mesh_index, int) or isinstance(mesh_index, bool): + raise TypeError("mesh_index must be an integer") + if mesh_index < 0 or mesh_index >= len(area_names): + raise IndexError( + f"mesh_index {mesh_index} is outside the {len(area_names)} " + "compiled 2D flow area(s)" + ) + mesh_name = area_names[mesh_index] + elif mesh_name not in hdf_metadata: + raise ValueError( + f"2D flow area {mesh_name!r} is not present exactly once in " + f"{hdf_path.name}" + ) + + if cell_size is None: + cell_size = _read_cell_size_from_text( + geom_path, + mesh_name=mesh_name, + mesh_index=mesh_index, + ) + if cell_size is None: + hdf_spacing = hdf_metadata[mesh_name].get("spacing_dx") + if isinstance(hdf_spacing, (int, float)) and hdf_spacing > 0: + cell_size = float(hdf_spacing) + if cell_size is None: + raise ValueError( + f"Base mesh spacing is unavailable for 2D flow area {mesh_name!r}" + ) + cell_size = _normalize_positive_value(cell_size, "cell_size") + return _audit_domain_containment_hdf( + hdf_path, + mesh_name, + cell_size, + ) + @staticmethod @log_call def set_breakline_spacing( @@ -2792,6 +3310,183 @@ def set_refinement_region_spacing( f"→ {hdf_path.name}" ) + @staticmethod + @log_call + def replace_refinement_regions( + geom_number: Union[str, Number, Path], + regions: Sequence[Mapping[str, Any]], + *, + expected_existing_names: Optional[Sequence[str]] = None, + hecras_dir: Optional[Union[str, Path]] = None, + ras_object=None, + create_backup: bool = True, + ) -> Optional[Path]: + """Atomically replace the complete HDF refinement-region collection. + + Each region mapping requires ``polygon`` (or ``geometry``) and + ``spacing_dx``; ``spacing_dy`` defaults to ``spacing_dx`` and ``name`` + defaults to an empty string. The replacement is built and validated in + a same-directory temporary HDF before it is promoted over the original. + + ``expected_existing_names`` is an optimistic-concurrency guard. Names + are compared in HDF/FID order and may contain duplicates or empty + strings because HEC-RAS permits both. Passing an empty ``regions`` + sequence removes the refinement-region group. + """ + import h5py + import numpy as np + + from .GeomParser import GeomParser + + geom_text_path = _resolve_geom_text_path(geom_number, ras_object) + hdf_path = _ensure_hdf( + geom_text_path, + hecras_dir=hecras_dir, + ras_object=ras_object, + # Refinement-region replacement does not consume computation + # cells. A freshly imported, intentionally unmeshed HDF is valid + # as long as perimeter and breakline identity are current. + ignore_seed_count=True, + ) + rr_group_key = "Geometry/2D Flow Area Refinement Regions" + attr_key = f"{rr_group_key}/Attributes" + info_key = f"{rr_group_key}/Polygon Info" + parts_key = f"{rr_group_key}/Polygon Parts" + points_key = f"{rr_group_key}/Polygon Points" + + with h5py.File(str(hdf_path), "r") as hf: + if attr_key in hf: + existing_names = [ + row.decode("utf-8", errors="replace").strip() + for row in hf[attr_key]["Name"] + ] + else: + existing_names = [] + if expected_existing_names is not None: + expected_existing = [str(name) for name in expected_existing_names] + if existing_names != expected_existing: + raise ValueError( + "Existing refinement-region collection changed before replacement: " + f"expected {expected_existing!r}, found {existing_names!r}" + ) + + normalized: List[dict] = [] + for index, raw in enumerate(regions): + if not isinstance(raw, Mapping): + raise TypeError(f"regions[{index}] must be a mapping") + polygon = raw.get("polygon", raw.get("geometry")) + if polygon is None: + raise ValueError(f"regions[{index}] is missing polygon or geometry") + coords = _normalise_polygon_coords(polygon) + if len(coords) < 3: + raise ValueError( + f"regions[{index}] polygon must have at least 3 vertices" + ) + if not np.isfinite(coords).all(): + raise ValueError(f"regions[{index}] polygon has non-finite coordinates") + if not np.allclose(coords[0], coords[-1]): + coords = np.vstack([coords, coords[:1]]) + + spacing_dx = float(raw.get("spacing_dx")) + spacing_dy = float(raw.get("spacing_dy", spacing_dx)) + if ( + not np.isfinite(spacing_dx) + or not np.isfinite(spacing_dy) + or spacing_dx <= 0 + or spacing_dy <= 0 + ): + raise ValueError( + f"regions[{index}] spacing must be finite and positive" + ) + normalized.append( + { + "name": _truncate_ras_name(str(raw.get("name", ""))), + "coords": np.asarray(coords, dtype=np.float64), + "spacing_dx": spacing_dx, + "spacing_dy": spacing_dy, + } + ) + + attr_dtype = np.dtype( + [("Name", "S32"), ("Spacing dx", " perimeter, breaklines 5. Generate seeds: RegenerateMeshPoints (primary) or GeneratePoints 6. Fix loop (matches TryAutoFix tier ordering): - - Tier 0: Remove short perimeter segments (pre-flight) + - Tier 0: Diagnose short perimeter segments while preserving + the exact authored perimeter - Tier 1: DuplicatePoints -> remove duplicate seed points - Tier 2: MaxFaces -> add midpoint seeds (before ratio escalation) - Tier 3: Perimeter errors -> remove bad vertices @@ -3640,8 +4336,6 @@ def generate( MeshResult with status, cell_count, face_count, fixes_applied, and the compiled geometry HDF path on success. """ - _load_dlls(hecras_dir) - ns = _imports() geom_path = _resolve_geom_text_path(geom_number, ras_object) result = MeshResult( @@ -3684,6 +4378,10 @@ def generate( ras_object=ras_object, mesh_name=mesh_name, recompile_via_rasexe=recompile_via_rasexe, + # Seed-count mismatch is expected here: this method replaces + # computation points. Perimeter and breakline identity remain + # mandatory content-current gates. + ignore_seed_count=True, ) # ── Step 1a: Auto-detect cell size from text if not provided ── @@ -3703,6 +4401,38 @@ def generate( f"Auto-detected cell size {cell_size} from HDF Spacing dx" ) + # ── Step 1b: Fail-closed mesh-feature containment gate ─── + # Breaklines, refinement regions, and SA/2D structures must + # remain one full base cell inside the exact new perimeter. + # External boundary-condition lines are audited separately. + domain_containment = _audit_domain_containment_hdf( + hdf_path, + mesh_name + or list(_read_mesh_metadata_from_hdf(hdf_path))[mesh_index], + float(cell_size), + ) + result.domain_containment = domain_containment + result.mesh_name = domain_containment.mesh_name + mesh_name = domain_containment.mesh_name + result.geom_hdf_path = str(hdf_path) + if not domain_containment: + preview = ", ".join( + f"{item.feature_type}:{item.feature_name or item.feature_index}" + for item in domain_containment.violations[:5] + ) + suffix = "..." if len(domain_containment.violations) > 5 else "" + result.error_message = ( + f"Pre-mesh domain containment failed for " + f"{len(domain_containment.violations)} feature(s): " + f"{preview}{suffix}" + ) + return result + + # Native dependencies are loaded only after the read-only spatial + # gate succeeds, so an invalid breakout cannot start meshing. + _load_dlls(hecras_dir) + ns = _imports() + # Only modify .g01 text breakline properties if explicitly provided. # Otherwise the geometry's existing per-breakline values are # the source of truth — don't overwrite them with defaults. @@ -3719,7 +4449,7 @@ def generate( _log_summary=False, ) - # ── Step 1b: Sync text → HDF ──────────────────────────────── + # ── Step 1c: Sync text → HDF ──────────────────────────────── # RegenerateMeshPoints reads spacing from the HDF, not from # .g01 text. Sync current cell size and per-breakline values from # text into the existing HDF workspace. @@ -3762,7 +4492,6 @@ def generate( # ── Step 4: Generate seeds via .NET ────────────────────────── # Always try RegenerateMeshPoints first — it uses the correct # grid origin from the HDF regardless of whether breaklines exist. - PointGenerator = ns["PointGenerator"] net_seeds_ok = False try: seeds_pm = _generate_seeds_via_net(str(hdf_path), ns, fid=fid) @@ -3793,17 +4522,19 @@ def generate( # Tier 0: Pre-simplify short perimeter segments pre_n = current_perim.Count - current_perim = _remove_short_perimeter_segments( + repaired_perim = _remove_short_perimeter_segments( current_perim, cell_size * min_face_length_ratio, ns ) - post_n = current_perim.Count + post_n = repaired_perim.Count if post_n < pre_n: - fix_msg = f"Tier0:short_seg_removal(-{pre_n - post_n})" - result.fixes_applied.append(fix_msg) - logger.debug(f"[{mesh_name}] Fix applied: {fix_msg}") - current_seeds_pm = _reseed_after_perimeter_fix( - text_path, hdf_path, current_perim, - cell_size, fid, mesh_name, ns, hecras_dir, + # A repaired in-memory perimeter would diverge from both the + # authored text and the HDF collection. Try the exact authored + # perimeter first; if HEC-RAS rejects it, report the need for a + # separately compiled geometry edit instead of silently + # changing the domain boundary. + logger.warning( + f"[{mesh_name}] Exact perimeter contains {pre_n - post_n} " + "short segment(s); preserving the authored perimeter" ) ratio_idx = 0 @@ -3826,6 +4557,10 @@ def generate( duplicate_points_val = int(MeshStatus.DuplicatePoints) except AttributeError: duplicate_points_val = None + try: + points_outside_val = int(MeshStatus.PointsOutsidePerimeter) + except AttributeError: + points_outside_val = None for iteration in range(max_iterations): ratio = ratios[min(ratio_idx, len(ratios) - 1)] @@ -3969,6 +4704,35 @@ def generate( ) break + # Remove only generated seed points that HEC-RAS identifies as + # outside. Never respond by mutating the authored perimeter. + if points_outside_val is not None and state_val == points_outside_val: + bad_indexes = _bad_seed_indexes(mesh, current_seeds_pm.Count) + if not bad_indexes: + bad_indexes = _seed_indexes_outside_perimeter( + current_seeds_pm, + current_perim, + ) + if bad_indexes: + current_seeds_pm, n_removed = _remove_seed_indexes( + current_seeds_pm, + bad_indexes, + ns, + ) + if n_removed > 0: + fix_msg = ( + "PointsOutsidePerimeter:seed-removal" + f"(-{n_removed}pts)" + ) + result.fixes_applied.append(fix_msg) + logger.debug(f"[{mesh_name}] Fix applied: {fix_msg}") + continue + result.error_message = ( + "Mesh reported PointsOutsidePerimeter, but no outside " + "generated seed points could be identified." + ) + break + # MaxFaces → add midpoint seeds at current ratio if ( state_val == max_faces_val diff --git a/ras_commander/geom/GeomMeshDataclasses.py b/ras_commander/geom/GeomMeshDataclasses.py index 892640ed3..b6c6b2bf7 100644 --- a/ras_commander/geom/GeomMeshDataclasses.py +++ b/ras_commander/geom/GeomMeshDataclasses.py @@ -8,7 +8,68 @@ from __future__ import annotations from dataclasses import dataclass, field -from typing import List, Optional +from typing import Any, Dict, List, Optional + + +@dataclass(frozen=True) +class DomainContainmentViolation: + """One mesh-owned feature that is not inside the eroded 2D domain.""" + + feature_type: str + feature_name: str + feature_index: int + geometry_type: str + reason: str + outside_length: float = 0.0 + outside_area: float = 0.0 + structure_type: str = "" + + def to_dict(self) -> Dict[str, Any]: + """Return machine-readable violation evidence.""" + return { + "feature_type": self.feature_type, + "feature_name": self.feature_name, + "feature_index": self.feature_index, + "geometry_type": self.geometry_type, + "reason": self.reason, + "outside_length": self.outside_length, + "outside_area": self.outside_area, + "structure_type": self.structure_type, + } + + +@dataclass(frozen=True) +class DomainContainmentResult: + """Pre-mesh containment evidence for one 2D flow area.""" + + mesh_name: str + geom_hdf_path: str + base_cell_spacing: float + inward_buffer_distance: float + admissible_geometry_type: str + checked_counts: Dict[str, int] + violations: List[DomainContainmentViolation] + + @property + def ok(self) -> bool: + return not self.violations + + def __bool__(self) -> bool: + return self.ok + + def to_dict(self) -> Dict[str, Any]: + """Return machine-readable audit evidence without geometry payloads.""" + return { + "mesh_name": self.mesh_name, + "geom_hdf_path": self.geom_hdf_path, + "base_cell_spacing": self.base_cell_spacing, + "inward_buffer_distance": self.inward_buffer_distance, + "admissible_geometry_type": self.admissible_geometry_type, + "checked_counts": dict(self.checked_counts), + "violation_count": len(self.violations), + "ok": self.ok, + "violations": [item.to_dict() for item in self.violations], + } @dataclass @@ -25,6 +86,7 @@ class MeshResult: error_message: str = "" geom_text_path: str = "" geom_hdf_path: str = "" + domain_containment: Optional[DomainContainmentResult] = None @property def ok(self) -> bool: diff --git a/ras_commander/geom/GeomParser.py b/ras_commander/geom/GeomParser.py index 9addee563..f63716007 100644 --- a/ras_commander/geom/GeomParser.py +++ b/ras_commander/geom/GeomParser.py @@ -36,7 +36,7 @@ import re from pathlib import Path -from typing import List, Optional, Tuple, Dict, Any, Union +from typing import List, Optional, Tuple, Union from datetime import datetime from ..LoggingConfig import get_logger @@ -567,7 +567,7 @@ def safe_write_geometry(geom_file: Path, logger.debug(f"Created backup: {backup_path}") # Step 2: Write to temp file - with open(temp_path, 'w', encoding='utf-8') as f: + with open(temp_path, 'w', encoding='utf-8', newline='') as f: f.writelines(modified_lines) # Step 3: Basic validation - check temp file has content @@ -576,14 +576,7 @@ def safe_write_geometry(geom_file: Path, # Step 4: Atomic rename temp -> original import os - if os.name == 'nt': # Windows - # Windows doesn't support atomic rename over existing file - # Remove original first, then rename - geom_file.unlink() - temp_path.rename(geom_file) - else: # Unix-like - # Atomic rename - temp_path.rename(geom_file) + os.replace(temp_path, geom_file) logger.debug(f"Successfully wrote geometry file: {geom_file}") return backup_path diff --git a/ras_commander/geom/GeomReferenceFeatures.py b/ras_commander/geom/GeomReferenceFeatures.py index 37ce081e7..8bfba6367 100644 --- a/ras_commander/geom/GeomReferenceFeatures.py +++ b/ras_commander/geom/GeomReferenceFeatures.py @@ -16,16 +16,82 @@ import math import shutil from pathlib import Path -from typing import Any, Callable, List, Optional, Union +from typing import Any, Callable, List, Mapping, Optional, Sequence, Union import numpy as np from ..Decorators import log_call from ..LoggingConfig import get_logger +from .GeomParser import GeomParser logger = get_logger(__name__) +def _reference_line_blocks(file_lines: List[str]) -> List[dict]: + """Return exact reference-line block extents from geometry text lines.""" + blocks: List[dict] = [] + idx = 0 + while idx < len(file_lines): + stripped = file_lines[idx].rstrip("\r\n") + if not stripped.startswith("Reference Line Name="): + idx += 1 + continue + + start = idx + name = stripped.split("=", 1)[1].strip() + storage_area = "" + idx += 1 + while idx < len(file_lines): + current = file_lines[idx].rstrip("\r\n") + if current.startswith("Reference Line Name="): + break + if current.startswith("Reference Line Storage Area="): + storage_area = current.split("=", 1)[1].strip() + idx += 1 + if current.startswith("Reference Line Text Position="): + break + blocks.append( + { + "name": name, + "storage_area": storage_area, + "start": start, + "end": idx, + } + ) + return blocks + + +def _reference_line_insert_index(file_lines: List[str]) -> int: + """Choose the canonical reference-line insertion point.""" + last_bc_line_idx = -1 + first_existing_refline_idx = -1 + first_ic_point_idx = -1 + first_lcmann_idx = -1 + for idx, line in enumerate(file_lines): + stripped = line.rstrip("\r\n") + if stripped.startswith("BC Line Text Position="): + last_bc_line_idx = idx + if ( + stripped.startswith("Reference Line Name=") + and first_existing_refline_idx == -1 + ): + first_existing_refline_idx = idx + if stripped.startswith("IC Point Name=") and first_ic_point_idx == -1: + first_ic_point_idx = idx + if stripped.startswith("LCMann ") and first_lcmann_idx == -1: + first_lcmann_idx = idx + + if first_existing_refline_idx >= 0: + return first_existing_refline_idx + if last_bc_line_idx >= 0: + return last_bc_line_idx + 1 + if first_ic_point_idx >= 0: + return first_ic_point_idx + if first_lcmann_idx >= 0: + return first_lcmann_idx + return len(file_lines) + + def _format_coord_line(values: List[float], width: int = 16) -> str: """Format coordinate values into fixed-width fields, 4 values per line.""" parts = [] @@ -800,11 +866,6 @@ def add_reference_lines( if not lines: raise ValueError("lines must contain at least one reference line") - # Create backup - backup_path = Path(str(geom_file) + ".bak") - shutil.copy2(geom_file, backup_path) - logger.debug(f"Created backup: {backup_path}") - # Read file with CRLF preservation with open(geom_file, "r", encoding="utf-8", errors="ignore", newline="") as f: file_lines = f.readlines() @@ -858,9 +919,11 @@ def add_reference_lines( insert_lines = [block_line + line_ending for block_line in new_blocks] file_lines[insert_idx:insert_idx] = insert_lines - # Write back - with open(geom_file, "w", encoding="utf-8", newline="") as f: - f.writelines(file_lines) + GeomParser.safe_write_geometry( + geom_file, + file_lines, + create_backup=True, + ) logger.debug( f"Inserted {len(lines)} reference line(s) into {geom_file.name} " @@ -868,6 +931,116 @@ def add_reference_lines( ) return len(lines) + @staticmethod + @log_call + def replace_reference_lines( + geom_file: Union[str, Path], + lines: Sequence[Mapping[str, Any]], + storage_area: str, + expected_existing_names: Optional[Sequence[str]] = None, + ) -> dict: + """Atomically replace one 2D area's complete reference-line collection. + + This is the safe mutation for reduced 2D domains: callers can remove + parent reference lines that no longer intersect the retained mesh and + optionally clip the survivors before replacing them. Reference lines + belonging to other 2D areas are preserved byte-for-byte. + + ``expected_existing_names`` is an optional ordered optimistic- + concurrency guard. An empty ``lines`` list removes all reference lines + associated with ``storage_area``. + """ + geom_path = Path(geom_file) + if not geom_path.is_file(): + raise FileNotFoundError(f"Geometry file not found: {geom_path}") + area = str(storage_area).strip() + if not area: + raise ValueError("storage_area must be non-empty") + + prepared: List[tuple[str, List[str]]] = [] + seen: set[str] = set() + for item in lines: + if not isinstance(item, Mapping): + raise ValueError("each entry in lines must be a mapping") + name = str(item.get("name", "")).strip() + if not name: + raise ValueError("each reference line must have a non-empty name") + if name in seen: + raise ValueError(f"duplicate reference line name: {name!r}") + seen.add(name) + if "coordinates" not in item: + raise ValueError( + f"Reference line {name!r} is missing coordinates" + ) + coords = np.asarray(item["coordinates"], dtype=np.float64) + if coords.ndim != 2 or coords.shape[1] != 2 or len(coords) < 2: + raise ValueError( + f"Reference line {name!r} needs at least 2 points as " + "an (N, 2) array" + ) + if not np.isfinite(coords).all(): + raise ValueError( + f"Reference line {name!r} coordinates must be finite" + ) + prepared.append((name, _build_reference_line_block(name, area, coords))) + + with open( + geom_path, + "r", + encoding="utf-8", + errors="ignore", + newline="", + ) as handle: + file_lines = handle.readlines() + blocks = _reference_line_blocks(file_lines) + target_blocks = [block for block in blocks if block["storage_area"] == area] + existing_names = [str(block["name"]) for block in target_blocks] + if ( + expected_existing_names is not None + and existing_names != [str(name) for name in expected_existing_names] + ): + raise ValueError( + f"Reference-line population changed for {area!r}: expected " + f"{list(expected_existing_names)!r}, observed {existing_names!r}" + ) + + for block in sorted(target_blocks, key=lambda value: value["start"], reverse=True): + del file_lines[int(block["start"]):int(block["end"])] + + insert_idx = _reference_line_insert_index(file_lines) + line_ending = "\r\n" if any( + line.endswith("\r\n") for line in file_lines + ) else "\n" + replacement_lines = [ + block_line + line_ending + for _, block in prepared + for block_line in block + ] + file_lines[insert_idx:insert_idx] = replacement_lines + backup_path = GeomParser.safe_write_geometry( + geom_path, + file_lines, + create_backup=True, + ) + inserted_names = [name for name, _ in prepared] + logger.info( + "Replaced %d reference line(s) with %d for %s in %s", + len(existing_names), + len(inserted_names), + area, + geom_path.name, + ) + return { + "geom_file": str(geom_path), + "storage_area": area, + "removed": existing_names, + "inserted": inserted_names, + "removed_count": len(existing_names), + "inserted_count": len(inserted_names), + "insert_index": insert_idx, + "backup_path": str(backup_path), + } + @staticmethod @log_call def add_reference_points( diff --git a/ras_commander/geom/GeomStorage.py b/ras_commander/geom/GeomStorage.py index 8bea5c401..8350d1543 100644 --- a/ras_commander/geom/GeomStorage.py +++ b/ras_commander/geom/GeomStorage.py @@ -17,6 +17,7 @@ - get_2d_flow_area_settings() - Read 2D flow area cell/face property settings - set_2d_flow_area_settings() - Write 2D flow area cell/face property settings - set_breaklines() - Write breakline blocks into a 2D flow area geometry file +- replace_breaklines() - Atomically replace the geometry-wide breakline collection Example Usage: >>> from ras_commander import GeomStorage @@ -39,6 +40,7 @@ ... ) """ +import math from datetime import datetime from pathlib import Path from typing import TYPE_CHECKING, Union, Optional, List, Sequence @@ -1688,6 +1690,8 @@ def _format_breakline_block( coords: List[tuple], cell_size_near: Optional[float] = None, cell_size_far: Optional[float] = None, + near_repeats: int = 0, + protection_radius: int = 0, ) -> List[str]: """Build a complete breakline text block as a list of lines.""" block = [f"BreakLine Name={name}\n"] @@ -1701,12 +1705,233 @@ def _format_breakline_block( if cell_size_far is not None else "BreakLine CellSize Max=\n" ) - block.append("BreakLine Near Repeats=0\n") - block.append("BreakLine Protection Radius=0\n") + block.append(f"BreakLine Near Repeats={near_repeats}\n") + block.append(f"BreakLine Protection Radius={protection_radius}\n") block.append(f"BreakLine Polyline= {len(coords)} \n") block.extend(GeomStorage._format_breakline_coord_lines(coords)) return block + @staticmethod + def _breakline_block_ranges(lines: List[str]) -> List[tuple[int, int]]: + """Return exact ``[start, end)`` ranges for breakline text blocks.""" + ranges: List[tuple[int, int]] = [] + for start_idx, line in enumerate(lines): + if not line.startswith("BreakLine Name="): + continue + + polyline_idx = None + point_count = None + for idx in range(start_idx + 1, len(lines)): + candidate = lines[idx] + if candidate.startswith("BreakLine Name="): + break + if candidate.startswith("BreakLine Polyline="): + try: + point_count = int(candidate.split("=", 1)[1].strip()) + except (IndexError, ValueError) as exc: + raise ValueError( + f"Malformed breakline point count at line {idx + 1}" + ) from exc + polyline_idx = idx + break + + if polyline_idx is None or point_count is None: + raise ValueError( + f"Breakline at line {start_idx + 1} has no polyline record" + ) + coordinate_lines = (point_count + 1) // 2 + end_idx = polyline_idx + 1 + coordinate_lines + if end_idx > len(lines): + raise ValueError( + f"Breakline at line {start_idx + 1} has truncated coordinates" + ) + ranges.append((start_idx, end_idx)) + return ranges + + @staticmethod + def _normalize_breakline_specs(breaklines: Sequence[dict]) -> List[dict]: + """Validate and normalize a complete replacement collection.""" + normalized: List[dict] = [] + names: set[str] = set() + for index, raw in enumerate(breaklines): + if not isinstance(raw, dict): + raise TypeError(f"breaklines[{index}] must be a dict") + name = str(raw.get("name", "")).strip() + GeomStorage._validate_flow_area_name(name) + key = name.casefold() + if key in names: + raise ValueError(f"Duplicate breakline name: {name!r}") + names.add(key) + + raw_coords = raw.get("coords") + if raw_coords is None: + raise ValueError(f"Breakline {name!r} is missing coords") + coords = [] + for point in raw_coords: + if len(point) < 2: + raise ValueError(f"Breakline {name!r} has an invalid coordinate") + x_coord, y_coord = float(point[0]), float(point[1]) + if not math.isfinite(x_coord) or not math.isfinite(y_coord): + raise ValueError(f"Breakline {name!r} has a non-finite coordinate") + coords.append((x_coord, y_coord)) + if len(coords) < 2: + raise ValueError(f"Breakline {name!r} must contain at least two points") + + cell_sizes = {} + for field_name in ("cell_size_near", "cell_size_far"): + raw_value = raw.get(field_name) + if raw_value is None: + cell_sizes[field_name] = None + continue + value = float(raw_value) + if not math.isfinite(value) or value <= 0: + raise ValueError( + f"Breakline {name!r} {field_name} must be finite and positive" + ) + cell_sizes[field_name] = value + + def bounded_integer(field_name: str, maximum: int) -> int: + raw_value = raw.get(field_name, 0) + value = int(raw_value) + if isinstance(raw_value, float) and not raw_value.is_integer(): + raise ValueError( + f"Breakline {name!r} {field_name} must be an integer" + ) + if value < 0 or value > maximum: + raise ValueError( + f"Breakline {name!r} {field_name} must be between 0 and {maximum}" + ) + return value + + near_repeats = bounded_integer("near_repeats", 255) + protection_radius = bounded_integer("protection_radius", 1) + + normalized.append( + { + "name": name, + "coords": coords, + **cell_sizes, + "near_repeats": near_repeats, + "protection_radius": protection_radius, + } + ) + return normalized + + @staticmethod + @log_call + def replace_breaklines( + geom_file: Union[str, Path], + flow_area_name: str, + breaklines: Sequence[dict], + *, + expected_existing_names: Optional[Sequence[str]] = None, + create_backup: bool = True, + ) -> Optional[Path]: + """Atomically replace every breakline block in a geometry text file. + + HEC-RAS stores ``BreakLine`` blocks as a geometry-wide collection; the + text format does not carry a 2D-area identifier on each block. + ``flow_area_name`` therefore selects the insertion anchor, while this + method intentionally replaces the complete collection. Callers working + with multiple 2D areas must supply the complete retained collection. + + ``expected_existing_names`` is an optimistic-concurrency guard. When + supplied, the existing names and order must match before any write. + Passing an empty ``breaklines`` sequence removes the collection. + """ + geom_file = Path(geom_file) + if not geom_file.exists(): + raise FileNotFoundError(f"Geometry file not found: {geom_file}") + GeomStorage._validate_flow_area_name(flow_area_name) + normalized = GeomStorage._normalize_breakline_specs(breaklines) + + with open(geom_file, "r", encoding="utf-8", errors="replace") as stream: + lines = stream.readlines() + existing_block = GeomStorage._find_storage_area_block(lines, flow_area_name) + if existing_block is None: + raise ValueError(f"Flow area not found: {flow_area_name}") + + ranges = GeomStorage._breakline_block_ranges(lines) + existing_names = [ + lines[start].split("=", 1)[1].strip() for start, _end in ranges + ] + if expected_existing_names is not None: + expected_existing = [str(name).strip() for name in expected_existing_names] + if existing_names != expected_existing: + raise ValueError( + "Existing breakline collection changed before replacement: " + f"expected {expected_existing!r}, found {existing_names!r}" + ) + + new_blocks: List[str] = [] + for breakline in normalized: + new_blocks.extend( + GeomStorage._format_breakline_block( + name=breakline["name"], + coords=breakline["coords"], + cell_size_near=breakline["cell_size_near"], + cell_size_far=breakline["cell_size_far"], + near_repeats=breakline["near_repeats"], + protection_radius=breakline["protection_radius"], + ) + ) + + if ranges: + range_by_start = {start: end for start, end in ranges} + replacement_lines: List[str] = [] + index = 0 + inserted = False + while index < len(lines): + end = range_by_start.get(index) + if end is not None: + if not inserted: + replacement_lines.extend(new_blocks) + inserted = True + index = end + continue + replacement_lines.append(lines[index]) + index += 1 + else: + _start, end_idx, _block = existing_block + insert_idx = end_idx + for idx in range(end_idx, len(lines)): + if lines[idx].startswith( + ( + "BC Line Name=", + "Connection=", + "LCMann Time=", + "Storage Area=", + "River Reach=", + ) + ): + insert_idx = idx + break + replacement_lines = lines[:insert_idx] + new_blocks + lines[insert_idx:] + + written_ranges = GeomStorage._breakline_block_ranges(replacement_lines) + written_names = [ + replacement_lines[start].split("=", 1)[1].strip() + for start, _end in written_ranges + ] + expected_names = [breakline["name"] for breakline in normalized] + if written_names != expected_names: + raise RuntimeError( + f"Breakline replacement validation failed: {written_names!r}" + ) + + backup_path = GeomParser.safe_write_geometry( + geom_file, + replacement_lines, + create_backup=create_backup, + ) + logger.info( + "Replaced %d breaklines with %d in %s", + len(existing_names), + len(expected_names), + geom_file.name, + ) + return backup_path + @staticmethod @log_call def set_breaklines( diff --git a/ras_commander/geom/__init__.py b/ras_commander/geom/__init__.py index e956c7214..94602d492 100644 --- a/ras_commander/geom/__init__.py +++ b/ras_commander/geom/__init__.py @@ -55,7 +55,6 @@ GeomCrossSection, ) from .ManningsFromLandCover import ManningsFromLandCover -from .GeomCrossSection import GeomCrossSection from .GeomStorage import GeomStorage from .GeomProjection import GeomProjection from .GeomLateral import GeomLateral @@ -70,7 +69,13 @@ from .GeomBcLines import GeomBcLines from .GeomMesh import GeomMesh from .GeomPipeNetwork import GeomPipeNetwork -from .GeomMeshDataclasses import MeshResult, BCConflict, BCFixResult +from .GeomMeshDataclasses import ( + BCConflict, + BCFixResult, + DomainContainmentResult, + DomainContainmentViolation, + MeshResult, +) __all__ = [ 'GeomParser', @@ -100,4 +105,6 @@ 'MeshResult', 'BCConflict', 'BCFixResult', + 'DomainContainmentResult', + 'DomainContainmentViolation', ] diff --git a/ras_commander/gui/workflows/mesh_regeneration.py b/ras_commander/gui/workflows/mesh_regeneration.py index f6c833de9..a84f87984 100644 --- a/ras_commander/gui/workflows/mesh_regeneration.py +++ b/ras_commander/gui/workflows/mesh_regeneration.py @@ -8,10 +8,12 @@ This is the critical workflow for Glenn's ras-agent pipeline. Workflow (single attempt): -1. Launch HEC-RAS with project -2. Open RASMapper → triggers HDF regeneration from modified .g## file -3. Save geometry (Ctrl+S) → finalizes the HDF -4. Close HEC-RAS +1. Capture the selected geometry's terrain/classification associations +2. Launch the project with its configured HEC-RAS version +3. Transactionally displace only the selected geometry HDF +4. Open RASMapper and save the exact registered geometry from modified .g## text +5. Validate the imported 2D perimeter, close the owned process tree, and restore + the captured associations before committing the replacement HDF Iterative workflow (with error correction): 1. Attempt mesh regeneration @@ -21,10 +23,12 @@ 5. Retry (up to max_iterations) """ +import os import re import time from pathlib import Path -from typing import Optional, List, Tuple +from typing import Any, Dict, List, Optional, Tuple, Union +from uuid import uuid4 # Win32 imports - Windows only try: @@ -169,7 +173,14 @@ def _check_mesh_valid(geom_hdf_path: Path, mesh_name: Optional[str] = None) -> d """ import h5py - result = {"valid": False, "n_cells": 0, "n_faces": 0, "error": ""} + result = { + "valid": False, + "n_cells": 0, + "n_cell_rows": 0, + "n_virtual_cells": 0, + "n_faces": 0, + "error": "", + } if not geom_hdf_path.exists(): result["error"] = f"HDF file does not exist: {geom_hdf_path}" @@ -214,15 +225,38 @@ def _check_mesh_valid(geom_hdf_path: Path, mesh_name: Optional[str] = None) -> d result["error"] = f"Empty dataset: {ds_name}" return result - n_cells = f[f"{mesh_base}/Cells Center Coordinate"].shape[0] + n_cell_rows = f[f"{mesh_base}/Cells Center Coordinate"].shape[0] n_faces = f[f"{mesh_base}/Faces FacePoint Indexes"].shape[0] + n_cells = n_cell_rows + attributes = f[base].get("Attributes") + if attributes is not None: + data = attributes[()] + fields = data.dtype.names or () + if "Name" in fields and "Cell Count" in fields: + for row in data: + raw_name = row["Name"] + name = ( + raw_name.decode("utf-8", errors="replace") + if isinstance(raw_name, bytes) + else str(raw_name) + ).rstrip("\x00").strip() + if name == mesh_name: + n_cells = int(row["Cell Count"]) + break result["n_cells"] = n_cells + result["n_cell_rows"] = n_cell_rows + result["n_virtual_cells"] = max(0, n_cell_rows - n_cells) result["n_faces"] = n_faces if n_cells < 3: result["error"] = f"Too few cells: {n_cells}" return result + if n_cells > n_cell_rows: + result["error"] = ( + f"Active cell count {n_cells} exceeds coordinate rows {n_cell_rows}" + ) + return result # Check face-cell connectivity (face perimeter connection test) if f"{mesh_base}/Faces Cell Indexes" in f: @@ -242,6 +276,603 @@ def _check_mesh_valid(geom_hdf_path: Path, mesh_name: Optional[str] = None) -> d return result +def _current_plan_geometry_number(ras_obj) -> str: + """Resolve the geometry referenced by the sole ``Current Plan`` record.""" + from ...RasUtils import RasUtils + + matches = re.findall( + r"^Current Plan=p(\d{2})\s*$", + ras_obj.prj_file.read_text(encoding="utf-8", errors="replace"), + flags=re.IGNORECASE | re.MULTILINE, + ) + if len(matches) != 1: + raise ValueError( + f"Expected one Current Plan in {ras_obj.prj_file.name}; found {len(matches)}" + ) + plan_number = RasUtils.normalize_ras_number(matches[0]) + plan_df = getattr(ras_obj, "plan_df", None) + if plan_df is None or "plan_number" not in plan_df: + raise ValueError("Initialized project has no plan_df plan_number column") + rows = plan_df[plan_df["plan_number"].astype(str).str.zfill(2) == plan_number] + if len(rows) != 1: + raise ValueError( + f"Current plan p{plan_number} resolved to {len(rows)} plan_df rows" + ) + row = rows.iloc[0] + value = row.get("geometry_number") + if value is None or str(value).strip().casefold() in {"", "nan", "none", ""}: + value = row.get("Geom File") + if value is None: + raise ValueError(f"Current plan p{plan_number} has no geometry reference") + return RasUtils.normalize_ras_number(str(value).lstrip("gG")) + + +def _resolve_geometry_target( + ras_obj, + *, + geom_number: Optional[Union[str, int]] = None, + geometry_name: Optional[str] = None, +) -> Dict[str, Any]: + """Resolve one registered geometry by number and RASMapper tree identity.""" + from ...RasMap import RasMap + from ...RasUtils import RasUtils + + ras_obj.check_initialized() + mapper_geometries = list(RasMap.list_geometries(ras_obj)) + if geom_number is None and geometry_name is not None: + named = [ + geometry + for geometry in mapper_geometries + if str(geometry.get("name", "")).strip().casefold() + == geometry_name.strip().casefold() + ] + if len(named) != 1: + raise ValueError( + f"Geometry name {geometry_name!r} resolved to {len(named)} " + "RASMapper geometry layers; supply geom_number" + ) + geom_number = named[0].get("geom_number") + if geom_number is None: + geom_number = _current_plan_geometry_number(ras_obj) + geom_num = RasUtils.normalize_ras_number(geom_number) + + geom_df = getattr(ras_obj, "geom_df", None) + if geom_df is None or "geom_number" not in geom_df: + raise ValueError("Initialized project has no geom_df geom_number column") + rows = geom_df[geom_df["geom_number"].astype(str).str.zfill(2) == geom_num] + if len(rows) != 1: + raise ValueError(f"Geometry g{geom_num} resolved to {len(rows)} geom_df rows") + geom_file = Path(ras_obj.project_folder) / f"{ras_obj.project_name}.g{geom_num}" + if not geom_file.is_file(): + raise FileNotFoundError(f"Geometry text file not found: {geom_file}") + + mapper_matches = [ + geometry + for geometry in mapper_geometries + if str(geometry.get("geom_number", "")).zfill(2) == geom_num + ] + if len(mapper_matches) != 1: + raise ValueError( + f"Geometry g{geom_num} resolved to {len(mapper_matches)} RASMapper layers" + ) + tree_name = str(mapper_matches[0].get("name", "")).strip() + if not tree_name: + raise ValueError(f"RASMapper geometry g{geom_num} has no display name") + duplicate_tree_names = [ + geometry + for geometry in mapper_geometries + if str(geometry.get("name", "")).strip().casefold() == tree_name.casefold() + ] + if len(duplicate_tree_names) != 1: + raise ValueError( + f"RASMapper tree name {tree_name!r} is not unique; cannot target g{geom_num} safely" + ) + if geometry_name is not None and tree_name.casefold() != geometry_name.strip().casefold(): + raise ValueError( + f"Geometry g{geom_num} is named {tree_name!r}, not {geometry_name!r}" + ) + + return { + "geom_number": geom_num, + "geometry_name": tree_name, + "geom_file": geom_file, + "geom_hdf": Path(f"{geom_file}.hdf"), + } + + +def _select_text_flow_area(geom_file: Path, flow_area_name: Optional[str]): + """Return one exact 2D flow-area polygon from geometry text.""" + from ...geom.GeomStorage import GeomStorage + + areas = GeomStorage.get_storage_area_polygons(geom_file, exclude_2d=False) + if areas is None or areas.empty: + raise ValueError(f"No storage-area polygons found in {geom_file.name}") + if "is_2d" in areas: + areas = areas[areas["is_2d"]] + if flow_area_name is None: + if len(areas) != 1: + raise ValueError( + f"Geometry {geom_file.name} has {len(areas)} 2D flow areas; " + "supply flow_area_name" + ) + row = areas.iloc[0] + else: + name_column = "name" if "name" in areas else "Name" + matches = areas[areas[name_column].astype(str) == flow_area_name] + if len(matches) != 1: + raise ValueError( + f"2D flow area {flow_area_name!r} resolved to {len(matches)} text polygons" + ) + row = matches.iloc[0] + name_column = "name" if "name" in areas else "Name" + polygon = row.geometry + if polygon is None or polygon.is_empty or not polygon.is_valid: + raise ValueError(f"2D flow area {row[name_column]!r} has an invalid text perimeter") + return str(row[name_column]), polygon + + +def _geometry_hdf_stats(ras_obj) -> Dict[str, tuple[int, int]]: + """Return size/mtime evidence for every registered geometry HDF.""" + stats: Dict[str, tuple[int, int]] = {} + geom_df = getattr(ras_obj, "geom_df", None) + if geom_df is None or "geom_number" not in geom_df: + return stats + for value in geom_df["geom_number"]: + number = str(value).zfill(2) + path = Path(ras_obj.project_folder) / f"{ras_obj.project_name}.g{number}.hdf" + if path.is_file(): + stat = path.stat() + stats[str(path.resolve())] = (stat.st_size, stat.st_mtime_ns) + return stats + + +def _read_hdf_flow_area_polygon(geom_hdf: Path, flow_area_name: str): + """Read a registered 2D-area polygon before or after mesh compilation. + + A newly rebuilt RAS Mapper geometry HDF contains the collection-level + ``Polygon *`` datasets before it contains the per-area ``Perimeter`` and + mesh datasets. Reading the collection therefore validates the important + text-to-HDF import boundary without requiring a mesh to exist yet. + """ + import h5py + import numpy as np + from shapely.geometry import Polygon + from shapely.ops import unary_union + + group_path = "Geometry/2D Flow Areas" + with h5py.File(geom_hdf, "r") as hdf: + if group_path not in hdf: + raise ValueError("Geometry HDF has no 2D flow-area collection") + group = hdf[group_path] + required = {"Attributes", "Polygon Info", "Polygon Parts", "Polygon Points"} + missing = sorted(required.difference(group.keys())) + if missing: + raise ValueError( + "Geometry HDF 2D flow-area collection is missing: " + + ", ".join(missing) + ) + + attributes = group["Attributes"][()] + if attributes.dtype.names is None or "Name" not in attributes.dtype.names: + raise ValueError("Geometry HDF 2D flow-area Attributes has no Name field") + + def decode(value) -> str: + if isinstance(value, (bytes, np.bytes_)): + return bytes(value).decode("utf-8", errors="replace").rstrip("\x00").strip() + return str(value).rstrip("\x00").strip() + + names = [decode(value) for value in attributes["Name"]] + matches = [index for index, name in enumerate(names) if name == flow_area_name] + if len(matches) != 1: + raise ValueError( + f"2D flow area {flow_area_name!r} resolved to {len(matches)} " + "collection polygons" + ) + row_index = matches[0] + + polygon_info = np.asarray(group["Polygon Info"][()]) + polygon_parts = np.asarray(group["Polygon Parts"][()]) + polygon_points = np.asarray(group["Polygon Points"][()], dtype=float) + if polygon_info.ndim != 2 or polygon_info.shape[1] < 4: + raise ValueError("Geometry HDF Polygon Info must be an Nx4 dataset") + if row_index >= polygon_info.shape[0]: + raise ValueError("Geometry HDF polygon row count does not match Attributes") + + point_start, point_count, part_start, part_count = ( + int(value) for value in polygon_info[row_index, :4] + ) + if point_count < 3 or part_count < 1: + raise ValueError(f"2D flow area {flow_area_name!r} has no polygon ring") + area_points = polygon_points[point_start : point_start + point_count] + if area_points.shape != (point_count, 2): + raise ValueError(f"2D flow area {flow_area_name!r} has truncated points") + + rings = [] + for raw_start, raw_count in polygon_parts[part_start : part_start + part_count, :2]: + ring_start, ring_count = int(raw_start), int(raw_count) + # HEC-RAS files in the field use absolute point indexes. Accept a + # relative index only when the absolute range is outside this row. + if not (point_start <= ring_start < point_start + point_count): + ring_start += point_start + ring = polygon_points[ring_start : ring_start + ring_count] + if ring.shape[0] >= 3: + rings.append(ring) + if not rings: + rings = [area_points] + + shell = Polygon(rings[0]) + if len(rings) == 1: + polygon = shell + else: + holes = [ring for ring in rings[1:] if shell.covers(Polygon(ring))] + islands = [Polygon(ring) for ring in rings[1:] if not shell.covers(Polygon(ring))] + polygon = Polygon(rings[0], holes) + if islands: + polygon = unary_union([polygon, *islands]) + if polygon.is_empty or not polygon.is_valid: + raise ValueError(f"2D flow area {flow_area_name!r} has an invalid polygon") + return polygon + + +def _perimeter_validation( + geom_hdf: Path, + flow_area_name: str, + expected_polygon, + *, + coordinate_tolerance: Optional[float] = None, +) -> dict: + """Compare an exact text perimeter with its compiled HDF polygon.""" + result = { + "valid": False, + "flow_area_name": flow_area_name, + "coordinate_tolerance": coordinate_tolerance, + "hausdorff_distance": None, + "symmetric_difference_area": None, + "text_area": float(expected_polygon.area), + "hdf_area": None, + "error": "", + } + if not geom_hdf.is_file(): + result["error"] = f"Geometry HDF does not exist: {geom_hdf.name}" + return result + try: + actual_polygon = _read_hdf_flow_area_polygon(geom_hdf, flow_area_name) + except (OSError, KeyError, TypeError, ValueError) as exc: + result["error"] = str(exc) + return result + + min_x, min_y, max_x, max_y = expected_polygon.bounds + scale = max(max_x - min_x, max_y - min_y, 1.0) + tolerance = ( + float(coordinate_tolerance) + if coordinate_tolerance is not None + else max(1e-6, scale * 1e-9) + ) + if tolerance < 0: + raise ValueError("coordinate_tolerance must be non-negative") + hausdorff = float(expected_polygon.boundary.hausdorff_distance(actual_polygon.boundary)) + symmetric_area = float(expected_polygon.symmetric_difference(actual_polygon).area) + area_tolerance = max( + tolerance * tolerance, + float(expected_polygon.length) * tolerance, + float(expected_polygon.area) * 1e-10, + ) + result.update( + { + "coordinate_tolerance": tolerance, + "hausdorff_distance": hausdorff, + "symmetric_difference_area": symmetric_area, + "area_tolerance": area_tolerance, + "hdf_area": float(actual_polygon.area), + } + ) + if hausdorff > tolerance or symmetric_area > area_tolerance: + result["error"] = ( + "Compiled HDF perimeter does not match geometry text: " + f"Hausdorff={hausdorff:g} (limit {tolerance:g}), " + f"symmetric area={symmetric_area:g} (limit {area_tolerance:g})" + ) + return result + result["valid"] = True + return result + + +def _prepare_geometry_refresh_context( + ras_obj, + *, + geom_number: Optional[Union[str, int]], + geometry_name: Optional[str], + flow_area_name: Optional[str], + coordinate_tolerance: Optional[float], +) -> Dict[str, Any]: + """Resolve all identities and semantic expectations before GUI work.""" + target = _resolve_geometry_target( + ras_obj, + geom_number=geom_number, + geometry_name=geometry_name, + ) + area_name, expected_polygon = _select_text_flow_area( + target["geom_file"], flow_area_name + ) + # ``Edit Geometry`` is a context-menu action on the registered geometry + # root, not on its ``2D Flow Areas`` child. The flow area remains an + # independent semantic selector for post-save validation. + target_path = [target["geometry_name"]] + pre_stats = _geometry_hdf_stats(ras_obj) + pre_association_paths = _capture_geometry_association_paths(target["geom_hdf"]) + pre_perimeter = _perimeter_validation( + target["geom_hdf"], + area_name, + expected_polygon, + coordinate_tolerance=coordinate_tolerance, + ) + return { + **target, + "flow_area_name": area_name, + "expected_polygon": expected_polygon, + "target_path": target_path, + "coordinate_tolerance": coordinate_tolerance, + "pre_hdf_stats": pre_stats, + "pre_perimeter_validation": pre_perimeter, + "pre_geometry_association_paths": pre_association_paths, + } + + +def _capture_geometry_association_paths(geom_hdf: Path) -> dict[str, Path]: + """Capture all existing, resolvable layer associations before HDF import.""" + target = Path(geom_hdf) + if not target.is_file(): + return {} + from ...geom.GeomMesh import GeomMesh + + association = GeomMesh.get_geometry_association(target) + paths: dict[str, Path] = {} + for key in ( + "terrain_hdf_path", + "landcover_hdf_path", + "infiltration_hdf_path", + "sediment_soils_hdf_path", + ): + value = association.get(key) + if not value: + continue + resolved = Path(value) + if not resolved.is_file(): + raise FileNotFoundError( + f"Cannot preserve {key}; associated artifact is missing: {resolved}" + ) + paths[key] = resolved + return paths + + +def _restore_geometry_association(context: dict) -> dict: + """Restore captured layer associations on the rebuilt exact geometry HDF.""" + expected = dict(context.get("pre_geometry_association_paths") or {}) + evidence = { + "restored": False, + "expected_paths": {key: str(value) for key, value in expected.items()}, + "observed": {}, + } + if not expected: + return evidence + + from ...geom.GeomMesh import GeomMesh + + target = Path(context["geom_hdf"]) + GeomMesh.set_geometry_association( + target, + ras_object=context.get("ras_object"), + validate=True, + **expected, + ) + observed = GeomMesh.get_geometry_association(target) + evidence["restored"] = True + evidence["observed"] = { + key: observed.get(key) + for key in expected + } + return evidence + + +def _validate_geometry_import(context: dict) -> dict: + """Validate the exact text-to-HDF import without requiring a mesh yet.""" + target_hdf = Path(context["geom_hdf"]) + perimeter = _perimeter_validation( + target_hdf, + context["flow_area_name"], + context["expected_polygon"], + coordinate_tolerance=context.get("coordinate_tolerance"), + ) + if not perimeter["valid"]: + raise RuntimeError(perimeter["error"]) + + ras_obj = context["ras_object"] + before_stats = context["pre_hdf_stats"] + after_stats = _geometry_hdf_stats(ras_obj) + target_key = str(target_hdf.resolve()) + other_changes = [] + for path in sorted(set(before_stats) | set(after_stats)): + if path == target_key: + continue + if before_stats.get(path) != after_stats.get(path): + other_changes.append(path) + if other_changes: + raise RuntimeError( + "RASMapper changed non-target geometry HDFs: " + + ", ".join(Path(path).name for path in other_changes) + ) + + target_before = before_stats.get(target_key) + target_after = after_stats.get(target_key) + pre_perimeter = context["pre_perimeter_validation"] + if not pre_perimeter["valid"] and target_before == target_after: + raise RuntimeError( + f"Target geometry {target_hdf.name} matched after save but its file evidence " + "did not change from the known-stale precondition" + ) + + return { + "geom_number": context["geom_number"], + "geometry_name": context["geometry_name"], + "geom_file": str(context["geom_file"]), + "geom_hdf": str(target_hdf), + "flow_area_name": context["flow_area_name"], + "pre_perimeter": pre_perimeter, + "post_perimeter": perimeter, + "other_geometry_hdfs_unchanged": True, + } + + +def _validate_geometry_refresh(context: dict) -> dict: + """Fail closed unless the exact geometry import and mesh are current.""" + result = _validate_geometry_import(context) + target_hdf = Path(context["geom_hdf"]) + mesh = _check_mesh_valid(target_hdf, mesh_name=context["flow_area_name"]) + if not mesh["valid"]: + raise RuntimeError(f"Compiled target mesh is invalid: {mesh['error']}") + result["mesh"] = mesh + return result + + +def _begin_geometry_hdf_transaction(context: dict) -> None: + """Temporarily remove the exact HDF so HEC-RAS imports geometry text.""" + target = Path(context["geom_hdf"]) + temporary_backup = target.with_name( + f".{target.name}.rascommander-{uuid4().hex}.bak" + ) + had_original = target.is_file() + if had_original: + os.replace(target, temporary_backup) + context["hdf_transaction"] = { + "target": target, + "temporary_backup": temporary_backup, + "had_original": had_original, + } + + +def _finish_geometry_hdf_transaction( + context: dict, + *, + success: bool, + keep_backup: bool, +) -> dict: + """Commit a validated HDF or restore the exact target on failure.""" + transaction = context["hdf_transaction"] + target = Path(transaction["target"]) + temporary_backup = Path(transaction["temporary_backup"]) + had_original = bool(transaction["had_original"]) + evidence = { + "target": str(target), + "had_original": had_original, + "rolled_back": False, + "backup": None, + } + + if not success: + if target.exists(): + target.unlink() + if had_original and temporary_backup.exists(): + os.replace(temporary_backup, target) + evidence["rolled_back"] = True + return evidence + + if not target.is_file(): + if had_original and temporary_backup.exists(): + os.replace(temporary_backup, target) + evidence["rolled_back"] = True + raise RuntimeError(f"HEC-RAS did not rebuild geometry HDF: {target}") + + if had_original and temporary_backup.exists(): + if keep_backup: + committed_backup = target.with_name(f"{target.name}.pre-rasmapper.bak") + counter = 1 + while committed_backup.exists(): + committed_backup = target.with_name( + f"{target.name}.pre-rasmapper.bak{counter}" + ) + counter += 1 + os.replace(temporary_backup, committed_backup) + evidence["backup"] = str(committed_backup) + else: + temporary_backup.unlink() + return evidence + + +def _capture_owned_process_tree(process) -> list: + """Capture only the process tree rooted at the workflow's Ras.exe.""" + if process is None: + return [] + try: + import psutil + + root = psutil.Process(int(process.pid)) + return [root, *root.children(recursive=True)] + except Exception as exc: # noqa: BLE001 + logger.debug("Could not capture owned HEC-RAS process tree: %s", exc) + return [] + + +def _supervise_owned_process_exit(process, owned_processes: list) -> dict: + """Gracefully wait, then terminate/kill only captured owned processes.""" + evidence = { + "root_pid": int(process.pid) if process is not None else None, + "observed_pids": sorted({int(proc.pid) for proc in owned_processes}), + "terminated_pids": [], + "killed_pids": [], + "survivor_pids": [], + } + if process is None: + return evidence + + try: + process.wait(timeout=10) + except Exception: # noqa: BLE001 + pass + + try: + import psutil + except ImportError: + if process.poll() is None: + process.terminate() + try: + process.wait(timeout=3) + except Exception: # noqa: BLE001 + process.kill() + process.wait(timeout=3) + return evidence + + if not owned_processes: + try: + owned_processes = [psutil.Process(int(process.pid))] + evidence["observed_pids"] = [int(process.pid)] + except (psutil.NoSuchProcess, psutil.AccessDenied): + owned_processes = [] + + alive = [] + for owned in owned_processes: + try: + if owned.is_running() and owned.status() != psutil.STATUS_ZOMBIE: + alive.append(owned) + except (psutil.NoSuchProcess, psutil.AccessDenied): + continue + for owned in reversed(alive): + try: + evidence["terminated_pids"].append(int(owned.pid)) + owned.terminate() + except (psutil.NoSuchProcess, psutil.AccessDenied): + continue + _gone, alive = psutil.wait_procs(alive, timeout=3) + for owned in alive: + try: + evidence["killed_pids"].append(int(owned.pid)) + owned.kill() + except (psutil.NoSuchProcess, psutil.AccessDenied): + continue + _gone, survivors = psutil.wait_procs(alive, timeout=3) + evidence["survivor_pids"] = sorted(int(owned.pid) for owned in survivors) + return evidence + + # --------------------------------------------------------------------------- # Perimeter simplification # --------------------------------------------------------------------------- @@ -308,6 +939,115 @@ class MeshRegenerationWorkflow: All methods are static and decorated with @log_call. """ + # ------------------------------------------------------------------ + # Exact text-to-HDF import + # ------------------------------------------------------------------ + + @staticmethod + @log_call + def refresh_geometry_hdf_from_text( + geometry_name: Optional[str] = None, + flow_area_name: Optional[str] = None, + ras_object=None, + timeout: int = 600, + *, + geom_number: Optional[Union[str, int]] = None, + coordinate_tolerance: Optional[float] = None, + keep_backup: bool = True, + ) -> WorkflowResult: + """Rebuild one exact geometry HDF from its current ``.g##`` text. + + HEC-RAS/RAS Mapper treats an existing geometry HDF as authoritative + during editing, so merely opening and saving does not import external + text changes. This operation transactionally displaces only the exact + selected HDF, opens the task-local project through its configured + HEC-RAS version, saves the exact registered geometry, and validates the + collection-level 2D perimeter before committing the replacement. + + The geometry defaults to the sole current plan's geometry. Supplying + both ``geom_number`` and ``geometry_name`` makes them an identity + cross-check. Other registered geometry HDF files must remain unchanged. + A failed import restores the original target HDF. Existing terrain, + land-cover, infiltration, and sediment associations are captured before + import and restored on the rebuilt HDF; a missing associated artifact or + failed association validation aborts and rolls back the transaction. + + This method compiles geometry features but does not create computation + cells. Call :meth:`GeomMesh.generate` after applying HDF-only feature + edits such as refinement-region replacement. + """ + from ...RasPrj import ras + + ras_obj = ras_object or ras + ras_obj.check_initialized() + context = _prepare_geometry_refresh_context( + ras_obj, + geom_number=geom_number, + geometry_name=geometry_name, + flow_area_name=flow_area_name, + coordinate_tolerance=coordinate_tolerance, + ) + context.update( + { + "ras_object": ras_obj, + "timeout": timeout, + "close_after": True, + "force_text_import": True, + } + ) + _begin_geometry_hdf_transaction(context) + + result = WorkflowExecutor.execute( + MeshRegenerationWorkflow._build_single_attempt_steps( + context, + require_mesh=False, + ), + context, + workflow_name=f"GeometryTextImport[g{context['geom_number']}]", + ) + cleanup_error = None + if not context.get("closed"): + try: + MeshRegenerationWorkflow._step_close(context) + except Exception as exc: # noqa: BLE001 + logger.warning("Geometry text import cleanup failed: %s", exc) + cleanup_error = exc + if cleanup_error is not None: + result.success = False + result.error = cleanup_error + result.steps_failed.append("Close owned HEC-RAS process tree") + + if result.success: + try: + association = _restore_geometry_association(context) + except Exception as exc: # noqa: BLE001 + result.success = False + result.error = exc + result.steps_failed.append("Restore geometry associations") + association = { + "restored": False, + "error": str(exc), + } + result.step_results["Restore geometry associations"] = association + + try: + transaction = _finish_geometry_hdf_transaction( + context, + success=result.success, + keep_backup=keep_backup, + ) + except Exception as exc: # noqa: BLE001 + result.success = False + result.error = exc + result.steps_failed.append("Commit geometry HDF transaction") + transaction = { + "target": str(context["geom_hdf"]), + "rolled_back": True, + "backup": None, + } + result.step_results["Geometry HDF transaction"] = transaction + return result + # ------------------------------------------------------------------ # Single-attempt mesh regeneration # ------------------------------------------------------------------ @@ -320,6 +1060,9 @@ def regenerate_mesh( ras_object=None, timeout: int = 600, close_after: bool = True, + *, + geom_number: Optional[Union[str, int]] = None, + coordinate_tolerance: Optional[float] = None, ) -> WorkflowResult: """ Single mesh regeneration attempt. @@ -328,25 +1071,58 @@ def regenerate_mesh( saves (triggering HDF regeneration), and checks the result. Args: - geometry_name: Name of the geometry. Auto-detected if None. + geometry_name: Exact RASMapper geometry-layer name. When supplied + with ``geom_number``, both identities must agree. flow_area_name: Name of the 2D flow area. Auto-detected if None. ras_object: Optional RasPrj object instance. timeout: Max seconds to wait for mesh generation. Default 600. close_after: If True, close RASMapper and HEC-RAS when done. + geom_number: Exact geometry number, such as ``"03"``. When both + geometry selectors are omitted, resolve the geometry referenced + by the sole current plan; never fall back to the first ``Geom + File=`` registration. + coordinate_tolerance: Optional project-unit limit for comparing the + text perimeter with the compiled HDF perimeter. Returns: - WorkflowResult. step_results['mesh_check'] contains validation dict. + WorkflowResult. ``step_results['Validate exact geometry HDF']`` + contains target identity, perimeter, mesh, and non-target-file + validation evidence. """ - context = { - 'geometry_name': geometry_name, - 'flow_area_name': flow_area_name, - 'ras_object': ras_object, - 'timeout': timeout, - 'close_after': close_after, - } + from ...RasPrj import ras + + ras_obj = ras_object or ras + ras_obj.check_initialized() + context = _prepare_geometry_refresh_context( + ras_obj, + geom_number=geom_number, + geometry_name=geometry_name, + flow_area_name=flow_area_name, + coordinate_tolerance=coordinate_tolerance, + ) + context.update( + { + "ras_object": ras_obj, + "timeout": timeout, + "close_after": close_after, + } + ) steps = MeshRegenerationWorkflow._build_single_attempt_steps(context) - return WorkflowExecutor.execute(steps, context, workflow_name="MeshRegeneration") + result = WorkflowExecutor.execute( + steps, + context, + workflow_name=f"MeshRegeneration[g{context['geom_number']}]", + ) + if close_after and not context.get("closed"): + try: + MeshRegenerationWorkflow._step_close(context) + except Exception as exc: # noqa: BLE001 + logger.warning("Mesh regeneration cleanup failed: %s", exc) + result.success = False + result.error = exc + result.steps_failed.append("Close owned HEC-RAS process tree") + return result # ------------------------------------------------------------------ # Iterative mesh regeneration with error correction @@ -361,6 +1137,11 @@ def regenerate_mesh_iterative( initial_cell_size: Optional[float] = None, cell_size_increase_factor: float = 1.3, simplify_tolerance_factor: float = 0.25, + *, + geom_number: Optional[Union[str, int]] = None, + geometry_name: Optional[str] = None, + flow_area_name: Optional[str] = None, + coordinate_tolerance: Optional[float] = None, ) -> WorkflowResult: """ Iterative mesh regeneration with face perimeter error correction. @@ -382,6 +1163,11 @@ def regenerate_mesh_iterative( initial_cell_size: Starting cell size. Read from .g## if None. cell_size_increase_factor: Multiply cell size by this on each retry. Default 1.3. simplify_tolerance_factor: simplify_tolerance = cell_size * this. Default 0.25. + geom_number: Exact geometry number. Defaults to the current plan's + geometry, never the first project geometry registration. + geometry_name: Optional exact RASMapper name cross-check. + flow_area_name: Exact 2D flow-area name when geometry has more than one. + coordinate_tolerance: Optional text/HDF perimeter comparison limit. Returns: WorkflowResult with: @@ -395,26 +1181,21 @@ def regenerate_mesh_iterative( ras_obj = ras_object or ras ras_obj.check_initialized() - # Locate geometry file and HDF - geom_file = _find_geometry_file(ras_obj) - if geom_file is None: - return WorkflowResult( - success=False, - error=RuntimeError("No geometry file found in project"), - ) - - geom_hdf = geom_file.parent / (geom_file.name + ".hdf") - area_name = _get_2d_area_name(geom_file) - - if area_name is None: - return WorkflowResult( - success=False, - error=RuntimeError(f"No 2D flow area found in {geom_file.name}"), - ) + target = _resolve_geometry_target( + ras_obj, + geom_number=geom_number, + geometry_name=geometry_name, + ) + geom_file = target["geom_file"] + geom_hdf = target["geom_hdf"] + area_name, area_polygon = _select_text_flow_area(geom_file, flow_area_name) # Read initial state cell_size = initial_cell_size or _read_cell_size(geom_file) or 500.0 - original_coords = _read_perimeter_coords(geom_file) + original_coords = [ + (float(x_coord), float(y_coord)) + for x_coord, y_coord in area_polygon.exterior.coords + ] current_coords = list(original_coords) if original_coords else None if current_coords is None: @@ -436,9 +1217,13 @@ def regenerate_mesh_iterative( # Step 1: Open RASMapper, save, close attempt = MeshRegenerationWorkflow.regenerate_mesh( + geometry_name=target["geometry_name"], + flow_area_name=area_name, ras_object=ras_obj, timeout=timeout, close_after=True, + geom_number=target["geom_number"], + coordinate_tolerance=coordinate_tolerance, ) if not attempt.success: @@ -506,9 +1291,18 @@ def regenerate_mesh_iterative( # ------------------------------------------------------------------ @staticmethod - def _build_single_attempt_steps(context: dict) -> list: + def _build_single_attempt_steps( + context: dict, + *, + require_mesh: bool = True, + ) -> list: """Build step sequence for a single regeneration attempt.""" steps = [ + WorkflowStep( + name="Verify HEC-RAS is closed", + action=MeshRegenerationWorkflow._step_verify_hecras_closed, + max_retries=1, + ), WorkflowStep( name="Launch HEC-RAS", action=MeshRegenerationWorkflow._step_launch_hecras, @@ -527,6 +1321,12 @@ def _build_single_attempt_steps(context: dict) -> list: max_retries=1, timeout=context.get('timeout', 600), ), + WorkflowStep( + name="Select exact geometry for editing", + action=MeshRegenerationWorkflow._step_select_geometry, + max_retries=2, + retry_delay=2.0, + ), WorkflowStep( name="Save geometry (trigger HDF regeneration)", action=MeshRegenerationWorkflow._step_save_geometry, @@ -539,6 +1339,19 @@ def _build_single_attempt_steps(context: dict) -> list: max_retries=1, timeout=context.get('timeout', 600), ), + WorkflowStep( + name=( + "Validate exact geometry HDF" + if require_mesh + else "Validate exact geometry import" + ), + action=( + MeshRegenerationWorkflow._step_validate_geometry + if require_mesh + else MeshRegenerationWorkflow._step_validate_geometry_import + ), + max_retries=1, + ), ] if context.get('close_after', True): @@ -547,11 +1360,18 @@ def _build_single_attempt_steps(context: dict) -> list: action=MeshRegenerationWorkflow._step_close, max_retries=2, retry_delay=1.0, - required=False, + required=True, )) return steps + @staticmethod + def _step_verify_hecras_closed(context: dict) -> None: + """Refuse to attach the one-shot operation to an unrelated GUI session.""" + from .xsec_update import RasMapperLayerCommandWorkflow + + RasMapperLayerCommandWorkflow._step_verify_hecras_closed(context) + @staticmethod def _step_launch_hecras(context: dict) -> None: """Launch HEC-RAS and store process/hwnd in context.""" @@ -571,7 +1391,7 @@ def _step_open_rasmapper(context: dict) -> None: try: win32gui.SetForegroundWindow(hwnd) - except: + except Exception: pass time.sleep(0.5) @@ -593,6 +1413,13 @@ def _step_wait_for_rasmapper(context: dict) -> None: context['rasmapper_hwnd'] = result[0] context['rasmapper_title'] = result[1] + @staticmethod + def _step_select_geometry(context: dict) -> None: + """Enter edit mode through the exact geometry's 2D Flow Areas node.""" + from .xsec_update import RasMapperLayerCommandWorkflow + + RasMapperLayerCommandWorkflow._step_start_geometry_editing(context) + @staticmethod def _step_save_geometry(context: dict) -> None: """Save geometry in RASMapper via Ctrl+S, triggering HDF regeneration.""" @@ -600,7 +1427,7 @@ def _step_save_geometry(context: dict) -> None: try: win32gui.SetForegroundWindow(rasmapper_hwnd) - except: + except Exception: pass time.sleep(0.5) @@ -626,12 +1453,23 @@ def _step_wait_for_save(context: dict) -> None: if not RasMapperElements.wait_for_rasmapper_idle(rasmapper_hwnd, timeout=timeout): logger.warning("RASMapper may still be processing, but continuing...") + @staticmethod + def _step_validate_geometry(context: dict) -> dict: + """Validate target identity, perimeter, mesh, and non-target isolation.""" + return _validate_geometry_refresh(context) + + @staticmethod + def _step_validate_geometry_import(context: dict) -> dict: + """Validate exact text import before computation-cell generation.""" + return _validate_geometry_import(context) + @staticmethod def _step_close(context: dict) -> None: """Close RASMapper and HEC-RAS.""" rasmapper_hwnd = context.get('rasmapper_hwnd') hecras_hwnd = context.get('hecras_hwnd') hecras_process = context.get('hecras_process') + owned_processes = _capture_owned_process_tree(hecras_process) if rasmapper_hwnd: Win32Primitives.close_window(rasmapper_hwnd) @@ -642,12 +1480,15 @@ def _step_close(context: dict) -> None: if hecras_hwnd: Win32Primitives.close_window(hecras_hwnd) - if hecras_process: - try: - hecras_process.wait(timeout=10) - except: - pass + cleanup = _supervise_owned_process_exit(hecras_process, owned_processes) + context["owned_process_cleanup"] = cleanup + if cleanup["survivor_pids"]: + raise RuntimeError( + "Owned HEC-RAS processes survived cleanup: " + + ", ".join(str(pid) for pid in cleanup["survivor_pids"]) + ) + context["closed"] = True logger.debug("RASMapper and HEC-RAS closed") @@ -655,26 +1496,19 @@ def _step_close(context: dict) -> None: # Helpers # --------------------------------------------------------------------------- -def _find_geometry_file(ras_obj) -> Optional[Path]: - """Find the current geometry file from the RAS project.""" +def _find_geometry_file( + ras_obj, + geom_number: Optional[Union[str, int]] = None, + geometry_name: Optional[str] = None, +) -> Optional[Path]: + """Resolve an exact or current-plan geometry without first-entry fallback.""" try: - # Read the .prj file to find current geometry - prj_text = ras_obj.prj_file.read_text(errors="replace") - match = re.search(r"Geom File=(\S+)", prj_text) - if match: - geom_ext = match.group(1).strip() - # The geometry file is project_name.g01 etc. - geom_file = ras_obj.project_folder / f"{ras_obj.project_name}.{geom_ext}" - if geom_file.exists(): - return geom_file - - # Fallback: find any .g## file - for g in sorted(ras_obj.project_folder.glob(f"{ras_obj.project_name}.g*")): - if g.suffix and g.suffix[1:].startswith("g") and not g.suffix.endswith(".hdf"): - return g - + return _resolve_geometry_target( + ras_obj, + geom_number=geom_number, + geometry_name=geometry_name, + )["geom_file"] except Exception as e: - logger.warning("Could not find geometry file") + logger.warning("Could not resolve exact geometry file") logger.debug("Geometry file discovery failure: %s", e) - - return None + return None diff --git a/ras_commander/hdf/HdfBndry.py b/ras_commander/hdf/HdfBndry.py index ae7a45c39..64b01277c 100644 --- a/ras_commander/hdf/HdfBndry.py +++ b/ras_commander/hdf/HdfBndry.py @@ -25,7 +25,7 @@ """ from pathlib import Path -from typing import Dict, List, Optional, Union, Any +from typing import Optional import h5py import numpy as np import pandas as pd @@ -33,9 +33,8 @@ from shapely.geometry import LineString, MultiLineString, Polygon, MultiPolygon, Point from .HdfBase import HdfBase from .HdfUtils import HdfUtils -from .HdfMesh import HdfMesh -from ..Decorators import standardize_input, log_call -from ..LoggingConfig import setup_logging, get_logger +from ..Decorators import standardize_input +from ..LoggingConfig import get_logger logger = get_logger(__name__) @@ -156,6 +155,10 @@ def get_breaklines(hdf_path: Path) -> gpd.GeoDataFrame: # Initialize lists to store valid breakline data valid_ids = [] valid_names = [] + valid_spacing_near = [] + valid_spacing_far = [] + valid_near_repeats = [] + valid_protection_radius = [] valid_geoms = [] # Track invalid breaklines for summary @@ -205,6 +208,27 @@ def get_breaklines(hdf_path: Path) -> gpd.GeoDataFrame: valid_ids.append(idx) valid_names.append(name) + fields = attributes.dtype.names or () + valid_spacing_near.append( + float(attributes["Cell Spacing Near"][idx]) + if "Cell Spacing Near" in fields + else None + ) + valid_spacing_far.append( + float(attributes["Cell Spacing Far"][idx]) + if "Cell Spacing Far" in fields + else None + ) + valid_near_repeats.append( + int(attributes["Near Repeats"][idx]) + if "Near Repeats" in fields + else 0 + ) + valid_protection_radius.append( + int(attributes["Protection Radius"][idx]) + if "Protection Radius" in fields + else 0 + ) valid_geoms.append(geom) except Exception as e: @@ -240,6 +264,10 @@ def get_breaklines(hdf_path: Path) -> gpd.GeoDataFrame: { "bl_id": valid_ids, "Name": valid_names, + "cell_spacing_near": valid_spacing_near, + "cell_spacing_far": valid_spacing_far, + "near_repeats": valid_near_repeats, + "protection_radius": valid_protection_radius, "geometry": valid_geoms }, geometry="geometry", diff --git a/ras_commander/hdf/HdfMesh.py b/ras_commander/hdf/HdfMesh.py index 5fb515294..dea276f84 100644 --- a/ras_commander/hdf/HdfMesh.py +++ b/ras_commander/hdf/HdfMesh.py @@ -172,10 +172,44 @@ def get_mesh_areas(hdf_path: Path) -> 'GeoDataFrame': mesh_area_names = HdfMesh.get_mesh_area_names(hdf_path) if not mesh_area_names: return GeoDataFrame() - mesh_area_polygons = [ - Polygon(hdf_file["Geometry/2D Flow Areas/{}/Perimeter".format(n)][()]) - for n in mesh_area_names - ] + group = hdf_file["Geometry/2D Flow Areas"] + mesh_area_polygons = [] + polygon_info = group.get("Polygon Info") + polygon_parts = group.get("Polygon Parts") + polygon_points = group.get("Polygon Points") + for index, name in enumerate(mesh_area_names): + perimeter_path = f"{name}/Perimeter" + if perimeter_path in group: + mesh_area_polygons.append(Polygon(group[perimeter_path][()])) + continue + + if polygon_info is None or polygon_points is None: + raise KeyError( + f"2D flow area {name!r} has neither Perimeter nor " + "collection-level Polygon datasets" + ) + point_start, point_count, part_start, part_count = ( + int(value) for value in polygon_info[index, :4] + ) + points = polygon_points[()] + if part_count <= 1 or polygon_parts is None: + mesh_area_polygons.append( + Polygon(points[point_start : point_start + point_count]) + ) + continue + + parts = polygon_parts[part_start : part_start + part_count, :2] + rings = [] + for raw_start, raw_count in parts: + ring_start, ring_count = int(raw_start), int(raw_count) + if not (point_start <= ring_start < point_start + point_count): + ring_start += point_start + ring = points[ring_start : ring_start + ring_count] + if len(ring) >= 3: + rings.append(ring) + if not rings: + raise ValueError(f"2D flow area {name!r} has no valid polygon ring") + mesh_area_polygons.append(Polygon(rings[0], rings[1:])) return GeoDataFrame( {"mesh_name": mesh_area_names, "geometry": mesh_area_polygons}, geometry="geometry", diff --git a/ras_commander/hdf/HdfResultsPlan.py b/ras_commander/hdf/HdfResultsPlan.py index 92c026f30..9e2c65241 100644 --- a/ras_commander/hdf/HdfResultsPlan.py +++ b/ras_commander/hdf/HdfResultsPlan.py @@ -32,7 +32,7 @@ All methods are static and designed to be used without class instantiation. """ -from typing import Dict, List, Optional, Tuple, Union +from typing import Dict, List, Optional, Tuple from pathlib import Path import h5py import pandas as pd @@ -41,7 +41,6 @@ from ..LoggingConfig import get_logger import numpy as np from datetime import datetime -from ..RasPrj import ras logger = get_logger(__name__) @@ -211,7 +210,7 @@ def get_runtime_data(hdf_path: Path) -> Optional[pd.DataFrame]: """ try: if hdf_path is None: - logger.error(f"Could not find HDF file for input") + logger.error("Could not find HDF file for input") return None with h5py.File(hdf_path, 'r') as hdf_file: @@ -325,7 +324,14 @@ def get_runtime_data(hdf_path: Path) -> Optional[pd.DataFrame]: @standardize_input(file_type='plan_hdf') def get_reference_timeseries(hdf_path: Path, reftype: str) -> pd.DataFrame: """ - Get reference line or point timeseries output from HDF file. + Get reference-line or reference-point time series as a tidy DataFrame. + + Modern HEC-RAS HDF files store a ``Name`` vector plus shared + ``(time, feature)`` matrices such as Flow, Velocity, and Water Surface. + This compatibility method delegates schema handling to + :class:`HdfResultsXsec` and converts its native xarray Dataset to one + row per time/feature pair. It therefore supports all numeric variables + present in the HDF instead of assuming one dataset per feature. Args: hdf_path (Path): Path to HEC-RAS plan HDF file @@ -333,31 +339,29 @@ def get_reference_timeseries(hdf_path: Path, reftype: str) -> pd.DataFrame: ras_object (RasPrj, optional): Specific RAS object to use. If None, uses the global ras instance. Returns: - pd.DataFrame: DataFrame containing reference timeseries data + pd.DataFrame: Tidy frame containing time, reference ID/name, mesh + name, and every numeric native result variable. Returns an + empty frame when the requested reference output is absent. """ try: - with h5py.File(hdf_path, 'r') as hdf_file: - base_path = "Results/Unsteady/Output/Output Blocks/Base Output/Unsteady Time Series" - ref_path = f"{base_path}/Reference {reftype.capitalize()}" - - if ref_path not in hdf_file: - logger.debug(f"Reference {reftype} data not found in HDF file") - return pd.DataFrame() - - ref_group = hdf_file[ref_path] - time_data = hdf_file[f"{base_path}/Time"][:] - - dfs = [] - for ref_name in ref_group.keys(): - ref_data = ref_group[ref_name][:] - df = pd.DataFrame(ref_data, columns=[ref_name]) - df['Time'] = time_data - dfs.append(df) + from .HdfResultsXsec import HdfResultsXsec - if not dfs: - return pd.DataFrame() + normalized = str(reftype).strip().lower() + if normalized == "lines": + dataset = HdfResultsXsec.get_ref_lines_timeseries(hdf_path) + elif normalized == "points": + dataset = HdfResultsXsec.get_ref_points_timeseries(hdf_path) + else: + raise ValueError("reftype must be 'lines' or 'points'") - return pd.concat(dfs, axis=1) + if not dataset.data_vars: + logger.debug( + "Reference %s data not found in %s", + normalized, + Path(hdf_path).name, + ) + return pd.DataFrame() + return dataset.to_dataframe().reset_index() except Exception as e: logger.error(f"Error reading reference {reftype} timeseries: {str(e)}") @@ -958,7 +962,7 @@ def get_compute_messages(hdf_path: Path) -> str: txt_contents = RasControl.get_comp_msgs(hdf_path) if txt_contents: logger.debug( - f"HDF file not found, successfully retrieved computation messages from .txt file" + "HDF file not found, successfully retrieved computation messages from .txt file" ) return txt_contents except Exception as e: @@ -976,7 +980,7 @@ def get_compute_messages(hdf_path: Path) -> str: txt_contents = RasControl.get_comp_msgs(hdf_path) if txt_contents: logger.debug( - f"HDF extraction failed, successfully retrieved computation messages from .txt file" + "HDF extraction failed, successfully retrieved computation messages from .txt file" ) return txt_contents except Exception as fallback_error: diff --git a/setup.py b/setup.py index 507fd1264..286a99c59 100644 --- a/setup.py +++ b/setup.py @@ -128,7 +128,11 @@ def run(self): 'azure-mgmt-compute>=30.0', ], # GUI automation and screenshot capture (Windows only) - 'gui': ['Pillow>=9.0', 'comtypes>=1.4.0; sys_platform == "win32"'], + 'gui': [ + 'Pillow>=9.0', + 'psutil>=5.6.6', + 'comtypes>=1.4.0; sys_platform == "win32"', + ], # USGS gauge data integration 'usgs': ['dataretrieval>=1.0'], # Precipitation enhancements diff --git a/tests/test_geom_bc_lines.py b/tests/test_geom_bc_lines.py index bd97d6b2a..487a4e182 100644 --- a/tests/test_geom_bc_lines.py +++ b/tests/test_geom_bc_lines.py @@ -19,7 +19,6 @@ import re from pathlib import Path -import numpy as np import pytest @@ -97,6 +96,24 @@ def test_add_single_bc_line(self, skeleton_geom): # Existing block survives unchanged. assert "BC Line Name=Existing" in text + def test_add_preserves_existing_backup_and_uses_next_number(self, skeleton_geom): + from ras_commander import GeomBcLines + + original_backup = skeleton_geom.with_suffix(".g01.bak") + original_backup.write_text("parent geometry backup", encoding="utf-8") + + result = GeomBcLines.add_bc_lines( + skeleton_geom, + lines=[{ + "name": "NewBC", + "storage_area": "Perimeter 1", + "coordinates": [(0, 0), (10, 0)], + }], + ) + + assert original_backup.read_text(encoding="utf-8") == "parent geometry backup" + assert Path(result["backup_path"]).name == "project.g01.bak1" + def test_added_block_groups_with_existing_bc_lines(self, skeleton_geom): """New BC line is inserted after the LAST `BC Line Text Position=` so all BC lines stay contiguous in the file (matches HEC-RAS diff --git a/tests/test_geom_mesh.py b/tests/test_geom_mesh.py index 2bc6b4d41..f4ce9af3d 100644 --- a/tests/test_geom_mesh.py +++ b/tests/test_geom_mesh.py @@ -28,6 +28,7 @@ _dedupe_seed_points = geom_mesh_module._dedupe_seed_points _bad_seed_indexes = geom_mesh_module._bad_seed_indexes _remove_seed_indexes = geom_mesh_module._remove_seed_indexes +_seed_indexes_outside_perimeter = geom_mesh_module._seed_indexes_outside_perimeter _safe_non_virtual_cell_count = geom_mesh_module._safe_non_virtual_cell_count HECRAS_INTEGRATION_ENV = "RAS_COMMANDER_RUN_HECRAS_INTEGRATION" @@ -337,13 +338,39 @@ def _mock_generate_success(monkeypatch, geom_text_path: Path, *, has_breaklines: "spacing_dy": 50.0, "cell_count": 0, "dataset_count": 2, - } + }, + { + "name": "SecondaryArea", + "spacing_dx": 50.0, + "spacing_dy": 50.0, + "cell_count": 0, + "dataset_count": 2, + }, ], ) hdf_mtime = geom_text_path.stat().st_mtime + 1.0 os.utime(hdf_path, (hdf_mtime, hdf_mtime)) monkeypatch.setattr(geom_mesh_module, "_load_dlls", lambda hecras_dir=None: None) + monkeypatch.setattr( + geom_mesh_module, + "_audit_domain_containment_hdf", + lambda hdf_path, mesh_name, base_cell_spacing: ( + geom_mesh_module.DomainContainmentResult( + mesh_name=mesh_name, + geom_hdf_path=str(hdf_path), + base_cell_spacing=float(base_cell_spacing), + inward_buffer_distance=float(base_cell_spacing), + admissible_geometry_type="Polygon", + checked_counts={ + "breakline": 0, + "refinement_region": 0, + "structure": 0, + }, + violations=[], + ) + ), + ) # Build mock .NET geometry chain mock_geom_obj = MagicMock() @@ -458,6 +485,184 @@ def _write_mesh_hdf(hdf_path: Path, areas): flow_areas.create_dataset("Attributes", data=np.array(records, dtype=dtype)) +def _write_2d_perimeter_pair(tmp_path: Path, hdf_coords=None): + text_coords = [(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)] + if hdf_coords is None: + hdf_coords = text_coords + geom_text = tmp_path / "perimeter.g01" + coordinate_lines = "".join( + f"{x_coord:16.6f}{y_coord:16.6f}\n" for x_coord, y_coord in text_coords + ) + geom_text.write_text( + "Geom Title=Perimeter Test\n" + "Storage Area=MainArea,5,5\n" + "Storage Area Surface Line= 5\n" + + coordinate_lines + + "Storage Area Type= 0\n" + "Storage Area Area=\n" + "Storage Area Min Elev=\n" + "Storage Area Is2D=-1\n" + "Storage Area Point Generation Data=50,50,,\n" + "Storage Area 2D Points= 3\n", + encoding="utf-8", + ) + hdf_path = geom_text.with_suffix(".g01.hdf") + _write_mesh_hdf( + hdf_path, + [{"name": "MainArea", "spacing_dx": 50.0, "spacing_dy": 50.0, "cell_count": 0}], + ) + with h5py.File(hdf_path, "a") as hdf: + group = hdf["Geometry/2D Flow Areas"] + group.create_dataset( + "Polygon Info", + data=np.array([[0, len(hdf_coords), 0, 1]], dtype=np.int32), + ) + group.create_dataset( + "Polygon Parts", + data=np.array([[0, len(hdf_coords)]], dtype=np.int32), + ) + group.create_dataset( + "Polygon Points", + data=np.asarray(hdf_coords, dtype=np.float64), + ) + return geom_text, hdf_path + + +def _write_containment_fixture( + tmp_path: Path, + *, + breaklines=(), + refinement_regions=(), + structures=(), + bc_lines=(), +): + """Write a one-area geometry/HDF pair for domain-containment tests.""" + perimeter = np.array( + [[0.0, 0.0], [1000.0, 0.0], [1000.0, 1000.0], [0.0, 1000.0], [0.0, 0.0]], + dtype=np.float64, + ) + geom_text = tmp_path / "containment.g01" + text = ( + "Geom Title=Containment Test\n" + "Storage Area=MainArea,500,500\n" + "Storage Area Surface Line= 5\n" + + "".join(f"{x:16.6f}{y:16.6f}\n" for x, y in perimeter) + + "Storage Area Type= 0\n" + "Storage Area Area=\n" + "Storage Area Min Elev=\n" + "Storage Area Is2D=-1\n" + "Storage Area Point Generation Data=,,100.000000,100.000000\n" + "Storage Area 2D Points= 0 \n" + ) + for name, coords in breaklines: + text += ( + f"BreakLine Name={name}\n" + "BreakLine CellSize Min=25.0\n" + "BreakLine CellSize Max=100.0\n" + "BreakLine Near Repeats=0\n" + "BreakLine Protection Radius=0\n" + f"BreakLine Polyline= {len(coords)} \n" + + "".join(f"{x:16.6f}{y:16.6f}\n" for x, y in coords) + ) + geom_text.write_text(text, encoding="utf-8") + hdf_path = geom_text.with_suffix(".g01.hdf") + _write_mesh_hdf( + hdf_path, + [{"name": "MainArea", "spacing_dx": 100.0, "spacing_dy": 100.0, "cell_count": 0}], + ) + + def add_polylines(hf, path, features, dtype, info_name="Polyline Info", points_name="Polyline Points"): + if not features: + return + group = hf.create_group(path) + records = [] + info = [] + parts = [] + points = [] + for index, feature in enumerate(features): + name, coords, *extra = feature + start = len(points) + points.extend(coords) + info.append((start, len(coords), index, 1)) + parts.append((0, len(coords))) + records.append((name.encode(), *(value.encode() for value in extra))) + group.create_dataset("Attributes", data=np.array(records, dtype=dtype)) + group.create_dataset(info_name, data=np.array(info, dtype=np.int32)) + group.create_dataset("Polyline Parts", data=np.array(parts, dtype=np.int32)) + group.create_dataset(points_name, data=np.asarray(points, dtype=np.float64)) + + with h5py.File(hdf_path, "a") as hf: + areas = hf["Geometry/2D Flow Areas"] + areas.create_dataset("Polygon Info", data=np.array([[0, 5, 0, 1]], dtype=np.int32)) + areas.create_dataset("Polygon Parts", data=np.array([[0, 5]], dtype=np.int32)) + areas.create_dataset("Polygon Points", data=perimeter) + add_polylines( + hf, + "Geometry/2D Flow Area Break Lines", + breaklines, + np.dtype([("Name", "S64")]), + ) + if refinement_regions: + group = hf.create_group("Geometry/2D Flow Area Refinement Regions") + records = [] + info = [] + parts = [] + points = [] + for index, (name, coords) in enumerate(refinement_regions): + start = len(points) + points.extend(coords) + info.append((start, len(coords), index, 1)) + parts.append((0, len(coords))) + records.append((name.encode(),)) + group.create_dataset( + "Attributes", + data=np.array(records, dtype=np.dtype([("Name", "S64")])), + ) + group.create_dataset("Polygon Info", data=np.array(info, dtype=np.int32)) + group.create_dataset("Polygon Parts", data=np.array(parts, dtype=np.int32)) + group.create_dataset("Polygon Points", data=np.asarray(points, dtype=np.float64)) + if structures: + group = hf.create_group("Geometry/Structures") + records = [] + info = [] + points = [] + for name, coords, structure_type, upstream, downstream in structures: + start = len(points) + points.extend(coords) + info.append((start, len(coords))) + records.append( + ( + name.encode(), + structure_type.encode(), + upstream.encode(), + downstream.encode(), + ) + ) + group.create_dataset( + "Attributes", + data=np.array( + records, + dtype=np.dtype( + [ + ("Name", "S64"), + ("Type", "S32"), + ("US SA/2D", "S64"), + ("DS SA/2D", "S64"), + ] + ), + ), + ) + group.create_dataset("Centerline Info", data=np.array(info, dtype=np.int32)) + group.create_dataset("Centerline Points", data=np.asarray(points, dtype=np.float64)) + add_polylines( + hf, + "Geometry/Boundary Condition Lines", + bc_lines, + np.dtype([("Name", "S64"), ("SA-2D", "S64"), ("Type", "S16")]), + ) + return geom_text, hdf_path + + def _install_fake_rasmapper_scripting(monkeypatch): """Install a stub RasMapperLib.Scripting module for helper tests.""" @@ -692,6 +897,34 @@ def test_ensure_hdf_accepts_matching_content(self, breakline_geom_text): assert _ensure_hdf(breakline_geom_text) == hdf_path + def test_generate_gate_ignores_old_seeds_but_not_stale_perimeter(self, tmp_path): + geom_text, hdf_path = _write_2d_perimeter_pair(tmp_path) + + strict = geom_mesh_module._mesh_hdf_consistency_issues(geom_text, hdf_path) + regeneration = geom_mesh_module._mesh_hdf_consistency_issues( + geom_text, + hdf_path, + ignore_seed_count=True, + ) + + assert any("Storage Area 2D Points=3" in issue for issue in strict) + assert regeneration == [] + + stale_coords = [ + (0.0, 0.0), + (20.0, 0.0), + (20.0, 20.0), + (0.0, 20.0), + (0.0, 0.0), + ] + _geom_text, stale_hdf = _write_2d_perimeter_pair(tmp_path, stale_coords) + issues = geom_mesh_module._mesh_hdf_consistency_issues( + geom_text, + stale_hdf, + ignore_seed_count=True, + ) + assert any("does not match HDF perimeter" in issue for issue in issues) + def test_ensure_hdf_recompiles_missing_geometry_when_opted_in( self, monkeypatch, breakline_geom_text ): @@ -786,6 +1019,215 @@ def test_consistency_passes_when_breakline_counts_match(self, tmp_path): assert _ensure_hdf(geom_text) == hdf_path +class TestDomainContainmentAudit: + """One-cell inward containment is a mandatory pre-mesh gate.""" + + def test_audits_breaklines_refinements_and_target_structures(self, tmp_path): + geom_text, _hdf_path = _write_containment_fixture( + tmp_path, + breaklines=[ + ("inside", [(100.0, 100.0), (900.0, 900.0)]), + ("crossing", [(50.0, 500.0), (500.0, 500.0)]), + ("outside", [(20.0, 20.0), (80.0, 80.0)]), + ], + refinement_regions=[ + ( + "margin-region", + [ + (50.0, 300.0), + (150.0, 300.0), + (150.0, 400.0), + (50.0, 400.0), + (50.0, 300.0), + ], + ) + ], + structures=[ + ( + "inside-connection", + [(200.0, 200.0), (800.0, 200.0)], + "SA/2D Area Connection", + "MainArea", + "", + ), + ( + "crossing-connection", + [(50.0, 600.0), (500.0, 600.0)], + "SA/2D Area Connection", + "MainArea", + "", + ), + ( + "other-area-connection", + [(-1000.0, -1000.0), (-900.0, -900.0)], + "SA/2D Area Connection", + "OtherArea", + "", + ), + ], + # External BC lines are intentionally outside this gate. + bc_lines=[ + ( + "external-bc", + [(-500.0, 0.0), (-500.0, 1000.0)], + "MainArea", + "External", + ) + ], + ) + + result = GeomMesh.audit_domain_containment( + geom_text, + mesh_name="MainArea", + ) + + assert result.base_cell_spacing == 100.0 + assert result.inward_buffer_distance == 100.0 + assert result.checked_counts == { + "breakline": 3, + "refinement_region": 1, + "structure": 2, + } + assert not result.ok + assert len(result.violations) == 4 + assert { + (item.feature_type, item.feature_name) + for item in result.violations + } == { + ("breakline", "crossing"), + ("breakline", "outside"), + ("refinement_region", "margin-region"), + ("structure", "crossing-connection"), + } + assert all( + item.reason == "outside_one_cell_inward_buffer" + for item in result.violations + ) + evidence = result.to_dict() + assert evidence["violation_count"] == 4 + assert evidence["ok"] is False + assert all(not hasattr(value, "geom_type") for value in evidence.values()) + + def test_boundary_condition_lines_do_not_fail_mesh_feature_gate(self, tmp_path): + geom_text, _hdf_path = _write_containment_fixture( + tmp_path, + bc_lines=[ + ( + "external-bc", + [(-500.0, 0.0), (-500.0, 1000.0)], + "MainArea", + "External", + ) + ], + ) + + result = GeomMesh.audit_domain_containment(geom_text) + + assert result.ok + assert result.checked_counts == { + "breakline": 0, + "refinement_region": 0, + "structure": 0, + } + + def test_generate_refuses_before_loading_native_mesher(self, monkeypatch, tmp_path): + geom_text, _hdf_path = _write_containment_fixture( + tmp_path, + breaklines=[("too-close", [(50.0, 500.0), (500.0, 500.0)])], + ) + loaded = [] + monkeypatch.setattr( + geom_mesh_module, + "_load_dlls", + lambda hecras_dir=None: loaded.append(hecras_dir), + ) + + result = GeomMesh.generate(geom_text, mesh_name="MainArea") + + assert not result.ok + assert "Pre-mesh domain containment failed" in result.error_message + assert result.domain_containment is not None + assert result.domain_containment.violations[0].feature_name == "too-close" + assert loaded == [] + + def test_malformed_structure_collection_fails_closed(self, tmp_path): + geom_text, hdf_path = _write_containment_fixture(tmp_path) + with h5py.File(hdf_path, "a") as hf: + structures = hf.create_group("Geometry/Structures") + structures.create_dataset( + "Attributes", + data=np.array([(b"broken",)], dtype=np.dtype([("Name", "S32")])), + ) + structures.create_dataset( + "Centerline Info", + data=np.array([[0, 2]], dtype=np.int32), + ) + + with pytest.raises(RuntimeError, match="Centerline Points"): + GeomMesh.audit_domain_containment(geom_text) + + def test_native_zero_count_structure_group_is_empty(self, tmp_path): + geom_text, hdf_path = _write_containment_fixture(tmp_path) + with h5py.File(hdf_path, "a") as hf: + structures = hf.create_group("Geometry/Structures") + structures.attrs["Bridge/Culvert Count"] = 0 + structures.attrs["Connection Count"] = 0 + structures.attrs["Inline Structure Count"] = 0 + structures.attrs["Lateral Structure Count"] = 0 + + result = GeomMesh.audit_domain_containment(geom_text) + + assert result.ok + assert result.checked_counts["structure"] == 0 + + def test_native_multipart_collection_accepts_absolute_part_indexes(self, tmp_path): + _geom_text, hdf_path = _write_containment_fixture(tmp_path) + with h5py.File(hdf_path, "a") as hf: + group = hf.create_group("Geometry/2D Flow Area Break Lines") + group.create_dataset( + "Attributes", + data=np.array( + [(b"single",), (b"multipart",)], + dtype=np.dtype([("Name", "S32")]), + ), + ) + group.create_dataset( + "Polyline Info", + data=np.array([[0, 2, 0, 1], [2, 4, 1, 2]], dtype=np.int32), + ) + group.create_dataset( + "Polyline Parts", + # The multipart feature uses absolute offsets 2 and 4. + data=np.array([[0, 2], [2, 2], [4, 2]], dtype=np.int32), + ) + group.create_dataset( + "Polyline Points", + data=np.array( + [ + [200.0, 200.0], + [300.0, 300.0], + [200.0, 400.0], + [300.0, 400.0], + [600.0, 700.0], + [700.0, 700.0], + ], + dtype=np.float64, + ), + ) + + with h5py.File(hdf_path, "r") as hf: + features = geom_mesh_module._hdf_polyline_features( + hf, + "Geometry/2D Flow Area Break Lines", + ) + + multipart = features[1]["geometry"] + assert multipart.geom_type == "MultiLineString" + assert len(multipart.geoms) == 2 + assert list(multipart.geoms[0].coords)[0] == (200.0, 400.0) + assert list(multipart.geoms[1].coords)[0] == (600.0, 700.0) + + class TestGeometryAssociation: """Test geometry HDF association API.""" @@ -967,6 +1409,17 @@ def test_remove_seed_indexes_preserves_order(self): (3.0, 3.0), ] + def test_seed_outside_filter_preserves_exact_perimeter(self): + seeds = MockPointMs() + seeds.Add(MockPointM(5.0, 5.0)) + seeds.Add(MockPointM(10.0, 5.0)) + seeds.Add(MockPointM(10.1, 5.0)) + perimeter = MockPolygon( + [(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0)] + ) + + assert _seed_indexes_outside_perimeter(seeds, perimeter) == {2} + def test_defaults_persist_geometry_with_existing_hdf(self, monkeypatch, breakline_geom_text): captured = _mock_generate_success( monkeypatch, @@ -1031,7 +1484,7 @@ def test_generate_internal_progress_logs_debug_not_info( ) def test_accepts_breakline_spacing_override(self, monkeypatch, breakline_geom_text): - captured = _mock_generate_success( + _mock_generate_success( monkeypatch, breakline_geom_text, has_breaklines=True, @@ -1049,7 +1502,7 @@ def test_accepts_breakline_spacing_override(self, monkeypatch, breakline_geom_te assert "BreakLine CellSize Min=75.000000" in text assert "BreakLine CellSize Max=125.000000" in text - def test_perimeter_fix_reports_external_hdf_requirement( + def test_short_segment_suggestion_preserves_authored_perimeter( self, monkeypatch, breakline_geom_text ): captured = _mock_generate_success( @@ -1104,9 +1557,9 @@ def fake_generate_seeds_via_net(hdf_path, ns, fid=0): result = GeomMesh.generate(breakline_geom_text, mesh_index=1) - assert result.ok is False - assert result.status == "exception" - assert "cannot generate .g##.hdf" in result.error_message + assert result.ok is True + assert result.status == "complete" + assert result.fixes_applied == [] assert seed_calls == [1] @pytest.mark.parametrize( @@ -1455,6 +1908,81 @@ def test_utf8_truncation_preserves_valid_bytes(self, refinement_region_hdf): stored.encode("utf-8") assert "\ufffd" not in stored + def test_replace_collection_is_atomic_and_preserves_other_hdf_content( + self, refinement_region_hdf + ): + geom_text, hdf_path = refinement_region_hdf + with h5py.File(str(hdf_path), "a") as hf: + hf.create_dataset("Geometry/Sentinel", data=np.array([7, 8, 9])) + + backup = GeomMesh.replace_refinement_regions( + geom_text, + [ + { + "name": "Clipped", + "polygon": [ + (25.0, 25.0), + (75.0, 25.0), + (75.0, 75.0), + (25.0, 75.0), + ], + "spacing_dx": 20.0, + "spacing_dy": 10.0, + } + ], + expected_existing_names=["North", "North", ""], + ) + + assert backup == hdf_path.with_suffix(".hdf.bak") + assert backup.is_file() + regions = GeomMesh.get_refinement_regions(geom_text) + assert regions == [ + { + "fid": 0, + "name": "Clipped", + "spacing_dx": 20.0, + "spacing_dy": 10.0, + } + ] + with h5py.File(str(hdf_path), "r") as hf: + assert np.array_equal(hf["Geometry/Sentinel"][:], [7, 8, 9]) + points = hf[ + "Geometry/2D Flow Area Refinement Regions/Polygon Points" + ][:] + assert len(points) == 5 + assert np.allclose(points[0], points[-1]) + + def test_replace_collection_guard_failure_preserves_hdf( + self, refinement_region_hdf + ): + geom_text, hdf_path = refinement_region_hdf + original = hdf_path.read_bytes() + + with pytest.raises(ValueError, match="collection changed"): + GeomMesh.replace_refinement_regions( + geom_text, + [], + expected_existing_names=["wrong"], + ) + + assert hdf_path.read_bytes() == original + assert not hdf_path.with_suffix(".hdf.bak").exists() + + def test_replace_collection_with_empty_removes_group( + self, refinement_region_hdf + ): + geom_text, hdf_path = refinement_region_hdf + + GeomMesh.replace_refinement_regions( + geom_text, + [], + expected_existing_names=["North", "North", ""], + create_backup=False, + ) + + with h5py.File(str(hdf_path), "r") as hf: + assert "Geometry/2D Flow Area Refinement Regions" not in hf + # ── Add refinement region tests ─────────────────────────────────── diff --git a/tests/test_geom_reference_features.py b/tests/test_geom_reference_features.py index c8074b22d..4774a8ef2 100644 --- a/tests/test_geom_reference_features.py +++ b/tests/test_geom_reference_features.py @@ -205,3 +205,94 @@ def test_generated_reference_lines_match_hdf_reference_line_schema(tmp_path): LineString(generated[1]["coordinates"]), tolerance=1.0e-9, ) + + +def _reference_line_fixture(tmp_path): + geom_file = tmp_path / "Model.g01" + geom_file.write_text( + "Geom Title=Test\r\n" + "Reference Line Name=Keep A \r\n" + "Reference Line Storage Area=Mesh 1 \r\n" + "Reference Line Start Position= 0 , 0 \r\n" + "Reference Line Middle Position= 5 , 0 \r\n" + "Reference Line End Position= 10 , 0 \r\n" + "Reference Line Arc= 2 \r\n" + " 0 0 10 0\r\n" + "Reference Line Text Position= 1 , 1 \r\n" + "Reference Line Name=Outside \r\n" + "Reference Line Storage Area=Mesh 1 \r\n" + "Reference Line Start Position= 20 , 0 \r\n" + "Reference Line Middle Position= 25 , 0 \r\n" + "Reference Line End Position= 30 , 0 \r\n" + "Reference Line Arc= 2 \r\n" + " 20 0 30 0\r\n" + "Reference Line Text Position= 1 , 1 \r\n" + "Reference Line Name=Other Area \r\n" + "Reference Line Storage Area=Mesh 2 \r\n" + "Reference Line Start Position= 0 , 1 \r\n" + "Reference Line Middle Position= 5 , 1 \r\n" + "Reference Line End Position= 10 , 1 \r\n" + "Reference Line Arc= 2 \r\n" + " 0 1 10 1\r\n" + "Reference Line Text Position= 1 , 1 \r\n" + "LCMann TimeDateStamp=01JAN2026 0000\r\n", + encoding="utf-8", + newline="", + ) + return geom_file + + +def test_replace_reference_lines_is_area_scoped_and_preserves_crlf(tmp_path): + geom_file = _reference_line_fixture(tmp_path) + + result = GeomReferenceFeatures.replace_reference_lines( + geom_file, + lines=[{"name": "Keep A", "coordinates": [(1.0, 0.0), (9.0, 0.0)]}], + storage_area="Mesh 1", + expected_existing_names=["Keep A", "Outside"], + ) + + assert result["removed"] == ["Keep A", "Outside"] + assert result["inserted"] == ["Keep A"] + parsed = GeomReferenceFeatures.get_reference_lines(geom_file) + assert [item["name"] for item in parsed] == ["Keep A", "Other Area"] + assert np.allclose(parsed[0]["coordinates"], [(1.0, 0.0), (9.0, 0.0)]) + assert b"\r\n" in geom_file.read_bytes() + assert b"\n" not in geom_file.read_bytes().replace(b"\r\n", b"") + + +def test_replace_reference_lines_preserves_existing_backup(tmp_path): + geom_file = _reference_line_fixture(tmp_path) + original = geom_file.read_bytes() + first_backup = geom_file.with_suffix(".g01.bak") + first_backup.write_bytes(b"pre-existing evidence") + + result = GeomReferenceFeatures.replace_reference_lines( + geom_file, + lines=[], + storage_area="Mesh 1", + expected_existing_names=["Keep A", "Outside"], + ) + + backup = Path(result["backup_path"]) + assert first_backup.read_bytes() == b"pre-existing evidence" + assert backup.name == "Model.g01.bak1" + assert backup.read_bytes() == original + assert [ + item["name"] for item in GeomReferenceFeatures.get_reference_lines(geom_file) + ] == ["Other Area"] + + +def test_replace_reference_lines_expected_names_fail_closed(tmp_path): + geom_file = _reference_line_fixture(tmp_path) + before = geom_file.read_bytes() + + with pytest.raises(ValueError, match="population changed"): + GeomReferenceFeatures.replace_reference_lines( + geom_file, + lines=[], + storage_area="Mesh 1", + expected_existing_names=["Keep A"], + ) + + assert geom_file.read_bytes() == before diff --git a/tests/test_geom_storage_2d_flow_area_writer.py b/tests/test_geom_storage_2d_flow_area_writer.py index d46536a69..ed176dbf4 100644 --- a/tests/test_geom_storage_2d_flow_area_writer.py +++ b/tests/test_geom_storage_2d_flow_area_writer.py @@ -11,9 +11,9 @@ sys.path.insert(0, str(repo_root)) shapely = pytest.importorskip("shapely") -from shapely.geometry import Polygon +from shapely.geometry import Polygon # noqa: E402 -from ras_commander.geom import GeomParser, GeomStorage +from ras_commander.geom import GeomParser, GeomStorage # noqa: E402 def _format_xy_rows(points, *, values_per_line): @@ -421,3 +421,108 @@ def test_surface_line_fields_never_exceed_16_chars(): for i in range(0, len(raw), 16): field = raw[i:i+16] assert len(field) <= 16, f"Field {field!r} exceeds 16 chars" + + +def _breakline_replacement_fixture(tmp_path): + area_coords = [ + (0.0, 0.0), + (100.0, 0.0), + (100.0, 100.0), + (0.0, 100.0), + (0.0, 0.0), + ] + return _write_geom_file( + tmp_path, + [ + "Geom Title=Breakline Replacement\n", + "Program Version=6.60\n", + "Storage Area=Mesh,50.0000000,50.0000000\n", + "Storage Area Surface Line= 5\n", + *_format_xy_rows(area_coords, values_per_line=2), + "Storage Area Type= 0\n", + "Storage Area Area=\n", + "Storage Area Min Elev=\n", + "Storage Area Is2D=-1\n", + "Storage Area Point Generation Data=10,10,,\n", + "Storage Area 2D Points= 0\n", + "Storage Area 2D PointsPerimeterTime=01Jan2026 00:00:00\n", + "Storage Area Mannings=0.04\n", + *GeomStorage._format_breakline_block( + "Retain", [(10.0, 10.0), (90.0, 90.0)], 5.0, 10.0 + ), + *GeomStorage._format_breakline_block( + "Remove", [(10.0, 90.0), (90.0, 10.0)], 5.0, 10.0 + ), + "BC Line Name=Boundary\n", + ], + ) + + +def test_replace_breaklines_atomically_replaces_complete_collection(tmp_path): + geom_file = _breakline_replacement_fixture(tmp_path) + + backup = GeomStorage.replace_breaklines( + geom_file, + "Mesh", + [ + { + "name": "Retain", + "coords": [(10.0, 10.0), (50.0, 50.0)], + "cell_size_near": 4.0, + "cell_size_far": 8.0, + "near_repeats": 2, + "protection_radius": 1, + }, + { + "name": "Clipped", + "coords": [(50.0, 50.0), (95.0, 50.0)], + "cell_size_near": 3.0, + "cell_size_far": 6.0, + }, + ], + expected_existing_names=["Retain", "Remove"], + ) + + text = geom_file.read_text(encoding="utf-8") + assert backup == geom_file.with_suffix(".g01.bak") + assert backup.is_file() + assert text.count("BreakLine Name=") == 2 + assert "BreakLine Name=Remove" not in text + assert text.index("BreakLine Name=Retain") < text.index("BreakLine Name=Clipped") + assert text.index("BreakLine Name=Clipped") < text.index("BC Line Name=Boundary") + assert "BreakLine CellSize Min=4.0" in text + assert "BreakLine Near Repeats=2" in text + assert "BreakLine Protection Radius=1" in text + + +def test_replace_breaklines_guard_failure_preserves_original(tmp_path): + geom_file = _breakline_replacement_fixture(tmp_path) + original = geom_file.read_bytes() + + with pytest.raises(ValueError, match="collection changed"): + GeomStorage.replace_breaklines( + geom_file, + "Mesh", + [], + expected_existing_names=["wrong"], + ) + + assert geom_file.read_bytes() == original + assert not geom_file.with_suffix(".g01.bak").exists() + + +def test_replace_breaklines_rejects_duplicate_names_before_write(tmp_path): + geom_file = _breakline_replacement_fixture(tmp_path) + original = geom_file.read_bytes() + + with pytest.raises(ValueError, match="Duplicate breakline name"): + GeomStorage.replace_breaklines( + geom_file, + "Mesh", + [ + {"name": "same", "coords": [(0.0, 0.0), (1.0, 1.0)]}, + {"name": "SAME", "coords": [(1.0, 0.0), (0.0, 1.0)]}, + ], + ) + + assert geom_file.read_bytes() == original diff --git a/tests/test_gui_mesh_regeneration_exact.py b/tests/test_gui_mesh_regeneration_exact.py new file mode 100644 index 000000000..ad55797ce --- /dev/null +++ b/tests/test_gui_mesh_regeneration_exact.py @@ -0,0 +1,451 @@ +"""Regression coverage for exact-geometry RASMapper mesh regeneration.""" + +from pathlib import Path +from types import SimpleNamespace + +import h5py +import numpy as np +import pandas as pd +import pytest +from shapely.geometry import Polygon + +from ras_commander.RasMap import RasMap +from ras_commander.hdf.HdfMesh import HdfMesh +from ras_commander.gui.workflows.mesh_regeneration import ( + MeshRegenerationWorkflow, + _begin_geometry_hdf_transaction, + _capture_geometry_association_paths, + _finish_geometry_hdf_transaction, + _perimeter_validation, + _prepare_geometry_refresh_context, + _resolve_geometry_target, + _restore_geometry_association, + _supervise_owned_process_exit, + _validate_geometry_import, + _validate_geometry_refresh, +) + + +def _fake_project(tmp_path): + project_file = tmp_path / "Model.prj" + project_file.write_text( + "Proj Title=Model\n" + "Current Plan=p09\n" + "Geom File=g01\nGeom File=g03\n" + "Plan File=p01\nPlan File=p09\n", + encoding="utf-8", + ) + for number in ("01", "03"): + project_file.with_suffix(f".g{number}").write_text( + f"Geom Title=Geometry {number}\n", encoding="utf-8" + ) + project = SimpleNamespace( + initialized=True, + prj_file=project_file, + project_folder=tmp_path, + project_name="Model", + plan_df=pd.DataFrame( + [ + {"plan_number": "01", "geometry_number": "01"}, + {"plan_number": "09", "geometry_number": "03"}, + ] + ), + geom_df=pd.DataFrame( + [ + {"geom_number": "01", "geom_title": "Geometry 01"}, + {"geom_number": "03", "geom_title": "Geometry 03"}, + ] + ), + ) + project.check_initialized = lambda: None + return project + + +def _mapper_geometries(): + return [ + {"name": "Geometry 01", "geom_number": "01"}, + {"name": "Geometry 03", "geom_number": "03"}, + ] + + +def _write_mesh_hdf(path, polygon, *, area_name="Mesh"): + attr_dtype = np.dtype([("Name", "S32")]) + with h5py.File(path, "w") as hdf: + hdf.create_dataset( + "Geometry/2D Flow Areas/Attributes", + data=np.array([(area_name.encode("utf-8"),)], dtype=attr_dtype), + ) + points = np.asarray(polygon.exterior.coords) + hdf.create_dataset( + "Geometry/2D Flow Areas/Polygon Info", + data=np.array([[0, len(points), 0, 1]], dtype=np.int32), + ) + hdf.create_dataset( + "Geometry/2D Flow Areas/Polygon Parts", + data=np.array([[0, len(points)]], dtype=np.int32), + ) + hdf.create_dataset("Geometry/2D Flow Areas/Polygon Points", data=points) + base = f"Geometry/2D Flow Areas/{area_name}" + hdf.create_dataset(f"{base}/Perimeter", data=np.asarray(polygon.exterior.coords)) + hdf.create_dataset( + f"{base}/Cells Center Coordinate", + data=np.array([[1.0, 1.0], [5.0, 5.0], [9.0, 9.0]]), + ) + hdf.create_dataset( + f"{base}/FacePoints Coordinate", + data=np.array([[0.0, 0.0], [10.0, 0.0], [10.0, 10.0]]), + ) + hdf.create_dataset( + f"{base}/Faces FacePoint Indexes", + data=np.array([[0, 1], [1, 2], [2, 0]], dtype=np.int32), + ) + + +def test_default_geometry_comes_from_current_plan_not_first_registration( + tmp_path, monkeypatch +): + project = _fake_project(tmp_path) + monkeypatch.setattr( + RasMap, "list_geometries", staticmethod(lambda _ras: _mapper_geometries()) + ) + + target = _resolve_geometry_target(project) + + assert target["geom_number"] == "03" + assert target["geometry_name"] == "Geometry 03" + assert target["geom_file"] == tmp_path / "Model.g03" + + +def test_number_and_name_must_identify_same_mapper_geometry(tmp_path, monkeypatch): + project = _fake_project(tmp_path) + monkeypatch.setattr( + RasMap, "list_geometries", staticmethod(lambda _ras: _mapper_geometries()) + ) + + with pytest.raises(ValueError, match="is named"): + _resolve_geometry_target( + project, + geom_number="03", + geometry_name="Geometry 01", + ) + + +def test_duplicate_mapper_tree_names_are_rejected(tmp_path, monkeypatch): + project = _fake_project(tmp_path) + duplicate_names = [ + {"name": "Duplicate", "geom_number": "01"}, + {"name": "Duplicate", "geom_number": "03"}, + ] + monkeypatch.setattr( + RasMap, "list_geometries", staticmethod(lambda _ras: duplicate_names) + ) + + with pytest.raises(ValueError, match="not unique"): + _resolve_geometry_target(project, geom_number="03") + + +def test_refresh_context_edits_exact_geometry_root(tmp_path, monkeypatch): + project = _fake_project(tmp_path) + polygon = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)]) + monkeypatch.setattr( + "ras_commander.gui.workflows.mesh_regeneration._resolve_geometry_target", + lambda *_args, **_kwargs: { + "geom_number": "03", + "geometry_name": "Geometry 03", + "geom_file": tmp_path / "Model.g03", + "geom_hdf": tmp_path / "Model.g03.hdf", + }, + ) + monkeypatch.setattr( + "ras_commander.gui.workflows.mesh_regeneration._select_text_flow_area", + lambda *_args, **_kwargs: ("Mesh", polygon), + ) + monkeypatch.setattr( + "ras_commander.gui.workflows.mesh_regeneration._geometry_hdf_stats", + lambda *_args, **_kwargs: {}, + ) + monkeypatch.setattr( + "ras_commander.gui.workflows.mesh_regeneration._perimeter_validation", + lambda *_args, **_kwargs: {"valid": False, "error": "missing"}, + ) + + context = _prepare_geometry_refresh_context( + project, + geom_number="03", + geometry_name="Geometry 03", + flow_area_name="Mesh", + coordinate_tolerance=None, + ) + + assert context["target_path"] == ["Geometry 03"] + + +def test_perimeter_validation_rejects_stale_compiled_geometry(tmp_path): + expected = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)]) + stale = Polygon([(0, 0), (20, 0), (20, 20), (0, 20)]) + hdf_path = tmp_path / "Model.g03.hdf" + _write_mesh_hdf(hdf_path, stale) + + result = _perimeter_validation(hdf_path, "Mesh", expected) + + assert result["valid"] is False + assert "does not match" in result["error"] + assert result["hdf_area"] == pytest.approx(400.0) + + +def test_refresh_validation_proves_exact_target_and_other_hdf_isolation(tmp_path): + polygon = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)]) + target_hdf = tmp_path / "Model.g03.hdf" + other_hdf = tmp_path / "Model.g01.hdf" + _write_mesh_hdf(target_hdf, polygon) + _write_mesh_hdf(other_hdf, polygon) + project = _fake_project(tmp_path) + other_stat = other_hdf.stat() + + result = _validate_geometry_refresh( + { + "geom_number": "03", + "geometry_name": "Geometry 03", + "geom_file": tmp_path / "Model.g03", + "geom_hdf": target_hdf, + "flow_area_name": "Mesh", + "expected_polygon": polygon, + "coordinate_tolerance": None, + "ras_object": project, + "pre_hdf_stats": { + str(other_hdf.resolve()): (other_stat.st_size, other_stat.st_mtime_ns), + str(target_hdf.resolve()): (0, 0), + }, + "pre_perimeter_validation": {"valid": False, "error": "stale"}, + } + ) + + assert result["geom_number"] == "03" + assert result["post_perimeter"]["valid"] is True + assert result["mesh"]["valid"] is True + assert result["other_geometry_hdfs_unchanged"] is True + + +def test_import_validation_does_not_require_computation_cells(tmp_path): + polygon = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)]) + target_hdf = tmp_path / "Model.g03.hdf" + other_hdf = tmp_path / "Model.g01.hdf" + _write_mesh_hdf(target_hdf, polygon) + _write_mesh_hdf(other_hdf, polygon) + with h5py.File(target_hdf, "a") as hdf: + del hdf["Geometry/2D Flow Areas/Mesh"] + project = _fake_project(tmp_path) + other_stat = other_hdf.stat() + + result = _validate_geometry_import( + { + "geom_number": "03", + "geometry_name": "Geometry 03", + "geom_file": tmp_path / "Model.g03", + "geom_hdf": target_hdf, + "flow_area_name": "Mesh", + "expected_polygon": polygon, + "coordinate_tolerance": None, + "ras_object": project, + "pre_hdf_stats": { + str(other_hdf.resolve()): (other_stat.st_size, other_stat.st_mtime_ns), + str(target_hdf.resolve()): (0, 0), + }, + "pre_perimeter_validation": {"valid": False, "error": "missing"}, + } + ) + + assert result["post_perimeter"]["valid"] is True + assert "mesh" not in result + + +def test_geometry_refresh_captures_and_restores_associations(tmp_path, monkeypatch): + hdf_path = tmp_path / "Model.g03.hdf" + terrain = tmp_path / "Terrain" / "Terrain.hdf" + landcover = tmp_path / "Land Classification" / "LandCover.hdf" + terrain.parent.mkdir(parents=True) + landcover.parent.mkdir(parents=True) + terrain.write_bytes(b"terrain") + landcover.write_bytes(b"landcover") + with h5py.File(hdf_path, "w") as hdf: + geometry = hdf.create_group("Geometry") + geometry.attrs["Terrain Filename"] = b".\\Terrain\\Terrain.hdf" + geometry.attrs["Terrain Layername"] = b"Terrain" + geometry.attrs["Land Cover Filename"] = ( + b".\\Land Classification\\LandCover.hdf" + ) + geometry.attrs["Land Cover Layername"] = b"LandCover" + + captured = _capture_geometry_association_paths(hdf_path) + assert captured == { + "terrain_hdf_path": terrain.resolve(), + "landcover_hdf_path": landcover.resolve(), + } + + calls = [] + + def fake_set(target, **kwargs): + calls.append((target, kwargs)) + return Path(target) + + monkeypatch.setattr( + "ras_commander.geom.GeomMesh.set_geometry_association", + fake_set, + ) + monkeypatch.setattr( + "ras_commander.geom.GeomMesh.get_geometry_association", + lambda _target: {key: str(value) for key, value in captured.items()}, + ) + evidence = _restore_geometry_association( + { + "geom_hdf": hdf_path, + "ras_object": SimpleNamespace(), + "pre_geometry_association_paths": captured, + } + ) + + assert evidence["restored"] is True + assert calls[0][0] == hdf_path + assert calls[0][1]["terrain_hdf_path"] == terrain.resolve() + assert calls[0][1]["landcover_hdf_path"] == landcover.resolve() + + +def test_geometry_refresh_rejects_missing_association_artifact(tmp_path): + hdf_path = tmp_path / "Model.g03.hdf" + with h5py.File(hdf_path, "w") as hdf: + geometry = hdf.create_group("Geometry") + geometry.attrs["Terrain Filename"] = b".\\Terrain\\Missing.hdf" + geometry.attrs["Terrain Layername"] = b"Missing" + + with pytest.raises(FileNotFoundError, match="Cannot preserve terrain_hdf_path"): + _capture_geometry_association_paths(hdf_path) + + +def test_mesh_area_reader_supports_fresh_unmeshed_geometry_hdf(tmp_path): + polygon = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)]) + hdf_path = tmp_path / "Model.g03.hdf" + _write_mesh_hdf(hdf_path, polygon) + with h5py.File(hdf_path, "a") as hdf: + del hdf["Geometry/2D Flow Areas/Mesh"] + + areas = HdfMesh.get_mesh_areas(hdf_path) + + assert areas["mesh_name"].tolist() == ["Mesh"] + assert areas.iloc[0].geometry.area == pytest.approx(100.0) + + +def test_geometry_hdf_transaction_rolls_back_failed_import(tmp_path): + target = tmp_path / "Model.g03.hdf" + target.write_bytes(b"original") + context = {"geom_hdf": target} + + _begin_geometry_hdf_transaction(context) + target.write_bytes(b"failed replacement") + evidence = _finish_geometry_hdf_transaction( + context, + success=False, + keep_backup=False, + ) + + assert target.read_bytes() == b"original" + assert evidence["rolled_back"] is True + + +def test_geometry_hdf_transaction_commits_with_optional_backup(tmp_path): + target = tmp_path / "Model.g03.hdf" + target.write_bytes(b"original") + context = {"geom_hdf": target} + + _begin_geometry_hdf_transaction(context) + target.write_bytes(b"replacement") + evidence = _finish_geometry_hdf_transaction( + context, + success=True, + keep_backup=True, + ) + + assert target.read_bytes() == b"replacement" + backup = tmp_path / "Model.g03.hdf.pre-rasmapper.bak" + assert backup.read_bytes() == b"original" + assert evidence["backup"] == str(backup) + + +def test_geometry_hdf_transaction_numbers_existing_backup(tmp_path): + target = tmp_path / "Model.g03.hdf" + target.write_bytes(b"original") + existing_backup = tmp_path / "Model.g03.hdf.pre-rasmapper.bak" + existing_backup.write_bytes(b"earlier parent") + context = {"geom_hdf": target} + + _begin_geometry_hdf_transaction(context) + target.write_bytes(b"replacement") + evidence = _finish_geometry_hdf_transaction( + context, + success=True, + keep_backup=True, + ) + + next_backup = tmp_path / "Model.g03.hdf.pre-rasmapper.bak1" + assert target.read_bytes() == b"replacement" + assert existing_backup.read_bytes() == b"earlier parent" + assert next_backup.read_bytes() == b"original" + assert evidence["backup"] == str(next_backup) + + +def test_owned_process_supervision_terminates_only_captured_tree(monkeypatch): + import psutil + + class FakePopen: + pid = 101 + + def wait(self, timeout): + raise TimeoutError + + def poll(self): + return None + + class FakeOwned: + def __init__(self, pid): + self.pid = pid + self.alive = True + + def is_running(self): + return self.alive + + def status(self): + return "running" + + def terminate(self): + self.alive = False + + def kill(self): + self.alive = False + + owned = [FakeOwned(101), FakeOwned(202)] + + def fake_wait_procs(processes, timeout): + gone = [process for process in processes if not process.alive] + alive = [process for process in processes if process.alive] + return gone, alive + + monkeypatch.setattr(psutil, "wait_procs", fake_wait_procs) + + result = _supervise_owned_process_exit(FakePopen(), owned) + + assert result["observed_pids"] == [101, 202] + assert result["terminated_pids"] == [202, 101] + assert result["survivor_pids"] == [] + + +def test_single_attempt_steps_select_before_save_and_validate_after(tmp_path): + steps = MeshRegenerationWorkflow._build_single_attempt_steps( + {"timeout": 60, "close_after": True} + ) + names = [step.name for step in steps] + + assert names.index("Select exact geometry for editing") < names.index( + "Save geometry (trigger HDF regeneration)" + ) + assert names.index("Wait for save to complete") < names.index( + "Validate exact geometry HDF" + ) diff --git a/tests/test_hdf_bndry_logging.py b/tests/test_hdf_bndry_logging.py index 5c125c630..e96151a96 100644 --- a/tests/test_hdf_bndry_logging.py +++ b/tests/test_hdf_bndry_logging.py @@ -76,6 +76,48 @@ def _write_reference_lines_without_type_hdf(path: Path) -> None: ) +def _write_valid_breakline_attributes_hdf(path: Path) -> None: + attributes_dtype = np.dtype( + [ + ("Name", "S32"), + ("Cell Spacing Near", " Date: Sun, 30 Aug 2026 23:57:43 -0400 Subject: [PATCH 2/5] Add geometry-backed 2D boundary location authoring --- ras_commander/RasUnsteady.py | 402 ++++++++++++++++-- .../test_rasunsteady_2d_boundary_location.py | 298 +++++++++++++ 2 files changed, 664 insertions(+), 36 deletions(-) create mode 100644 tests/test_rasunsteady_2d_boundary_location.py diff --git a/ras_commander/RasUnsteady.py b/ras_commander/RasUnsteady.py index fbb0bdfaf..eaad3a472 100644 --- a/ras_commander/RasUnsteady.py +++ b/ras_commander/RasUnsteady.py @@ -67,6 +67,7 @@ def my_function(): - get_inline_hydrograph_boundaries() - Extract inline table BCs with time series data - inspect_boundary_blocks() - Inventory exact Boundary Location blocks in an owned stage - delete_boundary() - Preview/remove one exact block from an owned staged project +- ensure_2d_boundary_location() - Create an empty 2D boundary block after validating geometry - update_dss_run_identifier() - Update DSS path F-part for new scenarios - set_boundary_dss_link() - Convert inline BC to DSS-linked (complete state transition) - set_boundary_inline_hydrograph() - Write inline hydrograph, convert DSS to inline @@ -107,8 +108,9 @@ def my_function(): - set_non_newtonian_method() - Set the Non-Newtonian method by integer or name """ -import os +import math import numbers +import os import tempfile from datetime import datetime from pathlib import Path @@ -7312,6 +7314,262 @@ def find_line(prefix: str) -> Optional[int]: ) return True + @staticmethod + @log_call + def ensure_2d_boundary_location( + unsteady_file: Union[str, Path], + geometry_file: Union[str, Path], + *, + area_2d: str, + bc_line: str, + ras_object: Optional[Any] = None, + ) -> Dict[str, Any]: + """Ensure one geometry-backed 2D ``Boundary Location=`` block exists. + + This method creates only the empty location block. Follow it with a + boundary-type writer such as :meth:`set_boundary_inline_hydrograph` + or :meth:`set_normal_depth_boundary`. The geometry text file is + validated first so an unsteady-flow file cannot reference a missing + BC line or attach an existing line to the wrong 2D Flow Area. + + Parameters + ---------- + unsteady_file : str or Path + Unsteady-flow file path or number resolvable through + ``ras_object``. + geometry_file : str or Path + Exact plain-text geometry file containing the authored BC line. + area_2d : str + Exact ``BC Line Storage Area=`` value. + bc_line : str + Exact ``BC Line Name=`` value. + ras_object : optional + Project object used for short-number resolution and DataFrame + refresh. + + Returns + ------- + dict + Includes ``created``, the canonical location text, insertion + index, geometry validation evidence, and before/after boundary + counts. + + Raises + ------ + FileNotFoundError + If the explicit geometry file does not exist. + ValueError + If names are empty or exceed HEC-RAS fixed-field limits, the + geometry does not contain exactly one matching BC line, an + existing unsteady location conflicts, or the file header cannot + be identified safely. + + Notes + ----- + The emitted eight-field layout matches HEC-RAS 6.x 2D projects: + field 5 contains the 2D Flow Area and field 7 contains the BC line. + The operation is idempotent and uses atomic same-volume replacement. + """ + area_name = str(area_2d).strip() + line_name = str(bc_line).strip() + for value, label, maximum in ( + (area_name, "area_2d", 16), + (line_name, "bc_line", 32), + ): + if not value: + raise ValueError(f"{label} must be non-empty") + if any(character in value for character in (",", "\r", "\n")): + raise ValueError(f"{label} cannot contain commas or newlines") + if len(value) > maximum: + raise ValueError( + f"{label} exceeds the HEC-RAS fixed-field limit of {maximum} characters" + ) + + geometry_path = Path(geometry_file) + if not geometry_path.is_file(): + raise FileNotFoundError(f"Geometry file not found: {geometry_path}") + with open( + geometry_path, + "r", + encoding="utf-8", + errors="ignore", + newline="", + ) as geometry_stream: + geometry_lines = geometry_stream.readlines() + + geometry_records: List[Tuple[str, str]] = [] + pending_name: Optional[str] = None + for line in geometry_lines: + stripped = line.rstrip("\r\n") + if stripped.startswith("BC Line Name="): + pending_name = stripped[len("BC Line Name="):].strip() + elif stripped.startswith("BC Line Storage Area=") and pending_name is not None: + geometry_records.append( + ( + pending_name, + stripped[len("BC Line Storage Area="):].strip(), + ) + ) + pending_name = None + + exact_geometry_records = [ + record + for record in geometry_records + if record == (line_name, area_name) + ] + same_name_records = [ + record for record in geometry_records if record[0] == line_name + ] + if len(exact_geometry_records) != 1: + if same_name_records: + attached = sorted({record[1] for record in same_name_records}) + raise ValueError( + f"Geometry BC line {line_name!r} is attached to {attached}, " + f"not exactly once to {area_name!r}" + ) + raise ValueError( + f"Geometry file {geometry_path.name} does not contain BC line " + f"{line_name!r} on 2D Flow Area {area_name!r}" + ) + + unsteady_path = RasUnsteady._resolve_unsteady_file_path( + unsteady_file, + ras_object=ras_object, + ) + with open( + unsteady_path, + "r", + encoding="utf-8", + errors="ignore", + newline="", + ) as unsteady_stream: + lines = unsteady_stream.readlines() + newline = RasUnsteady._detect_line_ending(lines) + + boundary_locations: List[Tuple[int, str, List[str]]] = [] + for index, line in enumerate(lines): + if not line.startswith("Boundary Location="): + continue + location = line[len("Boundary Location="):].rstrip("\r\n") + boundary_locations.append( + (index, location, [part.strip() for part in location.split(",")]) + ) + + same_line_locations = [ + item + for item in boundary_locations + if len(item[2]) >= 8 and item[2][7] == line_name + ] + exact_locations = [ + item + for item in same_line_locations + if item[2][5] == area_name + ] + if len(exact_locations) > 1: + raise ValueError( + f"Unsteady file {unsteady_path.name} contains duplicate locations " + f"for {area_name!r}/{line_name!r}" + ) + if same_line_locations and not exact_locations: + attached = sorted( + { + item[2][5] + for item in same_line_locations + if len(item[2]) >= 6 + } + ) + raise ValueError( + f"Unsteady BC line {line_name!r} is already attached to {attached}, " + f"not {area_name!r}" + ) + + canonical_fields = ( + ("", 16), + ("", 16), + ("", 8), + ("", 8), + ("", 16), + (area_name, 16), + ("", 16), + (line_name, 32), + ) + canonical_location = ",".join( + f"{value:<{width}}" for value, width in canonical_fields + ) + boundary_count_before = len(boundary_locations) + + if exact_locations: + result = { + "unsteady_file": str(unsteady_path.resolve()), + "geometry_file": str(geometry_path.resolve()), + "area_2d": area_name, + "bc_line": line_name, + "created": False, + "insert_index": exact_locations[0][0], + "location": exact_locations[0][1], + "geometry_match_count": 1, + "boundary_count_before": boundary_count_before, + "boundary_count_after": boundary_count_before, + "boundaries_df_refreshed": False, + } + logger.info( + "2D boundary location already exists in %s: %s/%s", + unsteady_path.name, + area_name, + line_name, + ) + return result + + if boundary_locations: + insert_index = boundary_locations[0][0] + else: + header_prefixes = ("Flow Title=", "Program Version=", "Use Restart=") + header_indexes = [ + index + for index, line in enumerate(lines) + if line.startswith(header_prefixes) + ] + if not header_indexes: + raise ValueError( + f"Could not identify a safe boundary insertion point in {unsteady_path.name}" + ) + insert_index = max(header_indexes) + 1 + + if insert_index > 0 and not lines[insert_index - 1].endswith(("\r\n", "\n", "\r")): + lines[insert_index - 1] += newline + location_line = f"Boundary Location={canonical_location}{newline}" + lines.insert(insert_index, location_line) + RasUnsteady._atomic_write_lines(unsteady_path, lines) + + boundaries_df_refreshed = False + ras_obj = ras_object or ras + if ras_obj is not None: + try: + ras_obj.boundaries_df = ras_obj.get_boundary_conditions() + boundaries_df_refreshed = True + except Exception as exc: + logger.debug("boundaries_df refresh skipped: %s", exc) + + logger.info( + "Created 2D boundary location in %s: %s/%s", + unsteady_path.name, + area_name, + line_name, + ) + return { + "unsteady_file": str(unsteady_path.resolve()), + "geometry_file": str(geometry_path.resolve()), + "area_2d": area_name, + "bc_line": line_name, + "created": True, + "insert_index": insert_index, + "location": canonical_location, + "geometry_match_count": 1, + "boundary_count_before": boundary_count_before, + "boundary_count_after": boundary_count_before + 1, + "boundaries_df_refreshed": boundaries_df_refreshed, + } + @staticmethod @log_call def set_boundary_inline_hydrograph( @@ -7321,7 +7579,9 @@ def set_boundary_inline_hydrograph( river: Optional[str] = None, reach: Optional[str] = None, station: Optional[str] = None, - ras_object: Optional[Any] = None + ras_object: Optional[Any] = None, + area_2d: Optional[str] = None, + bc_line: Optional[str] = None, ) -> bool: """ Write an inline hydrograph table to a boundary condition, converting from DSS if needed. @@ -7359,6 +7619,12 @@ def set_boundary_inline_hydrograph( River station to locate the boundary. ras_object : optional Custom RAS object to use instead of the global one + area_2d : str, optional + Exact 2D Flow Area name (field 5 of ``Boundary Location=``). + Supply together with ``bc_line`` instead of the 1D selector. + bc_line : str, optional + Exact 2D BC line name (field 7 of ``Boundary Location=``). + Supply together with ``area_2d``. Returns ------- @@ -7408,6 +7674,8 @@ def set_boundary_inline_hydrograph( See Also -------- + ensure_2d_boundary_location : Create a geometry-validated empty 2D + boundary block before assigning its first boundary type. set_boundary_dss_link : Convert inline boundary to DSS mode set_precipitation_hyetograph : Write an incremental-depth precipitation hyetograph """ @@ -7427,11 +7695,22 @@ def set_boundary_inline_hydrograph( table_keyword = SUPPORTED_TYPES[bc_type] + has_1d_selector = any(value is not None for value in (river, reach, station)) + has_2d_selector = any(value is not None for value in (area_2d, bc_line)) + if has_1d_selector and has_2d_selector: + raise ValueError( + "Provide either (river, reach, station) or (area_2d, bc_line), not both" + ) + if has_1d_selector and not (river and reach and station): + raise ValueError("1D selector requires river, reach, and station") + if has_2d_selector and not (area_2d and bc_line): + raise ValueError("2D selector requires area_2d and bc_line") + ras_obj = ras_object or ras if ras_obj is not None: try: ras_obj.check_initialized() - except: + except Exception: pass # Resolve unsteady file path @@ -7456,52 +7735,83 @@ def set_boundary_inline_hydrograph( ) # Extract values - hours = hydrograph_df['hour'].values - values = hydrograph_df['value'].values + try: + hours = np.asarray(hydrograph_df['hour'].values, dtype=float) + values = np.asarray(hydrograph_df['value'].values, dtype=float) + except (TypeError, ValueError) as exc: + raise ValueError("Hydrograph hour and value columns must be numeric") from exc num_values = len(values) if num_values < 2: raise ValueError("DataFrame must have at least 2 rows") + if not np.isfinite(hours).all() or not np.isfinite(values).all(): + raise ValueError("Hydrograph hours and values must be finite") # Calculate interval from hour column - interval_hours = float(hours[1] - hours[0]) - if interval_hours >= 1.0: - if interval_hours == int(interval_hours): - interval_str = f"{int(interval_hours)}HOUR" - else: - interval_min = int(interval_hours * 60) - interval_str = f"{interval_min}MIN" + intervals = np.diff(hours) + interval_hours = float(intervals[0]) + if interval_hours <= 0 or not np.allclose( + intervals, + interval_hours, + rtol=0.0, + atol=1e-10, + ): + raise ValueError("Hydrograph hours must be strictly increasing and evenly spaced") + interval_minutes = interval_hours * 60.0 + rounded_minutes = round(interval_minutes) + if rounded_minutes < 1 or not math.isclose( + interval_minutes, + rounded_minutes, + rel_tol=0.0, + abs_tol=1e-8, + ): + raise ValueError("Hydrograph interval must be a positive whole number of minutes") + if rounded_minutes % 60 == 0: + interval_str = f"{rounded_minutes // 60}HOUR" else: - interval_min = int(interval_hours * 60) - interval_str = f"{interval_min}MIN" + interval_str = f"{rounded_minutes}MIN" # Format values as fixed-width (8 chars each, 10 values per line) # Flow/Stage hydrographs use values only (not time-value pairs) - formatted_lines = [] + formatted_rows = [] for i in range(0, num_values, 10): row_values = values[i:i+10] formatted_row = ''.join(f'{v:8.2f}' if abs(v) < 1e7 else f'{v:8.1f}' for v in row_values) - formatted_lines.append(formatted_row + '\n') + formatted_rows.append(formatted_row) # Read the file - with open(unsteady_path, 'r', encoding='utf-8', errors='ignore') as f: + with open( + unsteady_path, + 'r', + encoding='utf-8', + errors='ignore', + newline='', + ) as f: lines = f.readlines() + newline = RasUnsteady._detect_line_ending(lines) + formatted_lines = [row + newline for row in formatted_rows] # Find the target boundary target_boundary_idx = None for idx, line in enumerate(lines): if line.startswith('Boundary Location='): - if river is not None or reach is not None or station is not None: - loc_line = line.replace('Boundary Location=', '') + if has_1d_selector or has_2d_selector: + loc_line = line[len('Boundary Location='):] parts = [p.strip() for p in loc_line.split(',')] - match = True - if river is not None and (len(parts) < 1 or parts[0] != river): - match = False - if reach is not None and (len(parts) < 2 or parts[1] != reach): - match = False - if station is not None and (len(parts) < 3 or parts[2] != station): - match = False + if has_1d_selector: + match = ( + len(parts) >= 3 + and parts[0] == river + and parts[1] == reach + and parts[2] == station + ) + else: + match = ( + len(parts) >= 8 + and parts[5] == area_2d + and parts[7] == bc_line + ) if match: target_boundary_idx = idx break @@ -7518,7 +7828,12 @@ def set_boundary_inline_hydrograph( break if target_boundary_idx is None: - loc_str = f"{river}/{reach}/{station}" if river else "first matching" + if has_1d_selector: + loc_str = f"{river}/{reach}/{station}" + elif has_2d_selector: + loc_str = f"{area_2d}/{bc_line}" + else: + loc_str = "first matching" logger.warning(f"Boundary not found for {bc_type}: {loc_str}") return False @@ -7590,11 +7905,11 @@ def set_boundary_inline_hydrograph( # Step 2: Update table header with new count if table_header_idx is not None: - lines[table_header_idx] = f'{table_keyword} {num_values} \n' + lines[table_header_idx] = f'{table_keyword} {num_values} {newline}' else: # Insert table header after Interval line (or after boundary location) insert_pos = interval_idx + 1 if interval_idx is not None else target_boundary_idx + 1 - lines.insert(insert_pos, f'{table_keyword} {num_values} \n') + lines.insert(insert_pos, f'{table_keyword} {num_values} {newline}') table_header_idx = insert_pos # Adjust indices after insertion if interval_idx is not None and interval_idx >= insert_pos: @@ -7624,10 +7939,10 @@ def set_boundary_inline_hydrograph( # Step 4: Update Interval line if interval_idx is not None: - lines[interval_idx] = f'Interval={interval_str}\n' + lines[interval_idx] = f'Interval={interval_str}{newline}' else: # Insert interval before table header - lines.insert(table_header_idx, f'Interval={interval_str}\n') + lines.insert(table_header_idx, f'Interval={interval_str}{newline}') # Adjust all indices after this insertion if dss_file_idx is not None and dss_file_idx >= table_header_idx: dss_file_idx += 1 @@ -7638,15 +7953,30 @@ def set_boundary_inline_hydrograph( # Step 5: Set Use DSS=False and clear DSS File/Path if use_dss_idx is not None: - lines[use_dss_idx] = 'Use DSS=False\n' + lines[use_dss_idx] = f'Use DSS=False{newline}' + else: + block_end = next( + ( + index + for index in range(target_boundary_idx + 1, len(lines)) + if lines[index].startswith('Boundary Location=') + ), + len(lines), + ) + lines.insert(block_end, f'Use DSS=False{newline}') if dss_file_idx is not None: - lines[dss_file_idx] = 'DSS File=\n' + lines[dss_file_idx] = f'DSS File={newline}' if dss_path_idx is not None: - lines[dss_path_idx] = 'DSS Path=\n' + lines[dss_path_idx] = f'DSS Path={newline}' # Write updated content back to file - with open(unsteady_path, 'w', encoding='utf-8') as f: - f.writelines(lines) + RasUnsteady._atomic_write_lines(unsteady_path, lines) + + if ras_obj is not None: + try: + ras_obj.boundaries_df = ras_obj.get_boundary_conditions() + except Exception as exc: + logger.debug("boundaries_df refresh skipped: %s", exc) peak_value = float(np.max(values)) logger.info( diff --git a/tests/test_rasunsteady_2d_boundary_location.py b/tests/test_rasunsteady_2d_boundary_location.py new file mode 100644 index 000000000..41b1d68d7 --- /dev/null +++ b/tests/test_rasunsteady_2d_boundary_location.py @@ -0,0 +1,298 @@ +from __future__ import annotations + +import shutil +from pathlib import Path + +import pandas as pd +import pytest + +from ras_commander import RasUnsteady + + +def _write_geometry( + path: Path, + *, + area_name: str = "Breakout Area", + line_name: str = "Breakout Inflow", + duplicate: bool = False, +) -> Path: + block = ( + f"BC Line Name={line_name:<40}\r\n" + f"BC Line Storage Area={area_name:<16}\r\n" + "BC Line Start Position= 0 , 0 \r\n" + "BC Line End Position= 0 , 10 \r\n" + "BC Line Arc= 2 \r\n" + " 0 0 0 10\r\n" + "BC Line Text Position= 1.79769313486232E+308 , 1.79769313486232E+308 \r\n" + ) + path.write_bytes(("Geom Title=breakout\r\n" + block * (2 if duplicate else 1)).encode()) + return path + + +def _existing_boundary(area_name: str, line_name: str) -> str: + fields = ( + ("", 16), + ("", 16), + ("", 8), + ("", 8), + ("", 16), + (area_name, 16), + ("", 16), + (line_name, 32), + ) + return "Boundary Location=" + ",".join( + f"{value:<{width}}" for value, width in fields + ) + + +def _write_unsteady(path: Path) -> Path: + text = ( + "Flow Title=breakout\r\n" + "Program Version=6.60\r\n" + "Use Restart= 0 \r\n" + f"{_existing_boundary('Parent Area', 'Parent Outflow')}\r\n" + "Friction Slope=0.0003\r\n" + "Precipitation Mode=Disable\r\n" + ) + path.write_bytes(text.encode()) + return path + + +def _boundary_block(text: str, area_name: str, line_name: str) -> str: + marker = _existing_boundary(area_name, line_name).rstrip() + start = text.index(marker) + next_boundary = text.find("Boundary Location=", start + len(marker)) + return text[start:] if next_boundary < 0 else text[start:next_boundary] + + +def test_ensure_location_then_author_inline_flow_hydrograph(tmp_path): + geometry = _write_geometry(tmp_path / "breakout.g02") + unsteady = _write_unsteady(tmp_path / "breakout.u02") + + created = RasUnsteady.ensure_2d_boundary_location( + unsteady, + geometry, + area_2d="Breakout Area", + bc_line="Breakout Inflow", + ) + + assert created["created"] is True + assert created["geometry_match_count"] == 1 + assert created["boundary_count_before"] == 1 + assert created["boundary_count_after"] == 2 + + hydrograph = pd.DataFrame( + {"hour": [0.0, 0.5, 1.0, 1.5, 2.0], "value": [10, 25, 50, 25, 10]} + ) + assert RasUnsteady.set_boundary_inline_hydrograph( + unsteady, + hydrograph, + bc_type="Flow Hydrograph", + area_2d="Breakout Area", + bc_line="Breakout Inflow", + ) + slope = RasUnsteady.set_flow_hydrograph_slope( + unsteady, + 0.0005, + area_2d="Breakout Area", + bc_line="Breakout Inflow", + ) + assert slope["new_eg_slope"] == 0.0005 + + raw = unsteady.read_bytes() + assert b"\n" not in raw.replace(b"\r\n", b"") + text = raw.decode() + block = _boundary_block(text, "Breakout Area", "Breakout Inflow") + assert "Interval=30MIN\r\n" in block + assert "Flow Hydrograph= 5 \r\n" in block + assert "Flow Hydrograph Slope= 0.0005 " in block + assert "Use DSS=False\r\n" in block + assert block.index(" 10.00") < block.index("Use DSS=False") + + before = unsteady.read_bytes() + existing = RasUnsteady.ensure_2d_boundary_location( + unsteady, + geometry, + area_2d="Breakout Area", + bc_line="Breakout Inflow", + ) + assert existing["created"] is False + assert unsteady.read_bytes() == before + + +def test_ensure_location_supports_empty_unsteady_boundary_collection(tmp_path): + geometry = _write_geometry(tmp_path / "breakout.g02") + unsteady = tmp_path / "breakout.u02" + unsteady.write_bytes( + b"Flow Title=breakout\r\nProgram Version=6.60\r\nUse Restart= 0 \r\n" + b"Precipitation Mode=Disable\r\n" + ) + + result = RasUnsteady.ensure_2d_boundary_location( + unsteady, + geometry, + area_2d="Breakout Area", + bc_line="Breakout Inflow", + ) + + lines = unsteady.read_text(encoding="utf-8").splitlines() + assert result["insert_index"] == 3 + assert lines[3].startswith("Boundary Location=") + assert lines[4] == "Precipitation Mode=Disable" + + +def test_ensure_location_then_author_normal_depth_outflow(tmp_path): + geometry = _write_geometry( + tmp_path / "breakout.g02", + line_name="Breakout Outflow", + ) + unsteady = _write_unsteady(tmp_path / "breakout.u02") + RasUnsteady.ensure_2d_boundary_location( + unsteady, + geometry, + area_2d="Breakout Area", + bc_line="Breakout Outflow", + ) + + result = RasUnsteady.set_normal_depth_boundary( + unsteady, + friction_slope=0.0004, + area_2d="Breakout Area", + bc_line="Breakout Outflow", + ) + + assert result["previous_bc_type"] is None + assert result["new_friction_slope"] == 0.0004 + assert result["lines_inserted"] == 1 + block = _boundary_block( + unsteady.read_text(encoding="utf-8"), + "Breakout Area", + "Breakout Outflow", + ) + assert "Friction Slope=0.0004\n" in block + + +@pytest.mark.parametrize( + ("geometry_area", "geometry_line", "requested_area", "requested_line", "match"), + [ + ("Other Area", "Breakout Inflow", "Breakout Area", "Breakout Inflow", "attached"), + ("Breakout Area", "Other Inflow", "Breakout Area", "Breakout Inflow", "does not contain"), + ], +) +def test_ensure_location_rejects_geometry_mismatch( + tmp_path, + geometry_area, + geometry_line, + requested_area, + requested_line, + match, +): + geometry = _write_geometry( + tmp_path / "breakout.g02", + area_name=geometry_area, + line_name=geometry_line, + ) + unsteady = _write_unsteady(tmp_path / "breakout.u02") + + with pytest.raises(ValueError, match=match): + RasUnsteady.ensure_2d_boundary_location( + unsteady, + geometry, + area_2d=requested_area, + bc_line=requested_line, + ) + + +def test_ensure_location_rejects_duplicate_geometry_records(tmp_path): + geometry = _write_geometry(tmp_path / "breakout.g02", duplicate=True) + unsteady = _write_unsteady(tmp_path / "breakout.u02") + + with pytest.raises(ValueError, match="not exactly once"): + RasUnsteady.ensure_2d_boundary_location( + unsteady, + geometry, + area_2d="Breakout Area", + bc_line="Breakout Inflow", + ) + + +def test_ensure_location_rejects_existing_name_on_other_area(tmp_path): + geometry = _write_geometry(tmp_path / "breakout.g02") + unsteady = _write_unsteady(tmp_path / "breakout.u02") + text = unsteady.read_text(encoding="utf-8").replace( + "Parent Outflow", + "Breakout Inflow", + ) + unsteady.write_text(text, encoding="utf-8") + + with pytest.raises(ValueError, match="already attached"): + RasUnsteady.ensure_2d_boundary_location( + unsteady, + geometry, + area_2d="Breakout Area", + bc_line="Breakout Inflow", + ) + + +@pytest.mark.parametrize( + "hours", + [ + [0.0, 1.0, 0.5], + [0.0, 0.5, 1.25], + [0.0, 0.001, 0.002], + ], +) +def test_inline_hydrograph_rejects_invalid_time_axis(tmp_path, hours): + geometry = _write_geometry(tmp_path / "breakout.g02") + unsteady = _write_unsteady(tmp_path / "breakout.u02") + RasUnsteady.ensure_2d_boundary_location( + unsteady, + geometry, + area_2d="Breakout Area", + bc_line="Breakout Inflow", + ) + + with pytest.raises(ValueError, match="Hydrograph"): + RasUnsteady.set_boundary_inline_hydrograph( + unsteady, + pd.DataFrame({"hour": hours, "value": [1.0, 2.0, 3.0]}), + area_2d="Breakout Area", + bc_line="Breakout Inflow", + ) + + +def test_inline_hydrograph_rejects_mixed_selector_groups(tmp_path): + unsteady = _write_unsteady(tmp_path / "breakout.u02") + + with pytest.raises(ValueError, match="not both"): + RasUnsteady.set_boundary_inline_hydrograph( + unsteady, + pd.DataFrame({"hour": [0.0, 1.0], "value": [1.0, 2.0]}), + river="River", + reach="Reach", + station="1", + area_2d="Breakout Area", + bc_line="Breakout Inflow", + ) + + +def test_real_bald_eagle_existing_2d_location_round_trips_read_only(tmp_path): + fixture = Path(__file__).parents[1] / "example_projects" / "BaldEagleCrkMulti2D" + source_geometry = fixture / "BaldEagleDamBrk.g09" + source_unsteady = fixture / "BaldEagleDamBrk.u03" + if not source_geometry.is_file() or not source_unsteady.is_file(): + pytest.skip("BaldEagleCrkMulti2D g09/u03 fixture is unavailable") + geometry = Path(shutil.copy2(source_geometry, tmp_path / source_geometry.name)) + unsteady = Path(shutil.copy2(source_unsteady, tmp_path / source_unsteady.name)) + before = unsteady.read_bytes() + + result = RasUnsteady.ensure_2d_boundary_location( + unsteady, + geometry, + area_2d="BaldEagleCr", + bc_line="Upstream Inflow", + ) + + assert result["created"] is False + assert result["geometry_match_count"] == 1 + assert unsteady.read_bytes() == before From 6b855c254514ec0ac73a63410a24befeff636ee4 Mon Sep 17 00:00:00 2001 From: "William Mark Katzenmeyer, P.E., C.F.M." Date: Mon, 31 Aug 2026 00:03:42 -0400 Subject: [PATCH 3/5] Add atomic 2D boundary collection replacement --- ras_commander/RasUnsteady.py | 228 ++++++++++++++++++ .../test_rasunsteady_2d_boundary_location.py | 77 ++++++ 2 files changed, 305 insertions(+) diff --git a/ras_commander/RasUnsteady.py b/ras_commander/RasUnsteady.py index eaad3a472..cecf760b3 100644 --- a/ras_commander/RasUnsteady.py +++ b/ras_commander/RasUnsteady.py @@ -68,6 +68,7 @@ def my_function(): - inspect_boundary_blocks() - Inventory exact Boundary Location blocks in an owned stage - delete_boundary() - Preview/remove one exact block from an owned staged project - ensure_2d_boundary_location() - Create an empty 2D boundary block after validating geometry +- replace_2d_boundary_locations() - Reset 2D boundary blocks from geometry-backed specs - update_dss_run_identifier() - Update DSS path F-part for new scenarios - set_boundary_dss_link() - Convert inline BC to DSS-linked (complete state transition) - set_boundary_inline_hydrograph() - Write inline hydrograph, convert DSS to inline @@ -7570,6 +7571,233 @@ def ensure_2d_boundary_location( "boundaries_df_refreshed": boundaries_df_refreshed, } + @staticmethod + @log_call + def replace_2d_boundary_locations( + unsteady_file: Union[str, Path], + geometry_file: Union[str, Path], + locations: List[Dict[str, str]], + *, + ras_object: Optional[Any] = None, + ) -> Dict[str, Any]: + """Replace every 2D boundary block with validated empty locations. + + Existing 2D Flow Area boundary blocks—including area-wide blocks with + no BC-line name—are removed with their complete type/data payloads. + River/reach and storage-area/structure boundary blocks are preserved. + Every replacement must identify exactly one BC line in the supplied + plain-text geometry. This is intentionally a reset operation: callers + assign new Flow Hydrograph, Normal Depth, or other types afterward. + + Parameters + ---------- + unsteady_file : str or Path + Unsteady-flow file path or number resolvable through + ``ras_object``. + geometry_file : str or Path + Exact plain-text geometry file containing all replacement lines. + locations : list of dict + Non-empty list with exact ``area_2d`` and ``bc_line`` strings. + BC-line names must be unique. + ras_object : optional + Project object used for short-number resolution and DataFrame + refresh. + + Returns + ------- + dict + Removed 2D locations, inserted locations, preserved non-2D block + count, insertion index, and DataFrame-refresh evidence. + """ + if not isinstance(locations, list) or not locations: + raise ValueError("locations must be a non-empty list of dicts") + + desired: List[Tuple[str, str]] = [] + for index, spec in enumerate(locations): + if not isinstance(spec, dict): + raise ValueError(f"locations[{index}] must be a dict") + area_name = str(spec.get("area_2d", "")).strip() + line_name = str(spec.get("bc_line", "")).strip() + for value, label, maximum in ( + (area_name, f"locations[{index}].area_2d", 16), + (line_name, f"locations[{index}].bc_line", 32), + ): + if not value: + raise ValueError(f"{label} must be non-empty") + if any(character in value for character in (",", "\r", "\n")): + raise ValueError(f"{label} cannot contain commas or newlines") + if len(value) > maximum: + raise ValueError( + f"{label} exceeds the HEC-RAS fixed-field limit of " + f"{maximum} characters" + ) + desired.append((area_name, line_name)) + if len({line_name for _, line_name in desired}) != len(desired): + raise ValueError("Replacement 2D BC-line names must be unique") + + geometry_path = Path(geometry_file) + if not geometry_path.is_file(): + raise FileNotFoundError(f"Geometry file not found: {geometry_path}") + with open( + geometry_path, + "r", + encoding="utf-8", + errors="ignore", + newline="", + ) as geometry_stream: + geometry_lines = geometry_stream.readlines() + geometry_records: List[Tuple[str, str]] = [] + pending_name: Optional[str] = None + for line in geometry_lines: + stripped = line.rstrip("\r\n") + if stripped.startswith("BC Line Name="): + pending_name = stripped[len("BC Line Name="):].strip() + elif stripped.startswith("BC Line Storage Area=") and pending_name is not None: + geometry_records.append( + ( + stripped[len("BC Line Storage Area="):].strip(), + pending_name, + ) + ) + pending_name = None + for area_name, line_name in desired: + match_count = geometry_records.count((area_name, line_name)) + if match_count != 1: + raise ValueError( + f"Geometry must contain exactly one BC line {line_name!r} " + f"on 2D Flow Area {area_name!r}; found {match_count}" + ) + + unsteady_path = RasUnsteady._resolve_unsteady_file_path( + unsteady_file, + ras_object=ras_object, + ) + with open( + unsteady_path, + "r", + encoding="utf-8", + errors="ignore", + newline="", + ) as unsteady_stream: + lines = unsteady_stream.readlines() + newline = RasUnsteady._detect_line_ending(lines) + + starts = [ + index + for index, line in enumerate(lines) + if line.startswith("Boundary Location=") + ] + removed_locations: List[Dict[str, str]] = [] + removed_ranges: List[Tuple[int, int]] = [] + preserved_block_count = 0 + for position, start in enumerate(starts): + end = starts[position + 1] if position + 1 < len(starts) else len(lines) + location = lines[start][len("Boundary Location="):].rstrip("\r\n") + fields = [part.strip() for part in location.split(",")] + is_2d = len(fields) >= 6 and bool(fields[5]) + malformed_2d = len(fields) >= 8 and bool(fields[7]) and not is_2d + if malformed_2d: + raise ValueError( + f"Malformed 2D boundary location in {unsteady_path.name}: {location!r}" + ) + if is_2d: + removed_ranges.append((start, end)) + removed_locations.append( + { + "area_2d": fields[5], + "bc_line": fields[7] if len(fields) >= 8 else "", + } + ) + else: + preserved_block_count += 1 + + removed_indexes = { + index + for start, end in removed_ranges + for index in range(start, end) + } + filtered_lines = [ + line for index, line in enumerate(lines) if index not in removed_indexes + ] + + remaining_boundary_index = next( + ( + index + for index, line in enumerate(filtered_lines) + if line.startswith("Boundary Location=") + ), + None, + ) + if remaining_boundary_index is not None: + insert_index = remaining_boundary_index + else: + header_prefixes = ("Flow Title=", "Program Version=", "Use Restart=") + header_indexes = [ + index + for index, line in enumerate(filtered_lines) + if line.startswith(header_prefixes) + ] + if not header_indexes: + raise ValueError( + f"Could not identify a safe boundary insertion point in {unsteady_path.name}" + ) + insert_index = max(header_indexes) + 1 + + if ( + insert_index > 0 + and not filtered_lines[insert_index - 1].endswith(("\r\n", "\n", "\r")) + ): + filtered_lines[insert_index - 1] += newline + replacement_lines = [] + for area_name, line_name in desired: + fields = ( + ("", 16), + ("", 16), + ("", 8), + ("", 8), + ("", 16), + (area_name, 16), + ("", 16), + (line_name, 32), + ) + location = ",".join(f"{value:<{width}}" for value, width in fields) + replacement_lines.append(f"Boundary Location={location}{newline}") + updated_lines = ( + filtered_lines[:insert_index] + + replacement_lines + + filtered_lines[insert_index:] + ) + RasUnsteady._atomic_write_lines(unsteady_path, updated_lines) + + boundaries_df_refreshed = False + ras_obj = ras_object or ras + if ras_obj is not None: + try: + ras_obj.boundaries_df = ras_obj.get_boundary_conditions() + boundaries_df_refreshed = True + except Exception as exc: + logger.debug("boundaries_df refresh skipped: %s", exc) + + inserted_locations = [ + {"area_2d": area_name, "bc_line": line_name} + for area_name, line_name in desired + ] + logger.info( + "Replaced %d 2D boundary block(s) with %d geometry-backed location(s) in %s", + len(removed_locations), + len(inserted_locations), + unsteady_path.name, + ) + return { + "unsteady_file": str(unsteady_path.resolve()), + "geometry_file": str(geometry_path.resolve()), + "removed_locations": removed_locations, + "inserted_locations": inserted_locations, + "preserved_non_2d_block_count": preserved_block_count, + "insert_index": insert_index, + "boundaries_df_refreshed": boundaries_df_refreshed, + } + @staticmethod @log_call def set_boundary_inline_hydrograph( diff --git a/tests/test_rasunsteady_2d_boundary_location.py b/tests/test_rasunsteady_2d_boundary_location.py index 41b1d68d7..1b195ebf4 100644 --- a/tests/test_rasunsteady_2d_boundary_location.py +++ b/tests/test_rasunsteady_2d_boundary_location.py @@ -29,6 +29,22 @@ def _write_geometry( return path +def _write_geometry_locations(path: Path, locations: list[tuple[str, str]]) -> Path: + blocks = [] + for area_name, line_name in locations: + blocks.append( + f"BC Line Name={line_name:<40}\r\n" + f"BC Line Storage Area={area_name:<16}\r\n" + "BC Line Start Position= 0 , 0 \r\n" + "BC Line End Position= 0 , 10 \r\n" + "BC Line Arc= 2 \r\n" + " 0 0 0 10\r\n" + "BC Line Text Position= 1.79769313486232E+308 , 1.79769313486232E+308 \r\n" + ) + path.write_bytes(("Geom Title=breakout\r\n" + "".join(blocks)).encode()) + return path + + def _existing_boundary(area_name: str, line_name: str) -> str: fields = ( ("", 16), @@ -172,6 +188,67 @@ def test_ensure_location_then_author_normal_depth_outflow(tmp_path): assert "Friction Slope=0.0004\n" in block +def test_replace_2d_locations_removes_complete_2d_blocks_and_preserves_non_2d(tmp_path): + geometry = _write_geometry_locations( + tmp_path / "breakout.g02", + [ + ("Breakout Area", "Breakout Inflow"), + ("Breakout Area", "Breakout Outflow"), + ], + ) + unsteady = tmp_path / "breakout.u02" + area_wide_fields = ["", "", "", "", "", "Parent Area", "", ""] + one_d_fields = ["River", "Reach", "100", "", "", "", "", ""] + gate_fields = ["", "", "", "", "Storage Area", "", "", ""] + unsteady.write_bytes( + ( + "Flow Title=breakout\r\n" + "Program Version=6.60\r\n" + "Use Restart= 0 \r\n" + f"{_existing_boundary('Parent Area', 'Parent Inflow')}\r\n" + "Interval=1HOUR\r\n" + "Flow Hydrograph= 2 \r\n" + " 1.00 2.00\r\n" + "Use DSS=False\r\n" + f"Boundary Location={','.join(area_wide_fields)}\r\n" + "Precipitation Hydrograph= 2 \r\n" + " 0.10 0.20\r\n" + f"Boundary Location={','.join(one_d_fields)}\r\n" + "Friction Slope=0.0003,0\r\n" + f"Boundary Location={','.join(gate_fields)}\r\n" + "Gate Name=Gate 1\r\n" + "Precipitation Mode=Enable\r\n" + ).encode() + ) + + result = RasUnsteady.replace_2d_boundary_locations( + unsteady, + geometry, + [ + {"area_2d": "Breakout Area", "bc_line": "Breakout Inflow"}, + {"area_2d": "Breakout Area", "bc_line": "Breakout Outflow"}, + ], + ) + + assert result["removed_locations"] == [ + {"area_2d": "Parent Area", "bc_line": "Parent Inflow"}, + {"area_2d": "Parent Area", "bc_line": ""}, + ] + assert result["preserved_non_2d_block_count"] == 2 + assert result["inserted_locations"] == [ + {"area_2d": "Breakout Area", "bc_line": "Breakout Inflow"}, + {"area_2d": "Breakout Area", "bc_line": "Breakout Outflow"}, + ] + text = unsteady.read_text(encoding="utf-8") + assert "Parent Inflow" not in text + assert "Precipitation Hydrograph=" not in text + assert " 1.00 2.00" not in text + assert "River,Reach,100" in text + assert "Gate Name=Gate 1" in text + assert text.count("Breakout Inflow") == 1 + assert text.count("Breakout Outflow") == 1 + + @pytest.mark.parametrize( ("geometry_area", "geometry_line", "requested_area", "requested_line", "match"), [ From 031a57cf69c83b01fe64b0dc2748bc68995324f9 Mon Sep 17 00:00:00 2001 From: "William Mark Katzenmeyer, P.E., C.F.M." Date: Mon, 31 Aug 2026 00:23:48 -0400 Subject: [PATCH 4/5] Require zero-based inline hydrograph timing --- ras_commander/RasUnsteady.py | 2 ++ tests/test_rasunsteady_2d_boundary_location.py | 1 + 2 files changed, 3 insertions(+) diff --git a/ras_commander/RasUnsteady.py b/ras_commander/RasUnsteady.py index cecf760b3..4b7535676 100644 --- a/ras_commander/RasUnsteady.py +++ b/ras_commander/RasUnsteady.py @@ -7974,6 +7974,8 @@ def set_boundary_inline_hydrograph( raise ValueError("DataFrame must have at least 2 rows") if not np.isfinite(hours).all() or not np.isfinite(values).all(): raise ValueError("Hydrograph hours and values must be finite") + if not math.isclose(float(hours[0]), 0.0, rel_tol=0.0, abs_tol=1e-12): + raise ValueError("Hydrograph hours must start at zero") # Calculate interval from hour column intervals = np.diff(hours) diff --git a/tests/test_rasunsteady_2d_boundary_location.py b/tests/test_rasunsteady_2d_boundary_location.py index 1b195ebf4..bde7f874d 100644 --- a/tests/test_rasunsteady_2d_boundary_location.py +++ b/tests/test_rasunsteady_2d_boundary_location.py @@ -314,6 +314,7 @@ def test_ensure_location_rejects_existing_name_on_other_area(tmp_path): @pytest.mark.parametrize( "hours", [ + [1.0, 1.5, 2.0], [0.0, 1.0, 0.5], [0.0, 0.5, 1.25], [0.0, 0.001, 0.002], From f32769816a1dc64e6d134e6d25f6b428cf57b34a Mon Sep 17 00:00:00 2001 From: "William Mark Katzenmeyer, P.E., C.F.M." Date: Mon, 31 Aug 2026 07:32:32 -0400 Subject: [PATCH 5/5] Harden exact 2D Linux execution preparation --- ras_commander/RasCmdr.py | 145 +++++++- ras_commander/RasMap.py | 172 +++++++++ ras_commander/RasProcess.py | 329 +++++++++++++++++- ras_commander/RasUnsteady.py | 221 ++++++++++++ ras_commander/_gdal_runtime.py | 19 +- ras_commander/geom/GeomPreprocessor.py | 138 +++++++- tests/test_geom_preprocessor.py | 66 ++++ tests/test_linux_robustness.py | 109 +++++- tests/test_ras_geometry_compute_unit.py | 98 ++++++ tests/test_rasmap_event_conditions.py | 86 +++++ tests/test_rasunsteady_disable_meteorology.py | 156 +++++++++ 11 files changed, 1505 insertions(+), 34 deletions(-) create mode 100644 tests/test_rasmap_event_conditions.py create mode 100644 tests/test_rasunsteady_disable_meteorology.py diff --git a/ras_commander/RasCmdr.py b/ras_commander/RasCmdr.py index 35eff7465..8188457e5 100644 --- a/ras_commander/RasCmdr.py +++ b/ras_commander/RasCmdr.py @@ -2145,13 +2145,40 @@ def compute_plan_linux( f"See examples/510_linux_execution.ipynb for the complete workflow." ) - # Set num_cores if specified + effective_num_cores = None + # Set num_cores in both the text plan and the compiled execution HDF. + # Native RasUnsteady reads the latter; changing only .p## leaves the + # Phase-1 core count in force. if num_cores is not None: try: - RasPlan.set_num_cores(plan_path, num_cores=num_cores, ras_object=ras_obj) - logger.info(f"Set number of cores to {num_cores} for plan: {plan_num_str}") + effective_num_cores = RasCmdr._effective_linux_core_count(num_cores) + RasPlan.set_num_cores( + plan_path, + num_cores=effective_num_cores, + ras_object=ras_obj, + ) + hdf_core_evidence = None + if not layout["needs_c_file"]: + hdf_core_evidence = RasCmdr._set_linux_hdf_num_cores( + tmp_hdf, + effective_num_cores, + ) + logger.info( + "Configured %d native solver core(s) for plan %s%s", + effective_num_cores, + plan_num_str, + ( + f" in {len(hdf_core_evidence['updated_attributes'])} " + "compiled-HDF attribute(s)" + if hdf_core_evidence is not None + else "" + ), + ) except Exception as e: - logger.error(f"Error setting number of cores: {e}") + raise RuntimeError( + f"Could not configure native solver cores for plan " + f"{plan_num_str}: {e}" + ) from e if run_via_wsl: return RasCmdr._compute_plan_linux_via_wsl( @@ -2224,6 +2251,9 @@ def _create_io_link(source: Path, io_name: str): env = os.environ.copy() env["LD_LIBRARY_PATH"] = ld_path + if effective_num_cores is not None: + env["OMP_NUM_THREADS"] = str(effective_num_cores) + env["MKL_NUM_THREADS"] = str(effective_num_cores) log_path = project_dir / f"compute_linux_{plan_num_str}.log" success = False @@ -2384,6 +2414,87 @@ def _resolve_linux_layout(ras_exe_dir: Path) -> dict: "label": "canonical", } + @staticmethod + def _effective_linux_core_count(requested_cores: int) -> int: + """Cap a requested solver core count to the process affinity envelope.""" + if ( + isinstance(requested_cores, bool) + or not isinstance(requested_cores, Number) + or int(requested_cores) != requested_cores + or int(requested_cores) < 1 + ): + raise ValueError("num_cores must be a positive integer") + requested = int(requested_cores) + + available = None + affinity_reader = getattr(os, "sched_getaffinity", None) + if callable(affinity_reader): + try: + available = len(affinity_reader(0)) + except (OSError, TypeError): + available = None + if not available: + cpu_reader = getattr(os, "cpu_count", None) + if callable(cpu_reader): + available = cpu_reader() + if not available: + return requested + + effective = min(requested, int(available)) + if effective < requested: + logger.warning( + "Capped requested native solver cores from %d to %d based on " + "the process affinity envelope", + requested, + effective, + ) + return effective + + @staticmethod + def _set_linux_hdf_num_cores(tmp_hdf: Path, num_cores: int) -> dict: + """Write the effective core count into a canonical plan ``*.tmp.hdf``.""" + import h5py + import numpy as np + + tmp_hdf = Path(tmp_hdf) + if not tmp_hdf.name.casefold().endswith(".tmp.hdf"): + raise ValueError("Core control target must be a '*.tmp.hdf' file") + if not tmp_hdf.is_file(): + raise FileNotFoundError(f"Compiled plan HDF not found: {tmp_hdf}") + + updated_attributes = [] + with h5py.File(tmp_hdf, "r+") as hdf_file: + parameters = hdf_file.get("Plan Data/Plan Parameters") + if parameters is None: + raise ValueError("Compiled plan HDF lacks /Plan Data/Plan Parameters") + for attribute_name in ("1D Cores", "2D Cores (per mesh)"): + if attribute_name not in parameters.attrs: + continue + prior = np.asarray(parameters.attrs[attribute_name]) + replacement = np.full(prior.shape, num_cores, dtype=prior.dtype) + if prior.shape == (): + replacement = replacement[()] + parameters.attrs.modify(attribute_name, replacement) + updated_attributes.append( + { + "attribute": attribute_name, + "before": prior.tolist(), + "after": np.asarray( + parameters.attrs[attribute_name] + ).tolist(), + } + ) + if not updated_attributes: + raise ValueError( + "Compiled plan HDF contains no supported solver-core attributes" + ) + hdf_file.flush() + return { + "path": str(tmp_hdf), + "effective_cores": num_cores, + "updated_attributes": updated_attributes, + } + @staticmethod def _build_linux_ld_path(ras_exe_dir: Path, layout: dict) -> str: """Build LD_LIBRARY_PATH for a Linux RasUnsteady run, per layout (CLB-886).""" @@ -2439,14 +2550,38 @@ def _validate_linux_solve(log_path, result_hdf, plan_num_str: str): for marker in error_markers: if marker in low: return False, f"solver log reports failure ('{marker}')" + import re + + explicit_error = re.search( + r"(?im)^\s*(?:error\s*:|hdf_error\b)", + log_text, + ) + if explicit_error: + return False, ( + "solver log reports failure " + f"('{explicit_error.group(0).strip()}')" + ) + if "finished unsteady flow simulation" not in low: + return False, "solver log missing 'Finished Unsteady Flow Simulation' banner" try: import h5py with h5py.File(str(result_hdf), "r") as hf: results = hf.get("Results") if results is None: return False, "result HDF missing /Results group" - if results.get("Unsteady") is None: + unsteady = results.get("Unsteady") + if unsteady is None: return False, "result HDF missing /Results/Unsteady group" + populated_datasets = 0 + + def _count_populated(_name, item): + nonlocal populated_datasets + if isinstance(item, h5py.Dataset) and item.size > 0: + populated_datasets += 1 + + unsteady.visititems(_count_populated) + if populated_datasets == 0: + return False, "result HDF has no populated /Results/Unsteady datasets" except Exception as e: return False, f"result HDF unreadable or invalid: {e}" return True, "ok" diff --git a/ras_commander/RasMap.py b/ras_commander/RasMap.py index 41e349f1c..c0dbde12a 100644 --- a/ras_commander/RasMap.py +++ b/ras_commander/RasMap.py @@ -69,6 +69,7 @@ - add_wse_comparison_layers(): Batch add WSE comparison layers for existing/proposed plan pairs """ +import os import re import subprocess import warnings @@ -4933,6 +4934,177 @@ def add_terrain_layer( logger.info("%s terrain layer '%s' in .rasmap", action, layer_name) logger.debug("Terrain layer '%s' filename: %s", layer_name, rel_path_str) + @staticmethod + @log_call + def prune_event_condition_layers( + keep_filenames: Sequence[Union[str, Path]], + rasmap_path: Optional[Union[str, Path]] = None, + ras_object=None, + backup: bool = True, + ) -> Dict[str, Any]: + """Remove stale event-condition layers from a cloned project. + + HEC-RAS evaluates ``RASEventConditions`` layers during unsteady + preprocessing, including unchecked layers inherited from other plans + or source-project folders. A cloned breakout project can therefore + fail with ``Error processing event conditions`` even when its current + plan and unsteady-flow file are valid. + + This method retains only exact filename matches, removes unmatched + event-condition layers anywhere in the ``.rasmap`` tree, writes + atomically, and validates exact readback. It never edits the referenced + HDF files. Pass an empty sequence to remove every event-condition + layer. + + Args: + keep_filenames: Exact layer ``Filename`` values to retain. Paths + inside the project folder may be supplied as absolute paths; + they are normalized to ``.\\`` relative Windows form. + rasmap_path: Explicit ``.rasmap`` path. When omitted, resolve it + from ``ras_object``. + ras_object: Optional initialized project object. + backup: Create a durable sibling backup before mutation. + + Returns: + JSON-safe mutation evidence including removed/retained layers, + backup path, and exact readback counts. + + Raises: + FileNotFoundError: If the ``.rasmap`` does not exist. + FileExistsError: If the backup or atomic staging path exists. + ValueError: If XML is invalid or a requested retained filename is + absent before mutation. + RuntimeError: If exact readback validation fails. + """ + ras_obj = ras_object or ras + if rasmap_path is None: + resolved = RasMap.get_rasmap_path(ras_obj) + if resolved is None: + raise FileNotFoundError("Project .rasmap file was not found") + target = Path(resolved) + else: + target = Path(rasmap_path) + if not target.is_file(): + raise FileNotFoundError(f"RASMapper file not found: {target}") + + def normalize_filename(value: Union[str, Path]) -> str: + raw = str(value).strip() + candidate = Path(raw) + if candidate.is_absolute(): + try: + raw = ".\\" + str( + RasUtils.safe_resolve(candidate).relative_to( + RasUtils.safe_resolve(target.parent) + ) + ) + except ValueError: + raw = str(candidate) + return raw.replace("/", "\\").casefold() + + requested = {normalize_filename(value) for value in keep_filenames} + try: + tree = ET.parse(target) + except ET.ParseError as exc: + raise ValueError(f"Error parsing .rasmap XML: {exc}") from exc + root = tree.getroot() + + layers = [] + for parent in root.iter(): + for layer in list(parent): + if ( + layer.tag == "Layer" + and layer.get("Type") == "RASEventConditions" + ): + layers.append((parent, layer)) + available = { + normalize_filename(layer.get("Filename") or "") + for _parent, layer in layers + } + missing = sorted(requested - available) + if missing: + raise ValueError( + "Requested event-condition layer filename(s) are absent: " + + ", ".join(missing) + ) + + removed = [] + retained = [] + for parent, layer in layers: + record = { + "name": layer.get("Name") or "", + "filename": layer.get("Filename") or "", + "checked": layer.get("Checked"), + } + if normalize_filename(record["filename"]) in requested: + retained.append(record) + else: + parent.remove(layer) + removed.append(record) + + backup_path = None + staging = target.with_name(f".{target.name}.event-conditions.tmp") + if removed: + if staging.exists(): + raise FileExistsError( + f"Atomic .rasmap staging path already exists: {staging}" + ) + if backup: + backup_path = target.with_name( + f"{target.stem}.event-conditions.bak{target.suffix}" + ) + if backup_path.exists(): + raise FileExistsError( + f"Event-condition backup already exists: {backup_path}" + ) + shutil.copy2(target, backup_path) + try: + tree.write(staging, encoding="utf-8", xml_declaration=False) + ET.parse(staging) + os.replace(staging, target) + finally: + if staging.exists(): + staging.unlink() + + readback_root = ET.parse(target).getroot() + readback = [] + for layer in readback_root.iter("Layer"): + if layer.get("Type") == "RASEventConditions": + readback.append( + { + "name": layer.get("Name") or "", + "filename": layer.get("Filename") or "", + "checked": layer.get("Checked"), + } + ) + readback_filenames = { + normalize_filename(layer["filename"]) for layer in readback + } + if readback_filenames != requested: + raise RuntimeError( + "Event-condition readback mismatch: expected " + f"{sorted(requested)}, observed {sorted(readback_filenames)}" + ) + + evidence = { + "rasmap_path": str(target), + "backup_path": str(backup_path) if backup_path is not None else None, + "requested_filenames": sorted(requested), + "before_count": len(layers), + "removed_count": len(removed), + "retained_count": len(retained), + "readback_count": len(readback), + "removed": removed, + "retained": retained, + "readback": readback, + "changed": bool(removed), + } + logger.info( + "Pruned %d stale event-condition layer(s); retained %d", + len(removed), + len(retained), + ) + return evidence + # ── Calculated Layers ─────────────────────────────────────────────── @staticmethod diff --git a/ras_commander/RasProcess.py b/ras_commander/RasProcess.py index b31b102f8..9777fa099 100644 --- a/ras_commander/RasProcess.py +++ b/ras_commander/RasProcess.py @@ -1318,6 +1318,271 @@ def _run_rasprocess( cwd=str(working_dir) if working_dir else None, ) + @staticmethod + def _geometry_completion_semantics(geom_hdf_path: Path) -> Dict[str, Any]: + """Inspect 1D completion layers and 2D property-table readiness.""" + evidence: Dict[str, Any] = { + "readable": False, + "has_1d_geometry": False, + "has_2d_geometry": False, + "edge_lines_written": False, + "interpolation_surface_written": False, + "two_d_flow_areas": {}, + "success": False, + "error": None, + } + try: + import h5py + + with h5py.File(geom_hdf_path, "r") as hdf: + evidence["readable"] = True + evidence["edge_lines_written"] = ( + "Geometry/River Edge Lines" in hdf + ) + evidence["interpolation_surface_written"] = ( + "Geometry/Cross Section Interpolation Surfaces" in hdf + ) + cross_sections = hdf.get("Geometry/Cross Sections/Attributes") + evidence["has_1d_geometry"] = bool( + cross_sections is not None and cross_sections.shape[0] > 0 + ) + + collection = hdf.get("Geometry/2D Flow Areas") + attributes = ( + collection.get("Attributes") if collection is not None else None + ) + evidence["has_2d_geometry"] = bool( + attributes is not None and attributes.shape[0] > 0 + ) + if collection is not None: + for name, area in collection.items(): + if name == "Attributes" or not isinstance(area, h5py.Group): + continue + cells = area.get("Cells Center Coordinate") + faces = area.get("Faces FacePoint Indexes") + cell_info = area.get("Cells Volume Elevation Info") + cell_values = area.get("Cells Volume Elevation Values") + face_info = area.get("Faces Area Elevation Info") + face_values = area.get("Faces Area Elevation Values") + cell_rows = int(cells.shape[0]) if cells is not None else 0 + face_rows = int(faces.shape[0]) if faces is not None else 0 + area_ready = bool( + cell_rows > 0 + and face_rows > 0 + and cell_info is not None + and cell_info.shape[0] == cell_rows + and cell_values is not None + and cell_values.shape[0] > 0 + and face_info is not None + and face_info.shape[0] == face_rows + and face_values is not None + and face_values.shape[0] > 0 + ) + evidence["two_d_flow_areas"][name] = { + "cell_rows": cell_rows, + "face_rows": face_rows, + "cell_property_info_rows": ( + int(cell_info.shape[0]) + if cell_info is not None + else 0 + ), + "cell_property_value_rows": ( + int(cell_values.shape[0]) + if cell_values is not None + else 0 + ), + "face_property_info_rows": ( + int(face_info.shape[0]) + if face_info is not None + else 0 + ), + "face_property_value_rows": ( + int(face_values.shape[0]) + if face_values is not None + else 0 + ), + "ready": area_ready, + } + + one_d_ready = ( + not evidence["has_1d_geometry"] + or evidence["edge_lines_written"] + ) + two_d_ready = ( + not evidence["has_2d_geometry"] + or ( + bool(evidence["two_d_flow_areas"]) + and all( + item["ready"] + for item in evidence["two_d_flow_areas"].values() + ) + ) + ) + evidence["success"] = bool( + (evidence["has_1d_geometry"] or evidence["has_2d_geometry"]) + and one_d_ready + and two_d_ready + ) + except Exception as exc: + evidence["error"] = f"{type(exc).__name__}: {exc}" + return evidence + + @staticmethod + def _terminate_owned_process_tree(process: subprocess.Popen) -> Dict[str, Any]: + """Terminate only the process tree rooted at *process*.""" + evidence = { + "root_pid": int(process.pid), + "observed_pids": [], + "terminated_pids": [], + "killed_pids": [], + "survivor_pids": [], + } + try: + import psutil + + try: + root = psutil.Process(int(process.pid)) + owned = [root, *root.children(recursive=True)] + except (psutil.NoSuchProcess, psutil.AccessDenied): + owned = [] + evidence["observed_pids"] = sorted(int(item.pid) for item in owned) + alive = [] + for item in reversed(owned): + try: + if item.is_running() and item.status() != psutil.STATUS_ZOMBIE: + evidence["terminated_pids"].append(int(item.pid)) + item.terminate() + alive.append(item) + except (psutil.NoSuchProcess, psutil.AccessDenied): + continue + _gone, alive = psutil.wait_procs(alive, timeout=3) + for item in alive: + try: + evidence["killed_pids"].append(int(item.pid)) + item.kill() + except (psutil.NoSuchProcess, psutil.AccessDenied): + continue + _gone, survivors = psutil.wait_procs(alive, timeout=3) + evidence["survivor_pids"] = sorted( + int(item.pid) for item in survivors + ) + except Exception: + if process.poll() is None: + process.kill() + try: + process.wait(timeout=10) + except Exception: + pass + return evidence + + @staticmethod + def _run_complete_geometry_supervised( + rasprocess_path: Path, + args: List[str], + geom_hdf_path: Path, + timeout: int, + working_dir: Path, + semantic_stable_seconds: float, + ) -> Tuple[subprocess.CompletedProcess, str, Dict[str, Any]]: + """Run CompleteGeometry until exit or stable semantic HDF completion.""" + if IS_LINUX: + wine_config = RasProcess._get_wine_config() + if wine_config is None: + raise RuntimeError( + "Wine not configured on Linux. Call RasProcess.configure_wine() " + "or set WINEPREFIX." + ) + command, env, cwd = RasProcess._build_wine_command( + rasprocess_path, + args, + wine_config, + working_dir=working_dir, + ) + else: + command = [str(rasprocess_path), *args] + env = None + cwd = str(working_dir) if working_dir else None + + if semantic_stable_seconds < 0: + raise ValueError("semantic_stable_seconds must be non-negative") + + with tempfile.TemporaryFile( + mode="w+", encoding="utf-8", errors="replace" + ) as stdout_file, tempfile.TemporaryFile( + mode="w+", encoding="utf-8", errors="replace" + ) as stderr_file: + process = subprocess.Popen( + command, + cwd=cwd, + env=env, + stdout=stdout_file, + stderr=stderr_file, + text=True, + ) + deadline = time.monotonic() + timeout + ready_since = None + ready_signature = None + completion_mode = "process_exit" + cleanup = { + "root_pid": int(process.pid), + "observed_pids": [int(process.pid)], + "terminated_pids": [], + "killed_pids": [], + "survivor_pids": [], + } + while process.poll() is None: + semantics = RasProcess._geometry_completion_semantics( + geom_hdf_path + ) + if semantics["success"]: + try: + stat = geom_hdf_path.stat() + signature = (stat.st_size, stat.st_mtime_ns) + except OSError: + signature = None + if signature is not None and signature == ready_signature: + ready_since = ready_since or time.monotonic() + else: + ready_signature = signature + ready_since = time.monotonic() + if ( + ready_since is not None + and time.monotonic() - ready_since + >= semantic_stable_seconds + ): + completion_mode = "semantic_hdf_stable" + cleanup = RasProcess._terminate_owned_process_tree(process) + break + else: + ready_since = None + ready_signature = None + + if time.monotonic() >= deadline: + cleanup = RasProcess._terminate_owned_process_tree(process) + raise subprocess.TimeoutExpired(command, timeout) + time.sleep(1) + + if process.poll() is None: + process.wait(timeout=10) + stdout_file.seek(0) + stderr_file.seek(0) + stdout = stdout_file.read() + stderr = stderr_file.read() + if IS_LINUX and stderr: + stderr = "\n".join( + line + for line in stderr.splitlines() + if not line.startswith(("0", "wine: ")) + and "err:" not in line[:20] + ) + completed = subprocess.CompletedProcess( + command, + process.returncode, + stdout, + stderr, + ) + return completed, completion_mode, cleanup + @staticmethod def _resolve_path_for_rasprocess(path: Path) -> str: """ @@ -1459,6 +1724,7 @@ def compute_geometry( ras_object=None, ras_version: str = None, timeout: int = 1800, + semantic_stable_seconds: float = 10.0, ) -> Dict[str, Any]: """ Run HEC-RAS's headless geometry completion (RasProcess.exe CompleteGeometry). @@ -1516,6 +1782,11 @@ def compute_geometry( Specific HEC-RAS version for the RasProcess.exe lookup. timeout : int, optional Command timeout in seconds (default 1800). + semantic_stable_seconds : float, optional + Under Wine, terminate only the owned ``RasProcess.exe`` tree after + the completed HDF semantics and file size/mtime remain stable for + this many seconds (default 10). Native Windows continues to wait + for normal process exit. Returns ------- @@ -1573,36 +1844,51 @@ def compute_geometry( ) logger.info(f"Running RasProcess.exe CompleteGeometry on {geom_hdf_path.name}") - result = RasProcess._run_rasprocess( - rasprocess, args, timeout=timeout, working_dir=geom_hdf_path.parent - ) + completion_mode = "process_exit" + owned_process_cleanup = None + if _is_wine_helper_runtime(): + result, completion_mode, owned_process_cleanup = ( + RasProcess._run_complete_geometry_supervised( + rasprocess, + args, + geom_hdf_path, + timeout=timeout, + working_dir=geom_hdf_path.parent, + semantic_stable_seconds=semantic_stable_seconds, + ) + ) + else: + result = RasProcess._run_rasprocess( + rasprocess, + args, + timeout=timeout, + working_dir=geom_hdf_path.parent, + ) stdout = result.stdout or "" stderr = result.stderr or "" - # Confirm the artifacts landed in the geometry HDF. - edge_lines_written = False - interp_surface_written = False - try: - import h5py - - with h5py.File(geom_hdf_path, "r") as hdf: - edge_lines_written = "Geometry/River Edge Lines" in hdf - interp_surface_written = ( - "Geometry/Cross Section Interpolation Surfaces" in hdf - ) - except Exception as e: - logger.debug(f"Post-run HDF inspection failed: {e}") + semantic_validation = RasProcess._geometry_completion_semantics( + geom_hdf_path + ) + edge_lines_written = semantic_validation["edge_lines_written"] + interp_surface_written = semantic_validation[ + "interpolation_surface_written" + ] + process_completed = ( + result.returncode == 0 or completion_mode == "semantic_hdf_stable" + ) success = ( - result.returncode == 0 + process_completed and "Error:" not in stdout and not stderr.lstrip().startswith("Error:") - and edge_lines_written + and semantic_validation["success"] ) if not success: logger.warning( f"CompleteGeometry did not fully succeed (rc={result.returncode}, " - f"edge_lines={edge_lines_written}). stdout: {stdout.strip()[:400]}" + f"mode={completion_mode}, semantics=" + f"{semantic_validation['success']}). stdout: {stdout.strip()[:400]}" ) return { @@ -1613,6 +1899,9 @@ def compute_geometry( "stderr": stderr, "edge_lines_written": edge_lines_written, "interpolation_surface_written": interp_surface_written, + "completion_mode": completion_mode, + "semantic_validation": semantic_validation, + "owned_process_cleanup": owned_process_cleanup, "success": success, } @@ -1624,6 +1913,7 @@ def complete_geometry( ras_object=None, ras_version: str = None, timeout: int = 1800, + semantic_stable_seconds: float = 10.0, ) -> Dict[str, Any]: """Deprecated alias for :meth:`compute_geometry`. Use ``compute_geometry``.""" warnings.warn( @@ -1638,6 +1928,7 @@ def complete_geometry( ras_object=ras_object, ras_version=ras_version, timeout=timeout, + semantic_stable_seconds=semantic_stable_seconds, ) @staticmethod diff --git a/ras_commander/RasUnsteady.py b/ras_commander/RasUnsteady.py index 4b7535676..4af14851c 100644 --- a/ras_commander/RasUnsteady.py +++ b/ras_commander/RasUnsteady.py @@ -41,6 +41,7 @@ def my_function(): - extract_tables() - write_table_to_file() - get_met_precipitation_config() +- disable_meteorology() - set_precipitation_hyetograph() - set_constant_precipitation() - set_gridded_precipitation() @@ -52,6 +53,7 @@ def my_function(): Precipitation Functions: - get_met_precipitation_config() - Read Meteorological Data tab precipitation settings +- disable_meteorology() - Disable all inherited meteorological forcing in text/HDF state - set_precipitation_hyetograph() - Write hyetograph DataFrame to unsteady file - set_gridded_precipitation() - Configure GDAL raster precipitation - configure_gridded_dss_precipitation() - Configure gridded DSS precipitation @@ -3659,6 +3661,225 @@ def set_met_precipitation_mode( mode_value, ) + @staticmethod + @log_call + def disable_meteorology( + unsteady_file: Union[str, Path], + *, + compiled_plan_hdf: Optional[Union[str, Path]] = None, + ras_object: Optional[Any] = None, + ) -> Dict[str, Any]: + """Disable inherited meteorological forcing in an unsteady-flow clone. + + HEC-RAS can retain meteorology state in the ``.u##.hdf`` sidecar even + when the corresponding ``.u##`` text contains no active precipitation, + wind, or evapotranspiration record. A cloned breakout plan can then + fail in-band while reading an enabled HDF variable whose ``Values`` + dataset is absent. This method makes the disabled state explicit in + the text file and unsteady sidecar without deleting sidecar data. For + an owned compiled ``*.tmp.hdf``, it removes the meteorology group to + match HEC-RAS's canonical no-meteorology execution state. + + The optional ``compiled_plan_hdf`` is limited to a results-free + ``*.tmp.hdf`` execution artifact. Supplying it is useful when Phase-1 + Windows/Wine preprocessing has already assembled the plan and the + caller must align that owned task-local artifact before native Linux + execution. Result ``*.hdf`` files are rejected. + + Parameters + ---------- + unsteady_file : str or Path + Unsteady-flow number or owned ``.u##`` path. + compiled_plan_hdf : str or Path, optional + Owned ``*.tmp.hdf`` whose ``/Event Conditions/Meteorology`` groups + should also be disabled. The file must already exist and contain + ``/Event Conditions``. + ras_object : optional + Project object used to resolve an unsteady-flow number. + + Returns + ------- + dict + Mutation evidence including the updated text keys, sidecar groups + whose ``Enabled`` attribute was set to zero, and any compiled-plan + meteorology groups removed. + """ + import h5py + + unsteady_path = RasUnsteady._resolve_unsteady_file_path( + unsteady_file, + ras_object=ras_object, + ) + sidecar_hdf = Path(str(unsteady_path) + ".hdf") + compiled_path: Optional[Path] = None + if compiled_plan_hdf is not None: + compiled_path = Path(compiled_plan_hdf) + if not compiled_path.name.casefold().endswith(".tmp.hdf"): + raise ValueError( + "compiled_plan_hdf must be an owned results-free '*.tmp.hdf' file" + ) + if not compiled_path.is_file(): + raise FileNotFoundError( + f"Compiled plan HDF not found: {compiled_path}" + ) + + hdf_targets = [ + path + for path in (sidecar_hdf, compiled_path) + if path is not None and path.exists() + ] + for hdf_path in hdf_targets: + with h5py.File(hdf_path, "r") as hdf_file: + if "Event Conditions" not in hdf_file: + raise ValueError( + f"Meteorology target lacks /Event Conditions: {hdf_path.name}" + ) + + RasUnsteady._replace_met_precipitation_keys( + unsteady_path, + [ + ("Precipitation Mode", "Disable"), + ("Met BC=Precipitation|Mode", "None"), + ], + ) + + with open( + unsteady_path, + "r", + encoding="utf-8", + errors="replace", + newline="", + ) as file: + lines = file.readlines() + newline = RasUnsteady._detect_line_ending(lines) + had_terminal_newline = bool( + lines and lines[-1].endswith(("\r\n", "\n", "\r")) + ) + managed_prefixes = ( + "Wind Mode=", + "Met BC=Evapotranspiration|Mode=", + ) + insert_index = next( + ( + index + for index, line in enumerate(lines) + if line.startswith(managed_prefixes) + ), + None, + ) + if insert_index is None: + precip_indexes = [ + index + for index, line in enumerate(lines) + if line.startswith(("Precipitation Mode=", "Met BC=Precipitation|")) + ] + insert_index = ( + max(precip_indexes) + 1 + if precip_indexes + else RasUnsteady._get_default_met_insert_index(lines) + ) + + filtered_lines: List[str] = [] + removed_before_insert = 0 + for index, line in enumerate(lines): + if line.startswith(managed_prefixes): + if index < insert_index: + removed_before_insert += 1 + continue + filtered_lines.append(line) + insert_index = max(0, insert_index - removed_before_insert) + desired_lines = [ + f"Wind Mode=Disable{newline}", + f"Met BC=Evapotranspiration|Mode=Disable{newline}", + ] + if ( + insert_index == len(filtered_lines) + and filtered_lines + and not filtered_lines[-1].endswith(("\n", "\r")) + ): + filtered_lines[-1] = f"{filtered_lines[-1]}{newline}" + if insert_index == len(filtered_lines) and not had_terminal_newline: + desired_lines[-1] = desired_lines[-1][:-len(newline)] + RasUnsteady._atomic_write_lines( + unsteady_path, + filtered_lines[:insert_index] + + desired_lines + + filtered_lines[insert_index:], + ) + + hdf_evidence = [] + for hdf_path in hdf_targets: + if compiled_path is not None and hdf_path == compiled_path: + with h5py.File(hdf_path, "r+") as hdf_file: + met_path = "Event Conditions/Meteorology" + removed_groups = [] + if met_path in hdf_file: + removed_groups = list(hdf_file[met_path].keys()) + del hdf_file[met_path] + hdf_file.flush() + hdf_evidence.append( + { + "path": str(hdf_path), + "representation": "compiled_plan", + "meteorology_group_removed": bool(removed_groups), + "removed_groups": removed_groups, + } + ) + continue + + disabled_groups = [] + with h5py.File(hdf_path, "r+") as hdf_file: + meteorology = hdf_file.require_group( + "Event Conditions/Meteorology" + ) + meteorology.require_group("Precipitation") + meteorology.require_group("Evapotranspiration") + for variable_name, item in meteorology.items(): + if not isinstance(item, h5py.Group): + continue + prior_value = item.attrs.get("Enabled") + item.attrs["Enabled"] = np.uint8(0) + disabled_groups.append( + { + "variable": variable_name, + "enabled_before": ( + int(prior_value) if prior_value is not None else None + ), + "enabled_after": int(item.attrs["Enabled"]), + } + ) + hdf_file.flush() + hdf_evidence.append( + { + "path": str(hdf_path), + "representation": "unsteady_sidecar", + "disabled_groups": disabled_groups, + } + ) + + ras_obj = ras_object or ras + if ras_obj is not None and hasattr(ras_obj, "get_unsteady_entries"): + try: + ras_obj.unsteady_df = ras_obj.get_unsteady_entries() + except Exception as exc: + logger.debug("unsteady_df refresh skipped: %s", exc) + + logger.info( + "Disabled meteorology in %s and %d HDF artifact(s)", + unsteady_path.name, + len(hdf_evidence), + ) + return { + "unsteady_file": str(unsteady_path), + "text_state": { + "precipitation_mode": "Disable", + "precipitation_met_mode": "None", + "wind_mode": "Disable", + "evapotranspiration_mode": "Disable", + }, + "hdf_targets": hdf_evidence, + } + @staticmethod @log_call def set_hydrograph_fixed_start_time( diff --git a/ras_commander/_gdal_runtime.py b/ras_commander/_gdal_runtime.py index f1b75d0d4..98cdced5b 100644 --- a/ras_commander/_gdal_runtime.py +++ b/ras_commander/_gdal_runtime.py @@ -215,9 +215,22 @@ def ensure_python_gdal_junction( continue if result.returncode == 0: - logger.info("Created GDAL junction for HEC-RAS GDAL bridge") - logger.debug("Created GDAL junction: %s -> %s", gdal_junction, paths.gdal_root) - ok = python_gdal_bridge_is_usable(target_python_dir) and ok + usable = python_gdal_bridge_is_usable(target_python_dir) + if usable: + logger.info("Created GDAL junction for HEC-RAS GDAL bridge") + logger.debug( + "Created GDAL junction: %s -> %s", + gdal_junction, + paths.gdal_root, + ) + else: + logger.warning( + "GDAL junction command reported success but %s is not " + "usable; Wine-hosted callers should prepare this " + "task-local bridge from the native Linux host", + gdal_junction, + ) + ok = usable and ok continue logger.error( diff --git a/ras_commander/geom/GeomPreprocessor.py b/ras_commander/geom/GeomPreprocessor.py index 55c7d4530..84ab2d67d 100644 --- a/ras_commander/geom/GeomPreprocessor.py +++ b/ras_commander/geom/GeomPreprocessor.py @@ -272,10 +272,13 @@ def run_geometry_preprocessor( _watchdog.add_pid(process.pid) geom_only_artifacts = None + geometry_hdf_path = ( + project_folder / f"{project_name}.g{geometry_number}.hdf" + ) if geometry_only: geom_only_artifacts = [ project_folder / f"{project_name}.c{geometry_number}", - project_folder / f"{project_name}.g{geometry_number}.hdf", + geometry_hdf_path, project_folder / f"{project_name}.x{geometry_number}", ] @@ -287,6 +290,8 @@ def run_geometry_preprocessor( signals=flow_start_signals or GEOMETRY_PREPROCESSOR_FLOW_START_SIGNALS, geometry_only_artifacts=geom_only_artifacts, + geometry_hdf_path=geometry_hdf_path, + flow_type=flow_type, ) timed_out = monitor_result["timed_out"] @@ -370,6 +375,17 @@ def run_geometry_preprocessor( "No geometry preprocessor artifacts found " f"(.c{geometry_number}, .g{geometry_number}.hdf, .x{geometry_number}, or .b{plan_num})" ) + geometry_hdf_ready, geometry_hdf_reason = ( + GeomPreprocessor._geometry_hdf_readiness( + geometry_hdf_path, + flow_type=flow_type, + ) + ) + if geometry_only and not geometry_hdf_ready: + errors.append( + "Geometry HDF is not semantically ready: " + f"{geometry_hdf_reason}" + ) if not existing_message_paths and not (geometry_only and artifact_paths): errors.append("No compute messages were produced") @@ -562,14 +578,18 @@ def _monitor_compute_messages( max_wait: int, signals: List[str], geometry_only_artifacts: Optional[List[Path]] = None, + geometry_hdf_path: Optional[Path] = None, + flow_type: Optional[str] = None, ) -> dict: """Poll compute-message files until a flow-start signal, exit, or timeout. When *geometry_only_artifacts* is provided (geometry_only mode), the loop also checks whether preprocessing artifacts have been freshly - written and stabilized. Because ``Run UNet=0`` prevents any - flow-start signal from appearing, artifact stability is the only - reliable completion indicator in geometry-only mode. + written, stabilized, and the compiled geometry HDF is semantically + populated. Because ``Run UNet=0`` prevents any flow-start signal + from appearing, semantic artifact readiness is the reliable completion + indicator in geometry-only mode. An empty HDF placeholder must never + terminate HEC-RAS early. """ positions = {path: 0 for path in message_paths} signal_detected = None @@ -595,6 +615,18 @@ def _monitor_compute_messages( lower_chunk = chunk.lower() for signal, lower_signal in signal_patterns: if lower_signal in lower_chunk: + if geometry_hdf_path is not None: + ready, reason = GeomPreprocessor._geometry_hdf_readiness( + geometry_hdf_path, + flow_type=flow_type, + ) + if not ready: + logger.debug( + "Ignoring flow-start signal until geometry " + "HDF is ready: %s", + reason, + ) + continue signal_detected = signal return { "signal_detected": signal_detected, @@ -619,9 +651,21 @@ def _monitor_compute_messages( if fresh: try: latest_mtime = max(a.stat().st_mtime for a in fresh) - if time.time() - latest_mtime >= _ARTIFACT_STABLE: + hdf_ready = True + if geometry_hdf_path is not None: + hdf_ready, _reason = ( + GeomPreprocessor._geometry_hdf_readiness( + geometry_hdf_path, + flow_type=flow_type, + ) + ) + if ( + hdf_ready + and time.time() - latest_mtime >= _ARTIFACT_STABLE + ): signal_detected = ( - "Geometry preprocessing artifacts stable " + "Geometry preprocessing artifacts semantically " + "ready and stable " "(geometry_only mode)" ) return { @@ -638,6 +682,88 @@ def _monitor_compute_messages( return {"signal_detected": signal_detected, "timed_out": True} + @staticmethod + def _geometry_hdf_readiness( + geometry_hdf_path: Path, + flow_type: Optional[str] = None, + ) -> tuple[bool, str]: + """Return whether a compiled geometry HDF contains usable geometry. + + HEC-RAS creates a small HDF placeholder containing only an empty + ``/Geometry`` group before it imports the text geometry. File + existence, non-zero size, and a quiet modification time therefore do + not prove preprocessing is complete. + + For a classified 2D plan, require at least one flow area with non-empty + cell and face topology. For other plan types, require a recognized, + non-empty geometry feature collection. The check is read-only and + returns a diagnostic instead of raising while HEC-RAS still owns the + file. + """ + path = Path(geometry_hdf_path) + if not path.is_file(): + return False, f"missing file: {path}" + try: + if path.stat().st_size <= 0: + return False, f"empty file: {path}" + except OSError as exc: + return False, f"could not stat {path}: {exc}" + + try: + import h5py + + with h5py.File(path, "r") as hdf: + if not hdf.attrs.get("File Type"): + return False, "root File Type attribute is absent" + if "Geometry" not in hdf: + return False, "/Geometry group is absent" + + normalized_flow_type = str(flow_type or "").casefold() + if "2d" in normalized_flow_type: + collection_path = "Geometry/2D Flow Areas" + if collection_path not in hdf: + return False, "/Geometry/2D Flow Areas is absent" + collection = hdf[collection_path] + attributes = collection.get("Attributes") + if attributes is None or attributes.shape[0] == 0: + return False, "2D Flow Areas/Attributes is empty" + + ready_areas = [] + for name, item in collection.items(): + if name == "Attributes" or not isinstance(item, h5py.Group): + continue + cells = item.get("Cells Center Coordinate") + faces = item.get("Faces FacePoint Indexes") + if ( + cells is not None + and faces is not None + and cells.shape[0] > 0 + and faces.shape[0] > 0 + ): + ready_areas.append(name) + if not ready_areas: + return False, "no 2D flow area has non-empty cell/face topology" + return True, f"2D geometry ready: {', '.join(ready_areas)}" + + collection_paths = ( + "Geometry/Cross Sections/Attributes", + "Geometry/River Centerlines/Attributes", + "Geometry/Storage Areas/Attributes", + "Geometry/Structures/Attributes", + "Geometry/SA 2D Area Conn/Attributes", + "Geometry/2D Flow Areas/Attributes", + ) + populated = [] + for collection_path in collection_paths: + dataset = hdf.get(collection_path) + if dataset is not None and getattr(dataset, "shape", (0,))[0] > 0: + populated.append(collection_path) + if not populated: + return False, "no recognized geometry feature collection is populated" + return True, f"geometry ready: {', '.join(populated)}" + except Exception as exc: + return False, f"could not inspect geometry HDF: {type(exc).__name__}: {exc}" + @staticmethod def _wait_for_preprocess_child( tmp_hdf_path: Path, diff --git a/tests/test_geom_preprocessor.py b/tests/test_geom_preprocessor.py index 42374ead0..e2d760195 100644 --- a/tests/test_geom_preprocessor.py +++ b/tests/test_geom_preprocessor.py @@ -52,6 +52,26 @@ def _write_legacy_geometry_hdf( return path +def _write_ready_2d_geometry_hdf(path: Path) -> Path: + with h5py.File(path, "w") as hdf: + hdf.attrs["File Type"] = "HEC-RAS Geometry" + collection = hdf.create_group("Geometry/2D Flow Areas") + collection.create_dataset( + "Attributes", + data=np.array([(b"MainArea",)], dtype=[("Name", "S32")]), + ) + mesh = collection.create_group("MainArea") + mesh.create_dataset( + "Cells Center Coordinate", + data=np.array([[1.0, 2.0]], dtype=np.float64), + ) + mesh.create_dataset( + "Faces FacePoint Indexes", + data=np.array([[0, 1]], dtype=np.int32), + ) + return path + + def test_compute_message_paths_include_data_error_files(tmp_path): paths = GeomPreprocessor._compute_message_paths(tmp_path, "Model", "04") @@ -126,6 +146,52 @@ def test_preprocessor_artifacts_include_fresh_tmp_hdf_only(tmp_path): assert artifacts == [fresh_tmp_hdf] +def test_geometry_hdf_readiness_rejects_empty_placeholder(tmp_path): + hdf_path = tmp_path / "Model.g03.hdf" + with h5py.File(hdf_path, "w") as hdf: + hdf.create_group("Geometry") + + ready, reason = GeomPreprocessor._geometry_hdf_readiness( + hdf_path, + flow_type="unsteady_2d", + ) + + assert not ready + assert "File Type" in reason + + +def test_geometry_hdf_readiness_rejects_2d_area_without_topology(tmp_path): + hdf_path = tmp_path / "Model.g03.hdf" + with h5py.File(hdf_path, "w") as hdf: + hdf.attrs["File Type"] = "HEC-RAS Geometry" + collection = hdf.create_group("Geometry/2D Flow Areas") + collection.create_dataset( + "Attributes", + data=np.array([(b"MainArea",)], dtype=[("Name", "S32")]), + ) + collection.create_group("MainArea") + + ready, reason = GeomPreprocessor._geometry_hdf_readiness( + hdf_path, + flow_type="unsteady_2d", + ) + + assert not ready + assert "cell/face topology" in reason + + +def test_geometry_hdf_readiness_accepts_populated_2d_mesh(tmp_path): + hdf_path = _write_ready_2d_geometry_hdf(tmp_path / "Model.g03.hdf") + + ready, reason = GeomPreprocessor._geometry_hdf_readiness( + hdf_path, + flow_type="unsteady_2d", + ) + + assert ready + assert "MainArea" in reason + + def test_geometry_only_run_flags_disable_unsteady_flow(tmp_path): plan_path = tmp_path / "Model.p01" plan_path.write_text( diff --git a/tests/test_linux_robustness.py b/tests/test_linux_robustness.py index e6a4b4958..6e180bbc6 100644 --- a/tests/test_linux_robustness.py +++ b/tests/test_linux_robustness.py @@ -2,6 +2,7 @@ import os import h5py +import numpy as np import pytest from ras_commander.RasCmdr import RasCmdr @@ -18,7 +19,8 @@ def _make_hdf(path, *, results=True, unsteady=True): if results: r = hf.create_group("Results") if unsteady: - r.create_group("Unsteady") + u = r.create_group("Unsteady") + u.create_dataset("Output", data=[1.0]) # --- CLB-882: validate solve beyond exit-code 0 --- @@ -43,6 +45,55 @@ def test_validate_linux_solve_fails_on_in_band_error(tmp_path): assert "solver log reports failure" in reason +@pytest.mark.parametrize( + "message", + [ + " ERROR: READ_UN_MET_EVAPO_DATA: Evapotranspiration values not found\n", + "HDF_ERROR with the Geometry\n", + ], +) +def test_validate_linux_solve_fails_on_explicit_native_error_lines( + tmp_path, + message, +): + log = tmp_path / "compute_linux_01.log" + _write_log(log, message) + hdf = tmp_path / "p01.hdf" + _make_hdf(hdf, results=True, unsteady=True) + + ok, reason = RasCmdr._validate_linux_solve(log, hdf, "01") + + assert ok is False + assert "solver log reports failure" in reason + + +def test_validate_linux_solve_requires_finished_banner(tmp_path): + log = tmp_path / "compute_linux_01.log" + _write_log(log, "Starting Unsteady Flow Computations\n") + hdf = tmp_path / "p01.hdf" + _make_hdf(hdf, results=True, unsteady=True) + + ok, reason = RasCmdr._validate_linux_solve(log, hdf, "01") + + assert ok is False + assert "Finished Unsteady Flow Simulation" in reason + + +def test_validate_linux_solve_ignores_volume_accounting_error_label(tmp_path): + log = tmp_path / "compute_linux_01.log" + _write_log( + log, + "Overall Volume Accounting Error as percentage: 0.0001\n" + "Finished Unsteady Flow Simulation\n", + ) + hdf = tmp_path / "p01.hdf" + _make_hdf(hdf, results=True, unsteady=True) + + ok, reason = RasCmdr._validate_linux_solve(log, hdf, "01") + + assert ok is True, reason + + def test_validate_linux_solve_fails_when_no_results_group(tmp_path): log = tmp_path / "compute_linux_01.log" _write_log(log, "Finished Unsteady Flow Simulation\n") @@ -64,12 +115,68 @@ def test_validate_linux_solve_fails_when_results_but_no_unsteady(tmp_path): assert "Unsteady" in reason +def test_validate_linux_solve_fails_when_unsteady_has_no_populated_datasets(tmp_path): + log = tmp_path / "compute_linux_01.log" + _write_log(log, "Finished Unsteady Flow Simulation\n") + hdf = tmp_path / "p01.hdf" + with h5py.File(hdf, "w") as hf: + hf.create_group("Results/Unsteady") + + ok, reason = RasCmdr._validate_linux_solve(log, hdf, "01") + + assert ok is False + assert "no populated" in reason + + def test_validate_linux_solve_fails_on_unreadable_log(tmp_path): ok, reason = RasCmdr._validate_linux_solve(tmp_path / "missing.log", tmp_path / "x.hdf", "01") assert ok is False assert "log" in reason.lower() +def test_effective_linux_core_count_caps_to_affinity(monkeypatch): + monkeypatch.setattr( + os, + "sched_getaffinity", + lambda _pid: {0, 2, 4}, + raising=False, + ) + + assert RasCmdr._effective_linux_core_count(8) == 3 + assert RasCmdr._effective_linux_core_count(2) == 2 + + +@pytest.mark.parametrize("invalid", [True, 0, -1, 1.5, "4"]) +def test_effective_linux_core_count_rejects_invalid_values(invalid): + with pytest.raises(ValueError, match="positive integer"): + RasCmdr._effective_linux_core_count(invalid) + + +def test_set_linux_hdf_num_cores_updates_1d_and_each_2d_mesh(tmp_path): + tmp_hdf = tmp_path / "Model.p01.tmp.hdf" + with h5py.File(tmp_hdf, "w") as hdf: + parameters = hdf.require_group("Plan Data/Plan Parameters") + parameters.attrs["1D Cores"] = np.int32(2) + parameters.attrs["2D Cores (per mesh)"] = np.array([2, 1], dtype=np.int32) + + evidence = RasCmdr._set_linux_hdf_num_cores(tmp_hdf, 3) + + with h5py.File(tmp_hdf, "r") as hdf: + parameters = hdf["Plan Data/Plan Parameters"] + assert int(parameters.attrs["1D Cores"]) == 3 + assert list(parameters.attrs["2D Cores (per mesh)"]) == [3, 3] + assert len(evidence["updated_attributes"]) == 2 + assert evidence["effective_cores"] == 3 + + +def test_set_linux_hdf_num_cores_rejects_result_hdf(tmp_path): + result_hdf = tmp_path / "Model.p01.hdf" + _make_hdf(result_hdf) + + with pytest.raises(ValueError, match=r"\*\.tmp\.hdf"): + RasCmdr._set_linux_hdf_num_cores(result_hdf, 2) + + # --- CLB-883: native Linux install discovery --- def _make_native_root(tmp_path): diff --git a/tests/test_ras_geometry_compute_unit.py b/tests/test_ras_geometry_compute_unit.py index 76a06428b..5fe022d09 100644 --- a/tests/test_ras_geometry_compute_unit.py +++ b/tests/test_ras_geometry_compute_unit.py @@ -219,6 +219,104 @@ def test_rasprocess_compute_geometry_rejects_stderr_error(monkeypatch, tmp_path) assert result["success"] is False +def _write_complete_2d_geometry(path, include_face_values=True): + import h5py + + with h5py.File(path, "w") as hdf: + collection = hdf.create_group("Geometry/2D Flow Areas") + collection.create_dataset("Attributes", data=[1]) + area = collection.create_group("Area") + area.create_dataset( + "Cells Center Coordinate", + data=[[0.0, 0.0], [1.0, 1.0]], + ) + area.create_dataset("Faces FacePoint Indexes", data=[[0, 1]]) + area.create_dataset("Cells Volume Elevation Info", data=[[0, 1], [1, 1]]) + area.create_dataset( + "Cells Volume Elevation Values", + data=[[0.0, 0.0], [1.0, 1.0]], + ) + area.create_dataset("Faces Area Elevation Info", data=[[0, 1]]) + if include_face_values: + area.create_dataset( + "Faces Area Elevation Values", + data=[[0.0, 1.0, 1.0, 0.04]], + ) + return path + + +def test_geometry_completion_semantics_accepts_complete_2d_tables(tmp_path): + geom = _write_complete_2d_geometry(tmp_path / "model.g01.hdf") + + evidence = RasProcess._geometry_completion_semantics(geom) + + assert evidence["success"] is True + assert evidence["has_2d_geometry"] is True + assert evidence["has_1d_geometry"] is False + assert evidence["two_d_flow_areas"]["Area"]["ready"] is True + + +def test_geometry_completion_semantics_rejects_incomplete_2d_tables(tmp_path): + geom = _write_complete_2d_geometry( + tmp_path / "model.g01.hdf", + include_face_values=False, + ) + + evidence = RasProcess._geometry_completion_semantics(geom) + + assert evidence["success"] is False + assert evidence["two_d_flow_areas"]["Area"]["ready"] is False + + +def test_compute_geometry_accepts_supervised_wine_semantic_completion( + monkeypatch, + tmp_path, +): + geom = _write_complete_2d_geometry(tmp_path / "model.g01.hdf") + rasprocess_exe = tmp_path / "RasProcess.exe" + rasprocess_exe.touch() + cleanup = { + "root_pid": 42, + "observed_pids": [42], + "terminated_pids": [42], + "killed_pids": [], + "survivor_pids": [], + } + + monkeypatch.setattr( + RasProcess, + "find_rasprocess", + staticmethod(lambda version=None: rasprocess_exe), + ) + rasprocess_module = __import__( + "ras_commander.RasProcess", + fromlist=["_is_wine_helper_runtime"], + ) + monkeypatch.setattr( + rasprocess_module, + "_is_wine_helper_runtime", + lambda: True, + ) + monkeypatch.setattr( + RasProcess, + "_run_complete_geometry_supervised", + staticmethod( + lambda *args, **kwargs: ( + SimpleNamespace(returncode=15, stdout="", stderr=""), + "semantic_hdf_stable", + cleanup, + ) + ), + ) + + result = RasProcess.compute_geometry(geom, ras_version="6.6") + + assert result["success"] is True + assert result["completion_mode"] == "semantic_hdf_stable" + assert result["semantic_validation"]["success"] is True + assert result["owned_process_cleanup"] == cleanup + + def _diff_frames(left, channel, right): import geopandas as gpd from shapely.geometry import LineString diff --git a/tests/test_rasmap_event_conditions.py b/tests/test_rasmap_event_conditions.py new file mode 100644 index 000000000..2247cafd9 --- /dev/null +++ b/tests/test_rasmap_event_conditions.py @@ -0,0 +1,86 @@ +from pathlib import Path +import xml.etree.ElementTree as ET + +import pytest + +from ras_commander import RasMap + + +def _write_rasmap(path: Path) -> Path: + path.write_text( + """ + + + + + + + + + + +""", + encoding="utf-8", + ) + return path + + +def _event_filenames(path: Path) -> list[str]: + root = ET.parse(path).getroot() + return [ + layer.get("Filename") or "" + for layer in root.iter("Layer") + if layer.get("Type") == "RASEventConditions" + ] + + +def test_prune_event_conditions_is_recursive_atomic_and_backed_up(tmp_path): + rasmap = _write_rasmap(tmp_path / "Model.rasmap") + current = tmp_path / "Model.u09.hdf" + current.touch() + + evidence = RasMap.prune_event_condition_layers( + [current], + rasmap_path=rasmap, + ) + + assert evidence["before_count"] == 3 + assert evidence["removed_count"] == 2 + assert evidence["retained_count"] == 1 + assert evidence["readback_count"] == 1 + assert evidence["changed"] is True + assert _event_filenames(rasmap) == [r".\Model.u09.hdf"] + backup = tmp_path / "Model.event-conditions.bak.rasmap" + assert backup.is_file() + assert len(_event_filenames(backup)) == 3 + assert not (tmp_path / ".Model.rasmap.event-conditions.tmp").exists() + + +def test_prune_event_conditions_missing_requested_filename_is_non_mutating(tmp_path): + rasmap = _write_rasmap(tmp_path / "Model.rasmap") + before = rasmap.read_bytes() + + with pytest.raises(ValueError, match=r"filename\(s\) are absent"): + RasMap.prune_event_condition_layers( + [r".\Model.u77.hdf"], + rasmap_path=rasmap, + ) + + assert rasmap.read_bytes() == before + assert not (tmp_path / "Model.event-conditions.bak.rasmap").exists() + + +def test_prune_event_conditions_refuses_existing_backup(tmp_path): + rasmap = _write_rasmap(tmp_path / "Model.rasmap") + backup = tmp_path / "Model.event-conditions.bak.rasmap" + backup.write_text("existing", encoding="utf-8") + before = rasmap.read_bytes() + + with pytest.raises(FileExistsError, match="backup already exists"): + RasMap.prune_event_condition_layers( + [r".\Model.u09.hdf"], + rasmap_path=rasmap, + ) + + assert rasmap.read_bytes() == before + assert backup.read_text(encoding="utf-8") == "existing" diff --git a/tests/test_rasunsteady_disable_meteorology.py b/tests/test_rasunsteady_disable_meteorology.py new file mode 100644 index 000000000..fb9f3d653 --- /dev/null +++ b/tests/test_rasunsteady_disable_meteorology.py @@ -0,0 +1,156 @@ +from pathlib import Path + +import h5py +import numpy as np +import pytest + +from ras_commander import RasUnsteady + + +def _write_event_hdf(path: Path, *, include_event_conditions: bool = True) -> Path: + with h5py.File(path, "w") as hdf: + if not include_event_conditions: + hdf.create_group("Geometry") + return path + met = hdf.require_group("Event Conditions/Meteorology") + met.require_group("Evapotranspiration").attrs["Enabled"] = np.uint8(1) + met.require_group("Precipitation").attrs["Enabled"] = np.uint8(1) + met.require_group("Wind Speed").attrs["Enabled"] = np.uint8(1) + met["Evapotranspiration"].create_dataset("Values", data=[1.0, 2.0]) + return path + + +def test_disable_meteorology_updates_text_sidecar_and_compiled_tmp(tmp_path): + unsteady = tmp_path / "Model.u01" + unsteady.write_bytes( + b"Flow Title=Clone\r\n" + b"Program Version=6.60\r\n" + b"Precipitation Mode=Enable\r\n" + b"Met BC=Precipitation|Mode=Gridded\r\n" + b"Met BC=Precipitation|Gridded Source=DSS\r\n" + b"Wind Mode=Enable\r\n" + b"Met BC=Evapotranspiration|Mode=Point Gage\r\n" + b"Met BC=Evapotranspiration|Point Time Series=Station A\r\n" + b"Boundary Location=,,,,,Area,,BC,\r\n" + ) + sidecar = _write_event_hdf(Path(str(unsteady) + ".hdf")) + compiled = _write_event_hdf(tmp_path / "Model.p01.tmp.hdf") + + evidence = RasUnsteady.disable_meteorology( + unsteady, + compiled_plan_hdf=compiled, + ) + + raw = unsteady.read_bytes() + assert raw.count(b"Precipitation Mode=Disable\r\n") == 1 + assert raw.count(b"Met BC=Precipitation|Mode=None\r\n") == 1 + assert raw.count(b"Wind Mode=Disable\r\n") == 1 + assert raw.count(b"Met BC=Evapotranspiration|Mode=Disable\r\n") == 1 + assert b"Met BC=Precipitation|Gridded Source=" not in raw + assert b"Met BC=Evapotranspiration|Point Time Series=Station A\r\n" in raw + assert raw.endswith(b"\r\n") + assert len(evidence["hdf_targets"]) == 2 + + with h5py.File(sidecar, "r") as hdf: + met = hdf["Event Conditions/Meteorology"] + assert all( + int(item.attrs["Enabled"]) == 0 + for item in met.values() + if isinstance(item, h5py.Group) + ) + assert list(met["Evapotranspiration/Values"][()]) == [1.0, 2.0] + + with h5py.File(compiled, "r") as hdf: + assert "Event Conditions/Meteorology" not in hdf + compiled_evidence = next( + item + for item in evidence["hdf_targets"] + if item["representation"] == "compiled_plan" + ) + assert compiled_evidence["meteorology_group_removed"] is True + assert sorted(compiled_evidence["removed_groups"]) == [ + "Evapotranspiration", + "Precipitation", + "Wind Speed", + ] + + +def test_disable_meteorology_is_idempotent(tmp_path): + unsteady = tmp_path / "Model.u01" + unsteady.write_text( + "Flow Title=Clone\n" + "Program Version=6.60\n" + "Boundary Location=,,,,,Area,,BC,", + encoding="utf-8", + ) + + RasUnsteady.disable_meteorology(unsteady) + first = unsteady.read_bytes() + RasUnsteady.disable_meteorology(unsteady) + + assert unsteady.read_bytes() == first + lines = unsteady.read_text(encoding="utf-8").splitlines() + assert lines.count("Precipitation Mode=Disable") == 1 + assert lines.count("Met BC=Precipitation|Mode=None") == 1 + assert lines.count("Wind Mode=Disable") == 1 + assert lines.count("Met BC=Evapotranspiration|Mode=Disable") == 1 + + +def test_disable_meteorology_compiled_target_is_idempotent(tmp_path): + unsteady = tmp_path / "Model.u01" + unsteady.write_text( + "Flow Title=Clone\nProgram Version=6.60\n", + encoding="utf-8", + ) + compiled = _write_event_hdf(tmp_path / "Model.p01.tmp.hdf") + + first = RasUnsteady.disable_meteorology( + unsteady, + compiled_plan_hdf=compiled, + ) + second = RasUnsteady.disable_meteorology( + unsteady, + compiled_plan_hdf=compiled, + ) + + assert first["hdf_targets"][-1]["meteorology_group_removed"] is True + assert second["hdf_targets"][-1]["meteorology_group_removed"] is False + with h5py.File(compiled, "r") as hdf: + assert "Event Conditions/Meteorology" not in hdf + + +@pytest.mark.parametrize("filename", ["Model.p01.hdf", "Model.u01.hdf"]) +def test_disable_meteorology_rejects_non_tmp_compiled_target_without_mutation( + tmp_path, + filename, +): + unsteady = tmp_path / "Model.u01" + original = b"Flow Title=Clone\nProgram Version=6.60\n" + unsteady.write_bytes(original) + target = _write_event_hdf(tmp_path / filename) + + with pytest.raises(ValueError, match=r"\*\.tmp\.hdf"): + RasUnsteady.disable_meteorology( + unsteady, + compiled_plan_hdf=target, + ) + + assert unsteady.read_bytes() == original + + +def test_disable_meteorology_rejects_tmp_without_event_conditions(tmp_path): + unsteady = tmp_path / "Model.u01" + original = b"Flow Title=Clone\nProgram Version=6.60\n" + unsteady.write_bytes(original) + target = _write_event_hdf( + tmp_path / "Model.p01.tmp.hdf", + include_event_conditions=False, + ) + + with pytest.raises(ValueError, match="lacks /Event Conditions"): + RasUnsteady.disable_meteorology( + unsteady, + compiled_plan_hdf=target, + ) + + assert unsteady.read_bytes() == original