Terrain Modifications: High-Ground and Polygon Writer Validation¶
This notebook validates ras-commander terrain modification writers with a real HEC-RAS example project. It creates copied terrain sidecar modifications for high-ground lines, boundary-sampled polygon multipoint grading, and shape_z polygon grading, verifies raster readback, runs HEC-RAS through RasCmdr.compute_plan(), and compares maximum water surface elevation impacts.
Workflow¶
- Extract isolated copies of the BaldEagleCrkMulti2D example project.
- Add a
take_higherhigh-ground line modification with a 6 ft raise to the proposed copy. - Add a boundary-from-terrain polygon multipoint grading area with interior control points set to existing grade plus 5 ft and plus 3 ft.
- Add a second polygon example using
elevation_method='shape_z'with ring Z values set to sampled existing grade plus 4 ft. - Verify HDF and
.rasmapmetadata for both polygon writer paths. - Plot writer-computed and raster-read before/after terrain profiles.
- Optionally export one bounded 2x GeoTIFF with the high-ground modification disabled and enabled through native RAS Mapper.
- Map modified terrain at geometry cell centers.
- Run base, high-ground, multipoint polygon, and shape_z polygon plan 18 with geometry preprocessing forced.
- Compare maximum WSE results from the computed plan HDF files and identify benefit/rise polygons with
HdfBenefitAreas.
Python
# Development-source setup and imports
from pathlib import Path
import json
import logging
import os
import shutil
import sys
import warnings
import xml.etree.ElementTree as ET
import h5py
import matplotlib.pyplot as plt
from matplotlib.path import Path as MplPath
import matplotlib.tri as mtri
import numpy as np
import pandas as pd
candidates = [Path.cwd(), Path.cwd().parent]
REPO_ROOT = next(path for path in candidates if (path / 'ras_commander').exists())
if str(REPO_ROOT) not in sys.path:
sys.path.insert(0, str(REPO_ROOT))
from ras_commander import (
HdfBenefitAreas,
HdfMesh,
HdfResultsMesh,
RasCmdr,
RasExamples,
RasMap,
RasTerrain,
TerrainExportResult,
init_ras_project,
)
from ras_commander.terrain import RasTerrainMod, RasTerrainModification, RasTerrainModWriter
logging.getLogger('ras_commander').setLevel(logging.WARNING)
warnings.filterwarnings('ignore', category=UserWarning)
pd.options.display.max_columns = 20
NOTEBOOK_DIR = REPO_ROOT / 'examples'
DEFAULT_EXAMPLE_RUN_ROOT = NOTEBOOK_DIR / 'example_projects'
EXAMPLE_RUN_ROOT = Path(
os.environ.get('RAS_COMMANDER_EXAMPLE_RUN_ROOT', DEFAULT_EXAMPLE_RUN_ROOT)
).expanduser().resolve()
ARTIFACT_ROOT = EXAMPLE_RUN_ROOT / '316_terrain_modifications'
ARTIFACT_ROOT.mkdir(parents=True, exist_ok=True)
RUN_NATIVE_TERRAIN_EXPORT = False
if os.environ.get('RAS_COMMANDER_RUN_NATIVE_TERRAIN_EXPORT') == '1':
RUN_NATIVE_TERRAIN_EXPORT = True
OVERWRITE_NATIVE_TERRAIN_EXPORT = False
print(f'Repository: {REPO_ROOT.name}')
print(f'Run root: {EXAMPLE_RUN_ROOT.name}')
print(f'Artifacts: {ARTIFACT_ROOT.name}')
Text Only
Repository: G:\GH\ras-commander-pr-docs
Artifacts: examples\example_projects\316_terrain_modifications
Python
# Extract isolated base/proposed/polygon project copies
template = RasExamples.extract_project(
'BaldEagleCrkMulti2D',
output_path=ARTIFACT_ROOT,
suffix='template',
)
base_project = ARTIFACT_ROOT / 'BaldEagleCrkMulti2D_base'
proposed_project = ARTIFACT_ROOT / 'BaldEagleCrkMulti2D_proposed'
polygon_project = ARTIFACT_ROOT / 'BaldEagleCrkMulti2D_polygon_multipoint'
shape_z_project = ARTIFACT_ROOT / 'BaldEagleCrkMulti2D_polygon_shape_z'
for project in (base_project, proposed_project, polygon_project, shape_z_project):
if project.exists():
shutil.rmtree(project)
shutil.copytree(template, project)
base_ras = init_ras_project(base_project, '6.6', ras_object='new', load_results_summary=False)
proposed_ras = init_ras_project(proposed_project, '6.6', ras_object='new', load_results_summary=False)
polygon_ras = init_ras_project(polygon_project, '6.6', ras_object='new', load_results_summary=False)
shape_z_ras = init_ras_project(shape_z_project, '6.6', ras_object='new', load_results_summary=False)
plan_number = '18'
plan_row = base_ras.plan_df.loc[base_ras.plan_df['plan_number'] == plan_number].iloc[0]
geom_number = str(plan_row['Geom File'])
base_geom_hdf = Path(base_ras.geom_df.loc[base_ras.geom_df['geom_number'] == geom_number, 'hdf_path'].iloc[0])
proposed_geom_hdf = proposed_project / base_geom_hdf.name
polygon_geom_hdf = polygon_project / base_geom_hdf.name
shape_z_geom_hdf = shape_z_project / base_geom_hdf.name
terrain_value = proposed_ras.rasmap_df['terrain_hdf_path'].iloc[0]
proposed_terrain_hdf = Path(terrain_value[0] if isinstance(terrain_value, list) else terrain_value)
polygon_terrain_value = polygon_ras.rasmap_df['terrain_hdf_path'].iloc[0]
polygon_terrain_hdf = Path(
polygon_terrain_value[0]
if isinstance(polygon_terrain_value, list)
else polygon_terrain_value
)
shape_z_terrain_value = shape_z_ras.rasmap_df['terrain_hdf_path'].iloc[0]
shape_z_terrain_hdf = Path(
shape_z_terrain_value[0]
if isinstance(shape_z_terrain_value, list)
else shape_z_terrain_value
)
summary = pd.DataFrame([
{'condition': 'base', 'project': str(base_project), 'geom_hdf': str(base_geom_hdf)},
{'condition': 'high-ground proposed', 'project': str(proposed_project), 'terrain_hdf': str(proposed_terrain_hdf)},
{'condition': 'polygon multipoint proposed', 'project': str(polygon_project), 'terrain_hdf': str(polygon_terrain_hdf)},
{'condition': 'shape_z polygon example', 'project': str(shape_z_project), 'terrain_hdf': str(shape_z_terrain_hdf)},
])
summary
| condition | project | geom_hdf | terrain_hdf | |
|---|---|---|---|---|
| 0 | base | G:\GH\ras-commander-pr-docs\examples\example_p... | G:\GH\ras-commander-pr-docs\examples\example_p... | NaN |
| 1 | high-ground proposed | G:\GH\ras-commander-pr-docs\examples\example_p... | NaN | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 2 | polygon multipoint proposed | G:\GH\ras-commander-pr-docs\examples\example_p... | NaN | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 3 | shape_z polygon example | G:\GH\ras-commander-pr-docs\examples\example_p... | NaN | G:\GH\ras-commander-pr-docs\examples\example_p... |
Python
# Define a levee/road alignment across the BaldEagleCr 2D flow area
mesh_areas = HdfMesh.get_mesh_areas(base_geom_hdf)
target_area = mesh_areas.loc[mesh_areas['mesh_name'] == 'BaldEagleCr'].iloc[0]
minx, miny, maxx, maxy = target_area.geometry.bounds
line_y = miny + 0.55 * (maxy - miny)
profile_x = [minx + 0.15 * (maxx - minx), maxx - 0.15 * (maxx - minx)]
profile_y = [line_y, line_y]
high_ground_raise_ft = 6.0
high_ground_max_extent_ft = 600.0
RasTerrainMod.setup_gdal_bridge('6.6')
baseline_profile = RasTerrainMod.get_terrain_profile(
base_project / 'BaldEagleDamBrk.rasmap',
base_geom_hdf,
profile_x,
profile_y,
filter_tolerance=4.0,
)
profile_points = baseline_profile[['station', 'elevation']].to_numpy(dtype=float).copy()
profile_points[:, 1] += high_ground_raise_ft
mod_line = np.array([[profile_x[0], profile_y[0]], [profile_x[1], profile_y[1]]], dtype=float)
RasTerrainModWriter.add_high_ground_modification(
terrain_hdf_path=proposed_terrain_hdf,
rasmap_path=proposed_project / 'BaldEagleDamBrk.rasmap',
name='API High Ground Levee',
polyline_points=mod_line,
profile_points=profile_points,
top_width=120.0,
side_slope=4.0,
max_extent=high_ground_max_extent_ft,
elev_pt_tolerance=4.0,
)
mods = RasTerrainModWriter.list_modifications(proposed_terrain_hdf)
rasmap_text = (proposed_project / 'BaldEagleDamBrk.rasmap').read_text(encoding='utf-8')
print('Modification registered in .rasmap:', 'API High Ground Levee' in rasmap_text)
mods
Text Only
Modification registered in .rasmap: True
| name | type | subtype | priority | modification_type | modification_mode | width | left_slope | right_slope | max_extent | elev_pt_tolerance | profile_points | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | API High Ground Levee | Levee | Levee | 0 | 1 | take_higher | 120.0 | 4.0 | 4.0 | 600.0 | 4.0 | 260 |
Python
# Before/after terrain profile from the writer-side verification helper
profile_compare = RasTerrainModWriter.apply_modification_to_profile(
proposed_terrain_hdf,
'API High Ground Levee',
baseline_profile,
x_coords=profile_x,
y_coords=profile_y,
)
profile_stats = profile_compare['difference'].describe().to_frame('terrain_difference_ft')
display(profile_stats)
assert profile_compare['difference'].max() <= 8.0
assert profile_compare['difference'].max() >= 5.0
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(profile_compare['station'], profile_compare['existing_elevation'], label='Existing terrain')
ax.plot(profile_compare['station'], profile_compare['proposed_elevation'], label='High-ground surface')
ax.fill_between(
profile_compare['station'],
profile_compare['existing_elevation'],
profile_compare['proposed_elevation'],
where=profile_compare['difference'] > 0,
alpha=0.25,
)
ax.set_xlabel('Station (ft)')
ax.set_ylabel('Elevation (ft)')
ax.set_title(f'Before/after terrain profile ({high_ground_raise_ft:.0f} ft raise)')
ax.legend()
ax.grid(True, alpha=0.25)
profile_png = ARTIFACT_ROOT / 'terrain_profile_difference.png'
fig.tight_layout()
fig.savefig(profile_png, dpi=160)
profile_png
| terrain_difference_ft | |
|---|---|
| count | 2.600000e+02 |
| mean | 6.000000e+00 |
| std | 5.025337e-13 |
| min | 6.000000e+00 |
| 25% | 6.000000e+00 |
| 50% | 6.000000e+00 |
| 75% | 6.000000e+00 |
| max | 6.000000e+00 |
Text Only
WindowsPath('G:/GH/ras-commander-pr-docs/examples/example_projects/316_terrain_modifications/terrain_profile_difference.png')
Python
# Read the proposed terrain back through RasMapperLib and compare to writer expectation
readback_filter_tolerance = 0.01
readback_existing = RasTerrainMod.get_terrain_profile(
base_project / 'BaldEagleDamBrk.rasmap',
base_geom_hdf,
profile_x,
profile_y,
filter_tolerance=readback_filter_tolerance,
)
readback_proposed = RasTerrainMod.get_terrain_profile(
proposed_project / 'BaldEagleDamBrk.rasmap',
proposed_geom_hdf,
profile_x,
profile_y,
filter_tolerance=readback_filter_tolerance,
)
assert not readback_existing.empty, 'Base raster-read terrain profile is empty'
assert not readback_proposed.empty, 'Proposed raster-read terrain profile is empty'
writer_readback_compare = RasTerrainModWriter.apply_modification_to_profile(
proposed_terrain_hdf,
'API High Ground Levee',
readback_existing[['station', 'elevation']],
x_coords=profile_x,
y_coords=profile_y,
)
raster_proposed_elevation = np.interp(
writer_readback_compare['station'].to_numpy(dtype=float),
readback_proposed['station'].to_numpy(dtype=float),
readback_proposed['elevation'].to_numpy(dtype=float),
)
readback_validation = writer_readback_compare[[
'station',
'x',
'y',
'existing_elevation',
'proposed_elevation',
'difference',
]].copy()
readback_validation = readback_validation.rename(columns={
'proposed_elevation': 'writer_proposed_elevation',
'difference': 'writer_difference',
})
readback_validation['raster_read_elevation'] = raster_proposed_elevation
readback_validation['raster_writer_difference'] = (
readback_validation['raster_read_elevation']
- readback_validation['writer_proposed_elevation']
)
readback_validation['raster_change_ft'] = (
readback_validation['raster_read_elevation']
- readback_validation['existing_elevation']
)
readback_max_abs_diff = float(
readback_validation['raster_writer_difference'].abs().max()
)
readback_max_change = float(readback_validation['raster_change_ft'].max())
readback_min_change = float(readback_validation['raster_change_ft'].min())
readback_tolerance_ft = 0.01
readback_stats = readback_validation[
['writer_difference', 'raster_writer_difference', 'raster_change_ft']
].describe().T
display(readback_stats)
print(f'Max writer-vs-raster profile difference: {readback_max_abs_diff:.6f} ft')
print(f'High-ground raster terrain change range: {readback_min_change:.3f} to {readback_max_change:.3f} ft')
assert readback_max_abs_diff <= readback_tolerance_ft, (
f'Raster-read profile differs from writer-computed profile by '
f'{readback_max_abs_diff:.6f} ft, exceeding {readback_tolerance_ft:.3f} ft'
)
assert readback_max_change <= 10.0, (
f'High-ground raster terrain change exceeds 10 ft: {readback_max_change:.3f} ft'
)
assert readback_max_change >= 5.0, 'High-ground modification did not produce a measurable raise'
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(
readback_validation['station'],
readback_validation['existing_elevation'],
label='Existing raster',
color='0.35',
)
ax.plot(
readback_validation['station'],
readback_validation['writer_proposed_elevation'],
label='Writer-computed high ground',
color='tab:blue',
)
ax.plot(
readback_validation['station'],
readback_validation['raster_read_elevation'],
label='Proposed raster readback',
color='tab:orange',
linestyle='--',
)
ax.set_xlabel('Station (ft)')
ax.set_ylabel('Elevation (ft)')
ax.set_title('Writer-computed profile vs raster readback')
ax.legend()
ax.grid(True, alpha=0.25)
readback_profile_png = ARTIFACT_ROOT / 'terrain_profile_readback_validation.png'
fig.tight_layout()
fig.savefig(readback_profile_png, dpi=160)
readback_profile_png
| count | mean | std | min | 25% | 50% | 75% | max | |
|---|---|---|---|---|---|---|---|---|
| writer_difference | 6995.0 | 5.912489e+00 | 1.510301 | 0.633742 | 5.043287 | 6.0 | 6.728656 | 9.983736 |
| raster_writer_difference | 6995.0 | -1.294628e-08 | 0.000012 | -0.000059 | 0.000000 | 0.0 | 0.000000 | 0.000060 |
| raster_change_ft | 6995.0 | 5.912489e+00 | 1.510301 | 0.633746 | 5.043276 | 6.0 | 6.728654 | 9.983736 |
Text Only
Max writer-vs-raster profile difference: 0.000060 ft
High-ground raster terrain change range: 0.634 to 9.984 ft
WindowsPath('G:/GH/ras-commander-pr-docs/examples/example_projects/316_terrain_modifications/terrain_profile_readback_validation.png')
Native modification-aware terrain export (optional)¶
The writer checks above are analytical profile evidence. For a production GeoTIFF derivative, RasTerrain.export_rasmapper_terrain() calls RAS Mapper's native Export to Single Raster operation. This bounded example selects the original registered Terrain50 layer DataFrame-first, exports the same 2x grid with resampleVecMods disabled and enabled, and verifies the two-source inventory, managed receipt, pixel response, and unchanged controls.
Terrain50 intentionally combines 36.504512049933-ft and 20-ft sources. RAS Mapper owns reconciliation of those mixed resolutions at the explicit 40-ft output cell size with nearest-neighbor resampling. The 20-ft source is the authoritative grid anchor; no replacement terrain is created, and neither exported TIFF is registered back into the project.
Keep RUN_NATIVE_TERRAIN_EXPORT = False for routine documentation builds. To reproduce the committed native evidence on a machine with qualified HEC-RAS 6.6, set RAS_COMMANDER_RUN_NATIVE_TERRAIN_EXPORT=1 in the notebook kernel environment after reviewing the bounded extent. The written high-ground profile is 6 ft above sampled terrain, but the native TIN and mixed-source consolidation are evaluated on the output grid, so the qualified pixel-delta range is validated independently below rather than treating 6 ft as a raster ceiling.
Python
# Opt-in native RAS Mapper exports on the original mixed-resolution Terrain50
native_export_results = {}
native_export_receipts = {}
native_export_review_png = None
native_export_extent = (
2041660.41918676,
347030.9951331958,
2044060.41918676,
349430.9951331958,
)
native_export_root = EXAMPLE_RUN_ROOT / 'native_terrain_export'
if RUN_NATIVE_TERRAIN_EXPORT:
terrain_layers_before = RasMap.list_terrain_layers(
proposed_project,
ras_object=proposed_ras,
)
selected_terrain = terrain_layers_before.loc[
terrain_layers_before['name'].eq('Terrain50')
]
assert len(selected_terrain) == 1, (
'Expected exactly one case-sensitive Terrain50 registration'
)
terrain_name = str(selected_terrain.iloc[0]['name'])
display(selected_terrain[[
'name', 'filename', 'checked', 'type', 'resample_method', 'surface_on'
]])
for condition, rasterize_modifications in (
('modifications_off', False),
('modifications_on', True),
):
output_tif = native_export_root / f'terrain50_{condition}_2x.tif'
result = RasTerrain.export_rasmapper_terrain(
ras_project_path=proposed_project,
output_tif=output_tif,
terrain_name=terrain_name,
extent=native_export_extent,
downsample_factor=2,
rasterize_modifications=rasterize_modifications,
overwrite=OVERWRITE_NATIVE_TERRAIN_EXPORT,
timeout_seconds=240.0,
ras_object=proposed_ras,
)
assert isinstance(result, TerrainExportResult)
if not result:
raise RuntimeError(
f'{condition} export failed: {result.error}; '
f'receipt: {result.receipt_path}'
)
receipt = json.loads(result.receipt_path.read_text(encoding='utf-8'))
assert receipt['status'] == 'success'
assert receipt['result']['success'] is True
assert receipt['native_helper']['resample_to_one_rfi'] is True
assert receipt['native_helper']['resample_method'] == 'near'
assert receipt['result']['validation'] == result.validation
native_export_results[condition] = result
native_export_receipts[condition] = receipt
export_off = native_export_results['modifications_off']
export_on = native_export_results['modifications_on']
for result in (export_off, export_on):
assert result.native_cell_size == 20.0
assert result.output_cell_size == 40.0
assert result.downsample_factor == 2
assert (result.validation['columns'], result.validation['rows']) == (61, 61)
assert result.validation['band_count'] == 1
assert result.validation['data_type'] == 'Float32'
assert result.validation['crs_present'] is True
assert result.validation['sidecars'] == []
source_inventory = (
export_on.source_inventory.sort_values('index').reset_index(drop=True)
)
source_cells = source_inventory['cell_sizes'].map(lambda values: values[0])
np.testing.assert_allclose(
source_cells.to_numpy(dtype=float),
[36.504512049933, 20.0],
rtol=0.0,
atol=1e-9,
)
assert source_inventory['authoritative_grid'].tolist() == [False, True]
assert source_inventory['intersects_output'].all()
pd.testing.assert_frame_equal(
export_off.source_inventory.reset_index(drop=True),
export_on.source_inventory.reset_index(drop=True),
)
terrain_layers_after = RasMap.list_terrain_layers(
proposed_project,
ras_object=proposed_ras,
)
pd.testing.assert_frame_equal(terrain_layers_before, terrain_layers_after)
source_inventory_review = source_inventory.copy()
source_inventory_review['filename'] = (
source_inventory_review['filename'].map(lambda value: Path(value).name)
)
display(source_inventory_review)
display(pd.DataFrame([
{
'condition': condition,
'elapsed_seconds': result.elapsed_seconds,
'checksum': result.validation['checksum'],
'receipt': result.receipt_path.name,
}
for condition, result in native_export_results.items()
]))
else:
print('Native terrain export skipped; set RUN_NATIVE_TERRAIN_EXPORT = True to run it.')
print(f'Planned original terrain: Terrain50; bounded extent: {native_export_extent}')
print('Planned grid: 2x the authoritative 20-ft source = 40-ft cells.')
Text Only
name filename checked type resample_method \
0 Terrain50 .\Terrain\Terrain50.hdf True TerrainLayer near
surface_on
0 True
index filename priority columns rows \
0 0 Terrain50.baldeagledem.tif -1 8643 6902
1 1 Terrain50.dtm_20ft.tif -1 6705 4852
extent \
0 {'min_x': 1834327.1955903, 'min_y': 162918.805...
1 {'min_x': 1952354.5353037, 'min_y': 276821.135...
cell_sizes levels \
0 [36.504512049933, 73.009024099866, 146.0180481... 7
1 [20, 40, 80, 160, 320, 640] 6
intersects_output authoritative_grid
0 True False
1 True True
condition elapsed_seconds checksum \
0 modifications_off 3.348955 49168
1 modifications_on 1.842163 48517
receipt
0 terrain50_modifications_off_2x.tif.receipt.json
1 terrain50_modifications_on_2x.tif.receipt.json
Python
# Semantic off/on checks and review figure
native_export_stats = None
if native_export_results:
import rasterio
with rasterio.open(export_off.output_path) as baseline_dataset, rasterio.open(
export_on.output_path
) as modified_dataset:
baseline_terrain = baseline_dataset.read(1, masked=True)
modified_terrain = modified_dataset.read(1, masked=True)
assert baseline_terrain.shape == modified_terrain.shape == (61, 61)
assert baseline_dataset.transform == modified_dataset.transform
assert baseline_dataset.crs == modified_dataset.crs
x_centers = (
modified_dataset.transform.c
+ (np.arange(modified_dataset.width, dtype=float) + 0.5)
* modified_dataset.transform.a
)
y_centers = (
modified_dataset.transform.f
+ (np.arange(modified_dataset.height, dtype=float) + 0.5)
* modified_dataset.transform.e
)
x_grid, y_grid = np.meshgrid(x_centers, y_centers)
output_bounds = modified_dataset.bounds
valid_cells = (
~np.ma.getmaskarray(baseline_terrain)
& ~np.ma.getmaskarray(modified_terrain)
)
terrain_delta = np.ma.array(
modified_terrain - baseline_terrain,
mask=~valid_cells,
)
delta_values = terrain_delta.filled(0.0)
affected_cells = valid_cells & (delta_values > readback_tolerance_ft)
control_cells = valid_cells & (
np.abs(y_grid - line_y)
> high_ground_max_extent_ft + export_on.output_cell_size
)
valid_cell_count = int(valid_cells.sum())
affected_cell_count = int(affected_cells.sum())
affected_delta_min_ft = float(delta_values[affected_cells].min())
affected_delta_max_ft = float(delta_values[affected_cells].max())
minimum_delta_ft = float(delta_values[valid_cells].min())
control_cell_count = int(control_cells.sum())
control_max_abs_delta_ft = float(
np.abs(delta_values[control_cells]).max()
)
assert valid_cell_count == 3721
assert affected_cell_count == 264
np.testing.assert_allclose(
[affected_delta_min_ft, affected_delta_max_ft],
[0.15625, 9.625],
rtol=0.0,
atol=1e-6,
)
assert minimum_delta_ft >= 0.0, 'take_higher must not lower terrain cells'
assert control_cell_count == 1769
assert control_max_abs_delta_ft == 0.0
native_export_stats = pd.DataFrame([{
'valid_cells': valid_cell_count,
'raised_cells': affected_cell_count,
'minimum_raise_ft': affected_delta_min_ft,
'maximum_raise_ft': affected_delta_max_ft,
'control_cells': control_cell_count,
'control_max_abs_change_ft': control_max_abs_delta_ft,
}])
display(native_export_stats)
raster_extent = (
output_bounds.left,
output_bounds.right,
output_bounds.bottom,
output_bounds.top,
)
elevation_values = np.concatenate([
baseline_terrain.compressed(),
modified_terrain.compressed(),
])
elevation_min, elevation_max = np.nanpercentile(
elevation_values,
[2.0, 98.0],
)
fig, axes = plt.subplots(1, 3, figsize=(16, 5.5), constrained_layout=True)
for ax, values, title in (
(axes[0], baseline_terrain, 'Modifications disabled'),
(axes[1], modified_terrain, 'Native modifications enabled'),
):
image = ax.imshow(
values,
extent=raster_extent,
origin='upper',
cmap='terrain',
vmin=elevation_min,
vmax=elevation_max,
)
ax.plot(
[output_bounds.left, output_bounds.right],
[line_y, line_y],
color='white',
linewidth=1.0,
linestyle='--',
)
ax.set_title(title)
ax.set_xlabel('Project X (ft)')
ax.set_ylabel('Project Y (ft)')
ax.set_aspect('equal')
fig.colorbar(image, ax=ax, shrink=0.84, label='Elevation (ft)')
delta_image = axes[2].imshow(
terrain_delta,
extent=raster_extent,
origin='upper',
cmap='magma',
vmin=0.0,
vmax=affected_delta_max_ft,
)
axes[2].contour(
x_grid,
y_grid,
affected_cells.astype(int),
levels=[0.5],
colors=['cyan'],
linewidths=1.2,
)
control_rows, control_columns = np.where(control_cells)
control_stride = max(1, len(control_rows) // 16)
axes[2].scatter(
x_grid[control_rows[::control_stride], control_columns[::control_stride]],
y_grid[control_rows[::control_stride], control_columns[::control_stride]],
marker='x',
color='white',
s=28,
linewidths=1.0,
label='Unchanged control cells',
)
axes[2].plot(
[output_bounds.left, output_bounds.right],
[line_y, line_y],
color='cyan',
linewidth=1.0,
linestyle='--',
label='Written high ground',
)
axes[2].set_title('Modification-on minus modification-off')
axes[2].set_xlabel('Project X (ft)')
axes[2].set_ylabel('Project Y (ft)')
axes[2].set_aspect('equal')
axes[2].legend(loc='lower right', fontsize=8)
fig.colorbar(delta_image, ax=axes[2], shrink=0.84, label='Elevation change (ft)')
fig.suptitle(
'HEC-RAS 6.6 native mixed-source terrain export: 2x modification evidence\n'
f'{affected_cell_count} raised cells; {control_cell_count} controls unchanged',
fontsize=13,
)
native_export_review_png = (
native_export_root / 'terrain50_native_modification_evidence.png'
)
native_export_review_png.parent.mkdir(parents=True, exist_ok=True)
fig.savefig(native_export_review_png, dpi=180, bbox_inches='tight')
plt.show()
print(f'Native export review figure: {native_export_review_png.name}')
else:
print('Run the opt-in native export cell to generate semantic and visual evidence.')
Text Only
Native export review figure: terrain50_native_modification_evidence.png
valid_cells raised_cells minimum_raise_ft maximum_raise_ft \
0 3721 264 0.15625 9.625
control_cells control_max_abs_change_ft
0 1769 0.0
Python
# Polygon multipoint terrain modification writer validation
base_rasmap = base_project / 'BaldEagleDamBrk.rasmap'
polygon_name = 'API Multipoint Wetland Fill'
poly_cx = minx + 0.38 * (maxx - minx)
poly_cy = miny + 0.34 * (maxy - miny)
poly_w = 0.004 * (maxx - minx)
poly_h = 0.004 * (maxy - miny)
polygon_coords = np.array(
[
[poly_cx - 0.5 * poly_w, poly_cy - 0.5 * poly_h],
[poly_cx + 0.5 * poly_w, poly_cy - 0.5 * poly_h],
[poly_cx + 0.5 * poly_w, poly_cy + 0.5 * poly_h],
[poly_cx - 0.5 * poly_w, poly_cy + 0.5 * poly_h],
],
dtype=float,
)
control_xy = np.array(
[
[poly_cx - 0.18 * poly_w, poly_cy],
[poly_cx + 0.18 * poly_w, poly_cy],
],
dtype=float,
)
control_existing_profile = RasTerrainMod.get_terrain_profile(
base_rasmap,
base_geom_hdf,
control_xy[:, 0].tolist(),
control_xy[:, 1].tolist(),
filter_tolerance=0.0,
)
control_stations = np.r_[
0.0,
np.cumsum(np.hypot(np.diff(control_xy[:, 0]), np.diff(control_xy[:, 1]))),
]
control_existing_elevations = np.interp(
control_stations,
control_existing_profile['station'].to_numpy(dtype=float),
control_existing_profile['elevation'].to_numpy(dtype=float),
)
polygon_control_raises = np.array([5.0, 3.0], dtype=float)
polygon_control_elevations = control_existing_elevations + polygon_control_raises
polygon_control_points = [
{
'x': float(x),
'y': float(y),
'elevation': float(elevation),
'name': f'Wetland fill +{raise_ft:.0f} ft',
}
for (x, y), elevation, raise_ft in zip(
control_xy,
polygon_control_elevations,
polygon_control_raises,
)
]
RasTerrainModification.add_modification_polygon(
terrain_hdf_path=polygon_terrain_hdf,
name=polygon_name,
polygon_coords=polygon_coords,
elevation_method='boundary_from_terrain',
control_points=polygon_control_points,
rasmap_path=polygon_project / 'BaldEagleDamBrk.rasmap',
elev_pt_tolerance=15.0,
)
shape_z_name = 'API Shape Z Wetland Fill'
# Keep the shape_z writer example outside the active 2D mesh so it validates
# the alternate polygon writer path and HEC-RAS compute compatibility without
# introducing a second hydraulic perturbation.
shape_z_cx = minx + 0.75 * (maxx - minx)
shape_z_cy = miny + 0.30 * (maxy - miny)
shape_z_w = 0.004 * (maxx - minx)
shape_z_h = 0.004 * (maxy - miny)
shape_z_raise_ft = 4.0
shape_z_ring_xy = np.array(
[
[shape_z_cx - 0.5 * shape_z_w, shape_z_cy - 0.5 * shape_z_h],
[shape_z_cx + 0.5 * shape_z_w, shape_z_cy - 0.5 * shape_z_h],
[shape_z_cx + 0.5 * shape_z_w, shape_z_cy + 0.5 * shape_z_h],
[shape_z_cx - 0.5 * shape_z_w, shape_z_cy + 0.5 * shape_z_h],
],
dtype=float,
)
shape_z_closed_xy = np.vstack([shape_z_ring_xy, shape_z_ring_xy[0]])
shape_z_ring_profile = RasTerrainMod.get_terrain_profile(
base_rasmap,
base_geom_hdf,
shape_z_closed_xy[:, 0].tolist(),
shape_z_closed_xy[:, 1].tolist(),
filter_tolerance=0.0,
)
shape_z_ring_stations = np.r_[
0.0,
np.cumsum(
np.hypot(
np.diff(shape_z_closed_xy[:, 0]),
np.diff(shape_z_closed_xy[:, 1]),
)
),
]
shape_z_existing_ring_elevations = np.interp(
shape_z_ring_stations[:-1],
shape_z_ring_profile['station'].to_numpy(dtype=float),
shape_z_ring_profile['elevation'].to_numpy(dtype=float),
)
shape_z_polygon_coords = np.column_stack(
[shape_z_ring_xy, shape_z_existing_ring_elevations + shape_z_raise_ft]
)
shape_z_center_elevation = float(np.mean(shape_z_existing_ring_elevations))
shape_z_control_points = [
{
'x': float(shape_z_cx),
'y': float(shape_z_cy),
'elevation': shape_z_center_elevation + shape_z_raise_ft,
'name': f'shape_z center +{shape_z_raise_ft:.0f} ft',
}
]
RasTerrainModification.add_modification_polygon(
terrain_hdf_path=shape_z_terrain_hdf,
name=shape_z_name,
polygon_coords=shape_z_polygon_coords,
elevation_method='shape_z',
control_points=shape_z_control_points,
rasmap_path=shape_z_project / 'BaldEagleDamBrk.rasmap',
elev_pt_tolerance=10.0,
)
with h5py.File(polygon_terrain_hdf, 'r') as hdf:
mod = hdf[f'Modifications/{polygon_name}']
boundary_xyz = mod['Boundary Points'][:]
control_xyz = mod['Control Points/Elevation Points'][:]
profile_values = mod['Profile Values'][:]
polygon_attrs = mod['Attributes'][0]
assert np.all(np.isfinite(boundary_xyz)), 'Boundary elevations must be finite'
np.testing.assert_allclose(control_xyz[:, 2], polygon_control_elevations, atol=1e-6)
assert int(polygon_attrs['Generate Boundary Elevations']) == 1
assert int(polygon_attrs['Use ShapeFile Z Elevations']) == 0
root = ET.parse(polygon_project / 'BaldEagleDamBrk.rasmap').getroot()
polygon_layers = [
layer for layer in root.findall('.//Layer')
if layer.get('Name') == polygon_name
]
assert len(polygon_layers) == 1
assert polygon_layers[0].get('Type') == 'PolygonElevationModificationLayer'
with h5py.File(shape_z_terrain_hdf, 'r') as hdf:
shape_z_mod = hdf[f'Modifications/{shape_z_name}']
shape_z_boundary_xyz = shape_z_mod['Boundary Points'][:]
shape_z_profile_values = shape_z_mod['Profile Values'][:]
shape_z_attrs = shape_z_mod['Attributes'][0]
shape_z_expected = np.r_[shape_z_polygon_coords[:, 2], shape_z_polygon_coords[0, 2]]
np.testing.assert_allclose(shape_z_boundary_xyz[:, 2], shape_z_expected, atol=1e-6)
np.testing.assert_allclose(shape_z_profile_values[:, 1], shape_z_expected, atol=1e-6)
assert int(shape_z_attrs['Generate Boundary Elevations']) == 0
assert int(shape_z_attrs['Use ShapeFile Z Elevations']) == 1
shape_z_root = ET.parse(shape_z_project / 'BaldEagleDamBrk.rasmap').getroot()
shape_z_layers = [
layer for layer in shape_z_root.findall('.//Layer')
if layer.get('Name') == shape_z_name
]
assert len(shape_z_layers) == 1
assert shape_z_layers[0].get('Type') == 'PolygonElevationModificationLayer'
polygon_summary = pd.DataFrame({
'modification': [polygon_name, polygon_name, polygon_name, shape_z_name, shape_z_name],
'dataset': ['boundary vertices', 'interior controls', 'profile values', 'boundary vertices', 'profile values'],
'count': [
len(boundary_xyz),
len(control_xyz),
len(profile_values),
len(shape_z_boundary_xyz),
len(shape_z_profile_values),
],
'min_elevation_ft': [
float(boundary_xyz[:, 2].min()),
float(control_xyz[:, 2].min()),
float(profile_values[:, 1].min()),
float(shape_z_boundary_xyz[:, 2].min()),
float(shape_z_profile_values[:, 1].min()),
],
'max_elevation_ft': [
float(boundary_xyz[:, 2].max()),
float(control_xyz[:, 2].max()),
float(profile_values[:, 1].max()),
float(shape_z_boundary_xyz[:, 2].max()),
float(shape_z_profile_values[:, 1].max()),
],
})
display(polygon_summary)
fig, ax = plt.subplots(figsize=(8, 6))
mesh_areas.boundary.plot(ax=ax, linewidth=0.8, color='0.35')
closed_polygon = np.vstack([polygon_coords, polygon_coords[0]])
ax.plot(closed_polygon[:, 0], closed_polygon[:, 1], color='tab:blue', linewidth=2.0)
shape_z_closed_polygon = np.vstack([shape_z_polygon_coords[:, :2], shape_z_polygon_coords[0, :2]])
ax.plot(shape_z_closed_polygon[:, 0], shape_z_closed_polygon[:, 1], color='tab:green', linewidth=2.0)
sc_boundary = ax.scatter(
boundary_xyz[:, 0],
boundary_xyz[:, 1],
c=boundary_xyz[:, 2],
cmap='terrain',
marker='s',
s=55,
label='Sampled boundary',
)
sc_control = ax.scatter(
control_xyz[:, 0],
control_xyz[:, 1],
c=control_xyz[:, 2],
cmap='viridis',
marker='o',
edgecolor='black',
s=90,
label='Interior controls',
)
ax.scatter(
shape_z_boundary_xyz[:, 0],
shape_z_boundary_xyz[:, 1],
color='tab:green',
marker='^',
s=70,
label='shape_z ring vertices',
)
for point, raise_ft in zip(control_xyz, polygon_control_raises):
ax.annotate(f'+{raise_ft:.0f} ft', xy=(point[0], point[1]), xytext=(6, 6), textcoords='offset points')
ax.annotate(f'shape_z +{shape_z_raise_ft:.0f} ft ring', xy=(shape_z_cx, shape_z_cy), xytext=(6, 6), textcoords='offset points')
ax.set_aspect('equal')
ax.set_title('Polygon terrain modification inputs')
ax.set_xlabel('X (ft)')
ax.set_ylabel('Y (ft)')
ax.legend(loc='upper right')
fig.colorbar(sc_boundary, ax=ax, label='Boundary elevation (ft)')
polygon_points_png = ARTIFACT_ROOT / 'polygon_multipoint_control_points.png'
fig.tight_layout()
fig.savefig(polygon_points_png, dpi=160)
polygon_points_png
Text Only
2026-07-07 18:32:11 - rasterio._env - WARNING - CPLE_AppDefined in PROJ: proj_create_from_database: C:\Program Files (x86)\HEC\HEC-RAS\6.6\GDAL\common\data\proj.db lacks DATABASE.LAYOUT.VERSION.MAJOR / DATABASE.LAYOUT.VERSION.MINOR metadata. It comes from another PROJ installation.
2026-07-07 18:32:11 - rasterio._env - WARNING - CPLE_AppDefined in PROJ: proj_create_from_name: C:\Program Files (x86)\HEC\HEC-RAS\6.6\GDAL\common\data\proj.db lacks DATABASE.LAYOUT.VERSION.MAJOR / DATABASE.LAYOUT.VERSION.MINOR metadata. It comes from another PROJ installation.
| modification | dataset | count | min_elevation_ft | max_elevation_ft | |
|---|---|---|---|---|---|
| 0 | API Multipoint Wetland Fill | boundary vertices | 5 | 564.843750 | 573.062500 |
| 1 | API Multipoint Wetland Fill | interior controls | 2 | 568.876953 | 570.660278 |
| 2 | API Multipoint Wetland Fill | profile values | 5 | 564.843750 | 573.062500 |
| 3 | API Shape Z Wetland Fill | boundary vertices | 5 | 763.179932 | 765.740478 |
| 4 | API Shape Z Wetland Fill | profile values | 5 | 763.179932 | 765.740479 |
Text Only
WindowsPath('G:/GH/ras-commander-pr-docs/examples/example_projects/316_terrain_modifications/polygon_multipoint_control_points.png')
Python
# Read the polygon modification back through RasMapperLib and compare profile response
def interpolate_tin_surface(points_xyz, sample_x, sample_y):
surface_points = (
pd.DataFrame(points_xyz, columns=['x', 'y', 'elevation'])
.drop_duplicates(subset=['x', 'y'], keep='first')
)
triangulation = mtri.Triangulation(
surface_points['x'].to_numpy(dtype=float),
surface_points['y'].to_numpy(dtype=float),
)
interpolator = mtri.LinearTriInterpolator(
triangulation,
surface_points['elevation'].to_numpy(dtype=float),
)
values = interpolator(np.asarray(sample_x, dtype=float), np.asarray(sample_y, dtype=float))
return np.asarray(np.ma.filled(values, np.nan), dtype=float)
polygon_profile_x = np.linspace(poly_cx - 0.85 * poly_w, poly_cx + 0.85 * poly_w, 121).tolist()
polygon_profile_y = [poly_cy] * len(polygon_profile_x)
polygon_profile_tolerance = 0.0
polygon_existing_profile = RasTerrainMod.get_terrain_profile(
base_project / 'BaldEagleDamBrk.rasmap',
base_geom_hdf,
polygon_profile_x,
polygon_profile_y,
filter_tolerance=polygon_profile_tolerance,
)
polygon_readback_profile = RasTerrainMod.get_terrain_profile(
polygon_project / 'BaldEagleDamBrk.rasmap',
polygon_geom_hdf,
polygon_profile_x,
polygon_profile_y,
filter_tolerance=polygon_profile_tolerance,
)
assert not polygon_existing_profile.empty, 'Base polygon profile is empty'
assert not polygon_readback_profile.empty, 'Polygon readback profile is empty'
polygon_stations = polygon_readback_profile['station'].to_numpy(dtype=float)
polygon_existing_interp = np.interp(
polygon_stations,
polygon_existing_profile['station'].to_numpy(dtype=float),
polygon_existing_profile['elevation'].to_numpy(dtype=float),
)
profile_length = float(np.hypot(
polygon_profile_x[-1] - polygon_profile_x[0],
polygon_profile_y[-1] - polygon_profile_y[0],
))
polygon_sample_x = (
polygon_profile_x[0]
+ (polygon_stations / profile_length)
* (polygon_profile_x[-1] - polygon_profile_x[0])
)
polygon_sample_y = (
polygon_profile_y[0]
+ (polygon_stations / profile_length)
* (polygon_profile_y[-1] - polygon_profile_y[0])
)
inside_polygon = MplPath(boundary_xyz[:-1, :2]).contains_points(
np.column_stack([polygon_sample_x, polygon_sample_y]),
radius=1e-6,
)
inside_polygon_count = int(inside_polygon.sum())
assert inside_polygon_count > 3, (
f'Expected more than 3 profile points inside polygon; got {inside_polygon_count}'
)
expected_polygon_surface = interpolate_tin_surface(
np.vstack([boundary_xyz, control_xyz]),
polygon_sample_x,
polygon_sample_y,
)
polygon_validation = polygon_readback_profile[['station', 'elevation']].copy()
polygon_validation = polygon_validation.rename(columns={'elevation': 'raster_read_elevation'})
polygon_validation['x'] = polygon_sample_x
polygon_validation['y'] = polygon_sample_y
polygon_validation['existing_elevation'] = polygon_existing_interp
polygon_validation['expected_tin_elevation'] = expected_polygon_surface
polygon_validation['inside_polygon'] = inside_polygon
polygon_validation['readback_difference'] = (
polygon_validation['raster_read_elevation']
- polygon_validation['existing_elevation']
)
polygon_validation['raster_minus_expected_tin_ft'] = (
polygon_validation['raster_read_elevation']
- polygon_validation['expected_tin_elevation']
)
inside_validation = polygon_validation[
polygon_validation['inside_polygon']
& np.isfinite(polygon_validation['expected_tin_elevation'])
].copy()
assert not inside_validation.empty, 'Polygon raster profile had no finite points inside the footprint'
polygon_max_abs_change = float(inside_validation['readback_difference'].abs().max())
polygon_max_abs_tin_diff = float(inside_validation['raster_minus_expected_tin_ft'].abs().max())
assert polygon_max_abs_change <= 10.0, (
f'Polygon terrain change exceeds 10 ft: {polygon_max_abs_change:.3f} ft'
)
assert polygon_max_abs_change >= 2.0, 'Polygon modification did not produce a measurable terrain change'
polygon_raster_tin_tolerance_ft = 1.0
assert polygon_max_abs_tin_diff <= polygon_raster_tin_tolerance_ft, (
f'Polygon raster-read profile differs from independent TIN expectation by '
f'{polygon_max_abs_tin_diff:.6f} ft'
)
control_checks = []
for control_point in control_xyz:
distance = np.hypot(
polygon_validation['x'].to_numpy(dtype=float) - control_point[0],
polygon_validation['y'].to_numpy(dtype=float) - control_point[1],
)
nearest_idx = int(np.argmin(distance))
row = polygon_validation.iloc[nearest_idx]
control_checks.append({
'target_x': float(control_point[0]),
'target_y': float(control_point[1]),
'target_elevation_ft': float(control_point[2]),
'readback_station_ft': float(row['station']),
'readback_elevation_ft': float(row['raster_read_elevation']),
'readback_minus_target_ft': float(row['raster_read_elevation'] - control_point[2]),
'profile_distance_ft': float(distance[nearest_idx]),
})
control_check_df = pd.DataFrame(control_checks)
control_max_abs_diff = float(control_check_df['readback_minus_target_ft'].abs().max())
assert control_max_abs_diff <= 1.5, (
f'Nearest raster profile point differs from control elevations by '
f'{control_max_abs_diff:.6f} ft'
)
polygon_profile_stats = inside_validation[[
'readback_difference',
'raster_minus_expected_tin_ft',
]].describe().T
polygon_profile_stats['inside_profile_points'] = inside_polygon_count
polygon_profile_stats['max_abs_ft'] = [
polygon_max_abs_change,
polygon_max_abs_tin_diff,
]
display(control_check_df)
display(polygon_profile_stats)
print(f'Polygon profile points inside footprint: {inside_polygon_count}')
print(f'Control-point nearest raster max abs diff: {control_max_abs_diff:.6f} ft')
print(f'Polygon raster-vs-independent-TIN max abs diff: {polygon_max_abs_tin_diff:.6f} ft')
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(
polygon_validation['station'],
polygon_validation['existing_elevation'],
label='Existing raster',
color='0.35',
)
ax.plot(
polygon_validation['station'],
polygon_validation['expected_tin_elevation'],
label='Independent TIN expectation',
color='tab:blue',
)
ax.plot(
polygon_validation['station'],
polygon_validation['raster_read_elevation'],
label='Polygon raster readback',
color='tab:orange',
linestyle='--',
)
control_station = [
float(np.hypot(point[0] - polygon_profile_x[0], point[1] - polygon_profile_y[0]))
for point in control_xyz
]
ax.scatter(
control_station,
control_xyz[:, 2],
color='black',
zorder=5,
label='Written control elevations',
)
inside_rows = polygon_validation[polygon_validation['inside_polygon']]
if not inside_rows.empty:
ax.axvspan(
inside_rows['station'].min(),
inside_rows['station'].max(),
color='tab:blue',
alpha=0.08,
label='Polygon footprint',
)
ax.set_xlabel('Station (ft)')
ax.set_ylabel('Elevation (ft)')
ax.set_title('Polygon multipoint raster readback vs independent TIN')
ax.legend()
ax.grid(True, alpha=0.25)
polygon_profile_png = ARTIFACT_ROOT / 'polygon_profile_readback_validation.png'
fig.tight_layout()
fig.savefig(polygon_profile_png, dpi=160)
polygon_profile_png
| target_x | target_y | target_elevation_ft | readback_station_ft | readback_elevation_ft | readback_minus_target_ft | profile_distance_ft | |
|---|---|---|---|---|---|---|---|
| 0 | 2.032737e+06 | 337601.921003 | 570.660278 | 224.743301 | 570.660278 | 0.0 | 0.000006 |
| 1 | 2.032858e+06 | 337601.921003 | 568.876953 | 345.500885 | 568.876953 | 0.0 | 0.000002 |
| count | mean | std | min | 25% | 50% | 75% | max | inside_profile_points | max_abs_ft | |
|---|---|---|---|---|---|---|---|---|---|---|
| readback_difference | 157.0 | 3.568489 | 1.108918 | 0.743785 | 3.072388 | 3.403687e+00 | 4.497009 | 5.009216 | 157 | 5.009216 |
| raster_minus_expected_tin_ft | 157.0 | 0.000002 | 0.000017 | -0.000029 | -0.000014 | 1.113549e-07 | 0.000015 | 0.000031 | 157 | 0.000031 |
Text Only
Polygon profile points inside footprint: 157
Control-point nearest raster max abs diff: 0.000000 ft
Polygon raster-vs-independent-TIN max abs diff: 0.000031 ft
WindowsPath('G:/GH/ras-commander-pr-docs/examples/example_projects/316_terrain_modifications/polygon_profile_readback_validation.png')
Python
# Verify shape_z polygon writer HDF profile readback against ring Z values
shape_z_existing_closed = np.r_[
shape_z_existing_ring_elevations,
shape_z_existing_ring_elevations[0],
]
shape_z_readback = pd.DataFrame({
'station': shape_z_profile_values[:, 0].astype(float),
'existing_elevation': shape_z_existing_closed.astype(float),
'readback_elevation': shape_z_profile_values[:, 1].astype(float),
'expected_elevation': shape_z_expected.astype(float),
})
shape_z_readback['terrain_change_ft'] = (
shape_z_readback['readback_elevation']
- shape_z_readback['existing_elevation']
)
shape_z_readback['readback_minus_shape_z_ft'] = (
shape_z_readback['readback_elevation']
- shape_z_readback['expected_elevation']
)
shape_z_max_abs_diff = float(shape_z_readback['readback_minus_shape_z_ft'].abs().max())
shape_z_max_abs_change = float(shape_z_readback['terrain_change_ft'].abs().max())
assert shape_z_max_abs_diff <= 1e-6, (
f'shape_z HDF profile readback differs from ring Z by {shape_z_max_abs_diff:.6f} ft'
)
assert shape_z_max_abs_change <= 10.0, (
f'shape_z terrain change exceeds 10 ft: {shape_z_max_abs_change:.3f} ft'
)
assert shape_z_readback['terrain_change_ft'].min() >= 3.0
shape_z_stats = shape_z_readback[[
'terrain_change_ft',
'readback_minus_shape_z_ft',
]].describe().T
shape_z_stats['profile_points'] = len(shape_z_readback)
shape_z_stats['max_abs_ft'] = [
shape_z_max_abs_change,
shape_z_max_abs_diff,
]
display(shape_z_readback)
display(shape_z_stats)
print(f'shape_z HDF readback max abs diff: {shape_z_max_abs_diff:.6f} ft')
print(f'shape_z terrain change max abs diff: {shape_z_max_abs_change:.6f} ft')
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(
shape_z_readback['station'],
shape_z_readback['existing_elevation'],
marker='o',
label='Existing ring elevation',
color='0.35',
)
ax.plot(
shape_z_readback['station'],
shape_z_readback['readback_elevation'],
marker='o',
label='Written shape_z profile',
color='tab:green',
)
ax.set_xlabel('Station (ft)')
ax.set_ylabel('Elevation (ft)')
ax.set_title('shape_z polygon writer HDF readback')
ax.legend()
ax.grid(True, alpha=0.25)
shape_z_profile_png = ARTIFACT_ROOT / 'shape_z_polygon_readback_validation.png'
fig.tight_layout()
fig.savefig(shape_z_profile_png, dpi=160)
shape_z_profile_png
| station | existing_elevation | readback_elevation | expected_elevation | terrain_change_ft | readback_minus_shape_z_ft | |
|---|---|---|---|---|---|---|
| 0 | 0.000000 | 759.211426 | 763.211426 | 763.211426 | 4.0 | 0.000000e+00 |
| 1 | 335.437744 | 759.644409 | 763.644409 | 763.644409 | 4.0 | -2.643937e-07 |
| 2 | 537.896301 | 761.740478 | 765.740479 | 765.740478 | 4.0 | 1.785431e-07 |
| 3 | 873.334045 | 759.179932 | 763.179932 | 763.179932 | 4.0 | 0.000000e+00 |
| 4 | 1075.792603 | 759.211426 | 763.211426 | 763.211426 | 4.0 | 0.000000e+00 |
| count | mean | std | min | 25% | 50% | 75% | max | profile_points | max_abs_ft | |
|---|---|---|---|---|---|---|---|---|---|---|
| terrain_change_ft | 5.0 | 4.000000e+00 | 1.583569e-07 | 4.000000e+00 | 4.0 | 4.0 | 4.0 | 4.000000e+00 | 5 | 4.000000e+00 |
| readback_minus_shape_z_ft | 5.0 | -1.717012e-08 | 1.583569e-07 | -2.643937e-07 | 0.0 | 0.0 | 0.0 | 1.785431e-07 | 5 | 2.643937e-07 |
Text Only
shape_z HDF readback max abs diff: 0.000000 ft
shape_z terrain change max abs diff: 4.000000 ft
WindowsPath('G:/GH/ras-commander-pr-docs/examples/example_projects/316_terrain_modifications/shape_z_polygon_readback_validation.png')
Python
# Terrain difference map at geometry cell-center points
cell_points = HdfMesh.get_mesh_cell_points(base_geom_hdf)
cell_elev = pd.DataFrame({
'mesh_name': cell_points['mesh_name'],
'cell_id': cell_points['cell_id'],
'x': cell_points.geometry.x,
'y': cell_points.geometry.y,
})
cell_probe = RasTerrainModWriter.sample_modification_surface(
proposed_terrain_hdf,
'API High Ground Levee',
cell_elev[['x', 'y']],
)
cell_existing = np.interp(
cell_probe['line_station'],
baseline_profile['station'],
baseline_profile['elevation'],
)
cell_mod = RasTerrainModWriter.sample_modification_surface(
proposed_terrain_hdf,
'API High Ground Levee',
cell_elev[['x', 'y']],
existing_elevations=cell_existing,
)
terrain_map = pd.concat(
[cell_elev[['mesh_name', 'cell_id']], cell_mod],
axis=1,
)
terrain_map_stats = terrain_map['difference'].describe().to_frame('cell_terrain_difference_ft')
display(terrain_map_stats)
fig, ax = plt.subplots(figsize=(8, 6))
plot_df = terrain_map[terrain_map['difference'] > 0.01]
mesh_areas.boundary.plot(ax=ax, linewidth=0.8, color='0.35')
sc = ax.scatter(plot_df['x'], plot_df['y'], c=plot_df['difference'], s=8, cmap='viridis')
ax.plot(profile_x, profile_y, color='black', linewidth=1.5)
ax.set_aspect('equal')
ax.set_title('High-ground terrain difference at mesh cells')
ax.set_xlabel('X (ft)')
ax.set_ylabel('Y (ft)')
fig.colorbar(sc, ax=ax, label='Terrain raise (ft)')
terrain_png = ARTIFACT_ROOT / 'terrain_difference_map.png'
fig.tight_layout()
fig.savefig(terrain_png, dpi=160)
terrain_png
| cell_terrain_difference_ft | |
|---|---|
| count | 5198.000000 |
| mean | 0.049019 |
| std | 0.538554 |
| min | 0.000000 |
| 25% | 0.000000 |
| 50% | 0.000000 |
| 75% | 0.000000 |
| max | 6.000000 |
Text Only
WindowsPath('G:/GH/ras-commander-pr-docs/examples/example_projects/316_terrain_modifications/terrain_difference_map.png')
Python
# Run HEC-RAS base, high-ground, multipoint polygon, and shape_z polygon plans through ras-commander
def force_clean_plan_inputs(project_path):
for stale in project_path.glob('*.g11.hdf'):
stale.unlink()
for stale in project_path.glob('*.c18'):
stale.unlink()
compute_rows = []
for label, project_path, ras_obj in [
('base', base_project, base_ras),
('high-ground proposed', proposed_project, proposed_ras),
('polygon proposed', polygon_project, polygon_ras),
('shape_z polygon proposed', shape_z_project, shape_z_ras),
]:
force_clean_plan_inputs(project_path)
result = RasCmdr.compute_plan(
plan_number,
ras_object=ras_obj,
force_geompre=True,
force_rerun=True,
num_cores=2,
verify=True,
)
compute_rows.append({
'condition': label,
'success': bool(result),
'plan_hdf': str(project_path / 'BaldEagleDamBrk.p18.hdf'),
})
compute_summary = pd.DataFrame(compute_rows)
assert compute_summary['success'].all(), 'One or more HEC-RAS plan computations failed'
compute_summary
| condition | success | plan_hdf | |
|---|---|---|---|
| 0 | base | True | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 1 | high-ground proposed | True | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 2 | polygon proposed | True | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 3 | shape_z polygon proposed | True | G:\GH\ras-commander-pr-docs\examples\example_p... |
Interpreting Benefit/Rise Areas¶
This is a benefits evaluation demonstration, not a conclusion that the evaluated terrain change has net benefits. HdfBenefitAreas surfaces calculated spatial differences between the existing and proposed plan HDF results as benefit and adverse-impact polygons for professional review. The mapped downstream rise area in this high-ground example is an impact signal from the model results; it may reflect that the terrain change was not paired with detention/storage, or that the project footprint needs further refinement or mitigation review. Treat the polygons as an engineering screening layer that identifies where to investigate, not as a standalone design conclusion.
Python
# Compare maximum water surface elevations from result HDF files
base_plan_hdf = base_project / 'BaldEagleDamBrk.p18.hdf'
proposed_plan_hdf = proposed_project / 'BaldEagleDamBrk.p18.hdf'
polygon_plan_hdf = polygon_project / 'BaldEagleDamBrk.p18.hdf'
shape_z_plan_hdf = shape_z_project / 'BaldEagleDamBrk.p18.hdf'
base_wse = HdfResultsMesh.get_mesh_max_ws(base_plan_hdf)
# HdfBenefitAreas packages the common existing/proposed comparison into
# reusable benefit/rise polygons plus matched point differences.
benefit_results = HdfBenefitAreas.identify_benefit_areas(
existing_hdf_path=base_plan_hdf,
proposed_hdf_path=proposed_plan_hdf,
min_delta=0.01,
adjacency_method='polygon_edges',
)
def benefit_area_summary(results):
rows = []
for label in ['benefit_polygons', 'rise_polygons']:
gdf = results[label]
rows.append({
'result': label,
'features': len(gdf),
'cells': int(gdf['cell_count'].sum()) if not gdf.empty else 0,
'area_acres': float(gdf['area_acres'].sum()) if not gdf.empty else 0.0,
'min_wse_difference_ft': np.nan,
'max_wse_difference_ft': np.nan,
})
diff = results['difference_points']
rows.append({
'result': 'difference_points',
'features': len(diff),
'cells': int((diff['wse_difference'].abs() >= 0.01).sum()),
'area_acres': np.nan,
'min_wse_difference_ft': float(diff['wse_difference'].min()),
'max_wse_difference_ft': float(diff['wse_difference'].max()),
})
return pd.DataFrame(rows)
benefit_summary = benefit_area_summary(benefit_results)
display(benefit_summary)
# Use the API difference points for the high-ground map/summary; keep the manual
# helper below for the other terrain variants.
high_ground_wse = benefit_results['difference_points'].rename(columns={'geometry': 'geometry_base'})
high_ground_wse = high_ground_wse.set_geometry('geometry_base')
def wse_difference_from_base(plan_hdf, suffix):
scenario_wse = HdfResultsMesh.get_mesh_max_ws(plan_hdf)
result = base_wse.merge(
scenario_wse,
on=['mesh_name', 'cell_id'],
suffixes=('_base', f'_{suffix}'),
)
result['wse_difference'] = (
result[f'maximum_water_surface_{suffix}']
- result['maximum_water_surface_base']
)
return result
polygon_wse = wse_difference_from_base(polygon_plan_hdf, 'polygon')
shape_z_wse = wse_difference_from_base(shape_z_plan_hdf, 'shape_z')
wse_scenarios = [
('high-ground proposed', high_ground_wse),
('polygon proposed', polygon_wse),
('shape_z polygon proposed', shape_z_wse),
]
wse_summary = pd.DataFrame([
{
'scenario': scenario,
'min_wse_difference_ft': float(wse_df['wse_difference'].min()),
'max_wse_difference_ft': float(wse_df['wse_difference'].max()),
'mean_wse_difference_ft': float(wse_df['wse_difference'].mean()),
'max_abs_wse_difference_ft': float(wse_df['wse_difference'].abs().max()),
'cells_abs_gt_0_01_ft': int((wse_df['wse_difference'].abs() > 0.01).sum()),
}
for scenario, wse_df in wse_scenarios
])
display(wse_summary)
for scenario, wse_df in wse_scenarios:
print(
f'{scenario} cells with |WSE difference| > 0.01 ft:',
int((wse_df['wse_difference'].abs() > 0.01).sum()),
)
def plot_wse_difference(wse_df, title, output_path, overlay='line'):
fig, ax = plt.subplots(figsize=(8, 6))
wse_plot = wse_df[np.isfinite(wse_df['wse_difference'])].copy().set_geometry('geometry_base')
wse_plot.plot(
ax=ax,
column='wse_difference',
cmap='coolwarm',
markersize=7,
legend=True,
legend_kwds={'label': 'Max WSE difference (ft)'},
)
mesh_areas.boundary.plot(ax=ax, linewidth=0.8, color='0.25')
if overlay == 'line':
ax.plot(profile_x, profile_y, color='black', linewidth=1.5)
elif overlay == 'polygon':
ax.plot(closed_polygon[:, 0], closed_polygon[:, 1], color='black', linewidth=1.5)
elif overlay == 'shape_z_polygon':
ax.plot(shape_z_closed_polygon[:, 0], shape_z_closed_polygon[:, 1], color='black', linewidth=1.5)
ax.set_aspect('equal')
ax.set_title(title)
ax.set_xlabel('X (ft)')
ax.set_ylabel('Y (ft)')
fig.tight_layout()
fig.savefig(output_path, dpi=160)
return output_path
def plot_benefit_rise_areas(results, output_path):
fig, ax = plt.subplots(figsize=(8, 6))
mesh_areas.boundary.plot(ax=ax, linewidth=0.8, color='0.25')
handles = []
benefit_gdf = results['benefit_polygons']
rise_gdf = results['rise_polygons']
if not benefit_gdf.empty:
benefit_gdf.plot(ax=ax, facecolor='#2ca25f', edgecolor='#006d2c', alpha=0.45, linewidth=0.7)
handles.append(plt.Line2D([0], [0], color='#006d2c', linewidth=6, alpha=0.45, label='Benefit: WSE reduced'))
if not rise_gdf.empty:
rise_gdf.plot(ax=ax, facecolor='#de2d26', edgecolor='#a50f15', alpha=0.45, linewidth=0.7)
handles.append(plt.Line2D([0], [0], color='#a50f15', linewidth=6, alpha=0.45, label='Rise: WSE increased'))
ax.plot(profile_x, profile_y, color='black', linewidth=1.5)
if handles:
ax.legend(handles=handles, loc='upper right')
ax.set_aspect('equal')
ax.set_title('High-ground benefit/rise areas from HdfBenefitAreas')
ax.set_xlabel('X (ft)')
ax.set_ylabel('Y (ft)')
fig.tight_layout()
fig.savefig(output_path, dpi=160)
return output_path
wse_png = ARTIFACT_ROOT / 'wse_difference_map.png'
polygon_wse_png = ARTIFACT_ROOT / 'polygon_wse_difference_map.png'
shape_z_wse_png = ARTIFACT_ROOT / 'shape_z_wse_difference_map.png'
benefit_rise_png = ARTIFACT_ROOT / 'benefit_rise_area_map.png'
plot_wse_difference(high_ground_wse, 'HEC-RAS high-ground maximum WSE impact', wse_png, overlay='line')
plot_wse_difference(polygon_wse, 'HEC-RAS polygon maximum WSE impact', polygon_wse_png, overlay='polygon')
plot_wse_difference(shape_z_wse, 'HEC-RAS shape_z polygon maximum WSE impact', shape_z_wse_png, overlay='shape_z_polygon')
plot_benefit_rise_areas(benefit_results, benefit_rise_png)
pd.DataFrame({
'artifact': [
'wse_difference_map.png',
'polygon_wse_difference_map.png',
'shape_z_wse_difference_map.png',
'benefit_rise_area_map.png',
],
'path': [str(wse_png), str(polygon_wse_png), str(shape_z_wse_png), str(benefit_rise_png)],
})
| result | features | cells | area_acres | min_wse_difference_ft | max_wse_difference_ft | |
|---|---|---|---|---|---|---|
| 0 | benefit_polygons | 1 | 73 | 381.28621 | NaN | NaN |
| 1 | rise_polygons | 1 | 832 | 4690.80324 | NaN | NaN |
| 2 | difference_points | 5198 | 905 | NaN | -0.078979 | 0.246094 |
| scenario | min_wse_difference_ft | max_wse_difference_ft | mean_wse_difference_ft | max_abs_wse_difference_ft | cells_abs_gt_0_01_ft | |
|---|---|---|---|---|---|---|
| 0 | high-ground proposed | -0.078979 | 0.246094 | 8.088239e-03 | 0.246094 | 905 |
| 1 | polygon proposed | -0.000488 | 0.002930 | 5.094874e-05 | 0.002930 | 0 |
| 2 | shape_z polygon proposed | -0.000061 | 0.000000 | -5.871023e-08 | 0.000061 | 0 |
Text Only
high-ground proposed cells with |WSE difference| > 0.01 ft: 905
polygon proposed cells with |WSE difference| > 0.01 ft: 0
shape_z polygon proposed cells with |WSE difference| > 0.01 ft: 0
| artifact | path | |
|---|---|---|
| 0 | wse_difference_map.png | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 1 | polygon_wse_difference_map.png | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 2 | shape_z_wse_difference_map.png | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 3 | benefit_rise_area_map.png | G:\GH\ras-commander-pr-docs\examples\example_p... |
Python
# Final artifact list for human review
pd.DataFrame({
'artifact': [
'terrain_profile_difference.png',
'terrain_profile_readback_validation.png',
'polygon_multipoint_control_points.png',
'polygon_profile_readback_validation.png',
'shape_z_polygon_readback_validation.png',
'terrain_difference_map.png',
'wse_difference_map.png',
'polygon_wse_difference_map.png',
'shape_z_wse_difference_map.png',
'benefit_rise_area_map.png',
'base plan HDF',
'high-ground proposed plan HDF',
'polygon proposed plan HDF',
'shape_z proposed plan HDF',
'shape_z terrain HDF',
],
'path': [
str(profile_png),
str(readback_profile_png),
str(polygon_points_png),
str(polygon_profile_png),
str(shape_z_profile_png),
str(terrain_png),
str(wse_png),
str(polygon_wse_png),
str(shape_z_wse_png),
str(benefit_rise_png),
str(base_plan_hdf),
str(proposed_plan_hdf),
str(polygon_plan_hdf),
str(shape_z_plan_hdf),
str(shape_z_terrain_hdf),
],
})
| artifact | path | |
|---|---|---|
| 0 | terrain_profile_difference.png | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 1 | terrain_profile_readback_validation.png | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 2 | polygon_multipoint_control_points.png | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 3 | polygon_profile_readback_validation.png | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 4 | shape_z_polygon_readback_validation.png | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 5 | terrain_difference_map.png | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 6 | wse_difference_map.png | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 7 | polygon_wse_difference_map.png | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 8 | shape_z_wse_difference_map.png | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 9 | benefit_rise_area_map.png | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 10 | base plan HDF | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 11 | high-ground proposed plan HDF | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 12 | polygon proposed plan HDF | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 13 | shape_z proposed plan HDF | G:\GH\ras-commander-pr-docs\examples\example_p... |
| 14 | shape_z terrain HDF | G:\GH\ras-commander-pr-docs\examples\example_p... |