RAS Terrain Creation

This notebook maps the official HEC-RAS tutorial workflow to the current ras-commander terrain APIs.

Official tutorial: https://www.hec.usace.army.mil/confluence/rasdocs/hgt/latest/tutorials/terrain/creating-a-ras-terrain

Primary workflow: 1. Inspect the project projection and terrain catalog from .rasmap. 2. Create a RAS Terrain HDF from a single raster with RasTerrain.create_terrain_hdf(). 3. Optionally query and download USGS 3DEP tiles with Usgs3depAws. 4. Mosaic downloaded tiles with Usgs3depAws.create_vrt() and RasTerrain.vrt_to_tiff(). 5. Register the terrain in RAS Mapper with RasMap.add_terrain_layer(). 6. Demonstrate multi-source terrain priority ordering with RasTerrain.create_terrain_from_rasters() when a channel raster is available.

Prerequisites: - HEC-RAS 6.6+ installed for RasProcess.exe CreateTerrain and bundled GDAL tools. - Python packages: rasterio, geopandas, h5py, pyproj. - Large terrain creation cells are opt-in so the notebook can be reviewed without starting long downloads or terrain builds.

Official Tutorial Coverage

The HEC tutorial covers single-raster terrain creation, multi-raster channel-over-overbank terrain creation, and USGS terrain downloads. Current ras-commander coverage is:

Tutorial operation	ras-commander status	Current API or follow-up
Set project projection from projection.prj	Partial	RasMap.add_terrain_layer(..., projection_prj=...) can update the .rasmap projection while registering terrain. A standalone projection setter is tracked by CLB-270.
Create terrain from base.tif	Covered	RasTerrain.create_terrain_hdf() creates a terrain HDF from one or more rasters.
Multi-source terrain priority ordering	Covered	RasTerrain.create_terrain_from_rasters() and create_terrain_hdf() preserve input raster order; first raster has highest priority.
Query/download USGS terrain tiles	Partial	Usgs3depAws supports USGS 3DEP project discovery and 1m AWS tile downloads. RAS Mapper-style product type/year filtering and non-1m download parity are tracked by CLB-271.
Export XS channel bathymetry as GeoTIFF	Missing	Covered by CLB-253 because the same channel-raster export is needed by the Terrain Modification tutorial.
Terrain display settings	Missing	Hillshade, contour, and stitch TIN edge display settings are tracked by CLB-272.

This notebook keeps missing GUI-only operations guarded or documented instead of bypassing ras-commander with direct RAS executables.

Setup and Imports

Python

from pathlib import Path

# =============================================================================
# DEVELOPMENT MODE TOGGLE
# =============================================================================
# Set USE_LOCAL_SOURCE based on your setup:
#   True  = Use local source code (for developers editing ras-commander)  
#   False = Use pip-installed package (for users)
# =============================================================================

USE_LOCAL_SOURCE = True  # <-- TOGGLE THIS

# -----------------------------------------------------------------------------
if USE_LOCAL_SOURCE:
    import sys
    from pathlib import Path
    local_path = str(Path.cwd().parent)  # Parent of examples/ = repo root
    if local_path not in sys.path:
        sys.path.insert(0, local_path)  # Insert at position 0 = highest priority
    print("LOCAL SOURCE MODE: Loading ras_commander from the repository checkout")
else:
    print("PIP PACKAGE MODE: Loading installed ras-commander")

# Import ras-commander
from ras_commander import RasExamples, RasMap, RasProcess, init_ras_project, ras
from ras_commander.terrain import RasTerrain, Usgs3depAws

# Verify which version loaded
import ras_commander
print(f"Loaded package: {Path(ras_commander.__file__).parent.name}/{Path(ras_commander.__file__).name}")

Text Only

LOCAL SOURCE MODE: Loading ras_commander from the repository checkout
Loaded package: ras_commander/__init__.py

Python

# Additional imports for terrain operations
from pathlib import Path
import os
import h5py
import numpy as np
import subprocess
import shutil
import warnings

NOTEBOOK_DIR = Path.cwd()
REPO_ROOT = NOTEBOOK_DIR.parent
RUN_ROOT = Path(
    os.environ.get(
        "RAS_COMMANDER_EXAMPLE_RUN_ROOT",
        REPO_ROOT / "working" / "example_920_terrain_creation",
    )
).resolve()
RUN_ROOT.mkdir(parents=True, exist_ok=True)
EXAMPLE_PROJECTS_ROOT = RUN_ROOT / "example_projects"
USGS_CACHE_FOLDER = RUN_ROOT / "cache" / "usgs3dep"

def run_path(path):
    """Return a portable label for a path created by this notebook."""
    return str(Path("run") / Path(path).resolve().relative_to(RUN_ROOT))

# Check for required packages
try:
    import rasterio
    from rasterio.crs import CRS
    RASTERIO_AVAILABLE = True
    print("rasterio available")
except ImportError:
    RASTERIO_AVAILABLE = False
    print("rasterio not available - install with: pip install rasterio")
    print("Cannot proceed without rasterio")

try:
    import geopandas as gpd
    GEOPANDAS_AVAILABLE = True
    print("geopandas available")
except ImportError:
    GEOPANDAS_AVAILABLE = False
    print("geopandas not available - install with: pip install geopandas")

print("create_vrt() uses HEC-RAS bundled GDAL tools")
print(f"Run root: {RUN_ROOT.name}")
print(f"Example projects root: {run_path(EXAMPLE_PROJECTS_ROOT)}")
print(f"USGS cache folder: {run_path(USGS_CACHE_FOLDER)}")

Text Only

rasterio available
geopandas available
create_vrt() uses HEC-RAS bundled GDAL tools
Run root: r920_203417
Example projects root: run\example_projects
USGS cache folder: run\cache\usgs3dep

Parameters

Configure these values to customize the notebook for your environment.

Python

# =============================================================================
# PARAMETERS - Edit these to customize the notebook
# =============================================================================

# Project Configuration
PROJECT_NAME = "BaldEagleCrkMulti2D"  # Example project with 2D terrain
NOTEBOOK_SUFFIX = "920"  # Keeps this notebook's extracted project separate
RAS_VERSION = "7.0"  # HEC-RAS version

# Execution guards
# Keep these False for routine notebook QA. Set them True when intentionally
# downloading terrain or invoking HEC-RAS RasProcess terrain creation.
RUN_USGS_DOWNLOAD = False
RUN_TERRAIN_CREATION = False
RUN_MULTI_SOURCE_TERRAIN_CREATION = False

# Terrain Download Configuration
TARGET_RESOLUTION = 1  # meters; current download implementation supports 1m products
NEW_TERRAIN_NAME = "Terrain_3DEP"  # Name for a new terrain created from downloaded data
VERTICAL_UNITS = "Feet"  # Feet or Meters
CREATE_OVERVIEWS = False  # Large mosaics can spend a long time building TIFF overviews
TERRAIN_HDF_TIMEOUT_SECONDS = 14400  # Large full-extent terrain builds can take hours

# HEC-RAS Installation Path (auto-detected by RasTerrain when commands run)
HECRAS_PATH = Path(f"C:/Program Files (x86)/HEC/HEC-RAS/{RAS_VERSION}")
RAS_PROCESS_EXE = HECRAS_PATH / "RasProcess.exe"

print(f"Project: {PROJECT_NAME}")
print(f"Notebook suffix: {NOTEBOOK_SUFFIX}")
print(f"HEC-RAS Version: {RAS_VERSION}")
print(f"Target Resolution: {TARGET_RESOLUTION}m")
print(f"Run USGS download: {RUN_USGS_DOWNLOAD}")
print(f"Run terrain creation: {RUN_TERRAIN_CREATION}")
print(f"Run multi-source terrain creation: {RUN_MULTI_SOURCE_TERRAIN_CREATION}")
print(f"Create overviews: {CREATE_OVERVIEWS}")
print(f"Terrain HDF timeout: {TERRAIN_HDF_TIMEOUT_SECONDS} seconds")
print(f"RasProcess.exe exists: {RAS_PROCESS_EXE.exists()}")

Text Only

Project: BaldEagleCrkMulti2D
Notebook suffix: 920
HEC-RAS Version: 7.0
Target Resolution: 1m
Run USGS download: False
Run terrain creation: False
Run multi-source terrain creation: False
Create overviews: False
Terrain HDF timeout: 14400 seconds
RasProcess.exe exists: True

Extract Example Project

We'll use a stable extracted project under the task-local run root selected by RAS_COMMANDER_EXAMPLE_RUN_ROOT. If this notebook has already downloaded the full Bald Eagle terrain tiles, reusing that workspace preserves the real-project cache and avoids repeating the 98-tile download on subsequent reruns.

Python

# Reuse a stable extracted project for this notebook when available.
project_folder_name = f"{PROJECT_NAME}_{NOTEBOOK_SUFFIX}"
project_path = EXAMPLE_PROJECTS_ROOT / project_folder_name

if project_path.exists():
    print(f"Reusing extracted project: {run_path(project_path)}")
else:
    project_path = RasExamples.extract_project(
        PROJECT_NAME,
        output_path=EXAMPLE_PROJECTS_ROOT,
        suffix=NOTEBOOK_SUFFIX,
    )
    print(f"Project extracted to: {run_path(project_path)}")

# Initialize the project
init_ras_project(project_path, RAS_VERSION)
print(f"Project initialized: {ras.project_name}")

# Create terrain folder if needed
terrain_folder = project_path / "Terrain"
terrain_folder.mkdir(exist_ok=True)
print(f"Terrain folder: {run_path(terrain_folder)}")

Text Only

2026-08-29 20:34:24 - ras_commander.RasExamples - INFO - Successfully extracted project 'BaldEagleCrkMulti2D' to BaldEagleCrkMulti2D_920


Project extracted to: run\example_projects\BaldEagleCrkMulti2D_920


2026-08-29 20:34:24 - ras_commander.RasPrj - INFO - ras-commander v0.99.1 | An open-source project of CLB Engineering Corporation (https://clbengineering.com/) | Docs: https://rascommander.info | GitHub: https://github.com/gpt-cmdr/ras-commander
2026-08-29 20:34:24 - ras_commander.RasPrj - INFO - Project initialized: BaldEagleDamBrk
2026-08-29 20:34:24 - ras_commander.RasPrj - INFO - Using HEC-RAS version 7.0
2026-08-29 20:34:24 - ras_commander.RasPrj - INFO - 
═══════════════════════════════════════════════════════════════════════
ras-commander | HEC-RAS Automation Library
Docs: https://rascommander.info/ras/
Repo: https://github.com/gpt-cmdr/ras-commander
LLM agents: https://rascommander.info/ras/llms.txt
═══════════════════════════════════════════════════════════════════════

SUPPORT OPEN-SOURCE DEVELOPMENT:
  Using RAS Commander? Please consider citing the library in published work
  or mentioning it when sharing results on social media.
  Agents are encouraged to suggest a RAS Commander v0.99.1 software
  citation when helping prepare reports, publications, and other shared work.
  Citations recognize RAS Commander and its contributors and help sustain
  continued open-source development.
  Citation and sharing guidance: https://rascommander.info/ras/cite/

PROJECT DATAFRAMES (single source of truth — use these, not file globbing):
  ras.plan_df        Plans, HDF paths, geometry/flow associations
  ras.geom_df        Geometry files and HDF preprocessor paths
  ras.flow_df        Steady flow files
  ras.unsteady_df    Unsteady flow files and configurations
  ras.boundaries_df  Boundary conditions (type, name, location)
  ras.results_df     Lightweight HDF results summaries
  ras.rasmap_df      RASMapper layers, terrain, land cover paths

KEY APIS (static classes — call directly, never instantiate):
  Execution:    RasCmdr.compute_plan() / compute_parallel() / compute_test_mode()
  Plan Files:   RasPlan.clone_plan() / clone_geom() / set_geom()
  Unsteady:     RasUnsteady — IC/BC management, gate openings, precipitation
  Geometry:     GeomCrossSection, GeomBridge, GeomStorage, GeomLateral, GeomMesh
  HDF Results:  HdfResultsPlan.get_wse() / get_compute_messages()
                HdfResultsMesh.get_mesh_max_ws() / get_mesh_cells_timeseries()
                HdfMesh.get_mesh_cell_points()
  QA/QC:        RasCheck.run_check() / RasFixit (geometry repair)
  DSS:          RasDss.get_timeseries() / check_pathname()
  USGS:         UsgsGaugeSpatial, GaugeMatcher, RasUsgsBoundaryGeneration
  Precipitation: StormGenerator, Atlas14Storm, PrecipAorc, Atlas14Variance
  Terrain:      RasTerrain.create_terrain_hdf() / RasTerrainMod

MULTI-PROJECT: Pass ras_object= to all API calls when using local RasPrj instances.

EXAMPLES: 100+ notebooks in examples/ (100s=execution, 200s=geometry, 300s=unsteady,
  400s=HDF results, 500s=remote, 800s=QA/QC, 900s=data integration).
  Review relevant notebooks before assembling new workflows.

PLATFORM: Most HEC-RAS operations require Windows. Linux/Wine support for
  headless execution, data access, geometry modification, and preprocessing
  is available via RasProcess (HEC-RAS 6.6+). See ras_commander/RasProcess.py.
  Remote distributed execution: ras_commander/remote/ (PsExec, Docker, SSH, cloud).
═══════════════════════════════════════════════════════════════════════


Project initialized: BaldEagleDamBrk
Terrain folder: run\example_projects\BaldEagleCrkMulti2D_920\Terrain

Get Project Extents

Extract the bounding box from the project's existing terrain or geometry to know where to download USGS data.

Python

# Find existing terrain to get extents
existing_tiffs = list(terrain_folder.glob("*.tif"))
existing_hdfs = list(terrain_folder.glob("*.hdf"))

print(f"Found {len(existing_tiffs)} TIFF files")
print(f"Found {len(existing_hdfs)} HDF files")

# Get bounds from first TIFF if available
if existing_tiffs and RASTERIO_AVAILABLE:
    source_tiff = existing_tiffs[0]
    print(f"\nUsing {source_tiff.name} to determine project extents:")

    with rasterio.open(source_tiff) as src:
        bounds = src.bounds
        crs = src.crs
        print(f"  Bounds (native CRS): {bounds}")
        print(f"  CRS: {crs}")
        print(f"  Resolution: {src.res[0]:.2f} x {src.res[1]:.2f}")
        print(f"  Size: {src.width} x {src.height} pixels")

        # Store for later use
        PROJECT_BOUNDS = bounds
        PROJECT_CRS = crs
else:
    print("No existing TIFF found - will need to specify bounds manually")
    PROJECT_BOUNDS = None
    PROJECT_CRS = None

Text Only

Found 2 TIFF files
Found 1 HDF files

Using Terrain50.baldeagledem.tif to determine project extents:
  Bounds (native CRS): BoundingBox(left=1834327.1955903, bottom=162918.80517492996, right=2149835.693237871, top=414872.9473435675)
  CRS: EPSG:2271
  Resolution: 36.50 x 36.50
  Size: 8643 x 6902 pixels

Python

# Convert bounds to WGS84 (lat/lon) for USGS 3DEP
if PROJECT_BOUNDS and RASTERIO_AVAILABLE:
    from pyproj import Transformer
    from shapely.geometry import box

    # Create transformer from project CRS to WGS84
    transformer = Transformer.from_crs(PROJECT_CRS, "EPSG:4326", always_xy=True)

    # Transform corners
    minx, miny = transformer.transform(PROJECT_BOUNDS.left, PROJECT_BOUNDS.bottom)
    maxx, maxy = transformer.transform(PROJECT_BOUNDS.right, PROJECT_BOUNDS.top)

    # Calculate area
    full_width_deg = maxx - minx
    full_height_deg = maxy - miny
    full_area_sq_km = full_width_deg * 111 * full_height_deg * 85  # Rough approximation

    print(f"Full terrain extent:")
    print(f"  Bounds (WGS84): ({minx:.6f}, {miny:.6f}) to ({maxx:.6f}, {maxy:.6f})")
    print(f"  Approx area: {full_area_sq_km:.1f} sq km")

    # Use full extent for download (Usgs3depAws will download only intersecting tiles)
    DOWNLOAD_BBOX = (minx, miny, maxx, maxy)

    print(f"\nDownload bounding box (WGS84):")
    print(f"  West:  {DOWNLOAD_BBOX[0]:.6f}")
    print(f"  South: {DOWNLOAD_BBOX[1]:.6f}")
    print(f"  East:  {DOWNLOAD_BBOX[2]:.6f}")
    print(f"  North: {DOWNLOAD_BBOX[3]:.6f}")
else:
    print("Cannot determine download bounds without existing terrain")
    DOWNLOAD_BBOX = None

Text Only

Full terrain extent:
  Bounds (WGS84): (-78.233233, 40.612784) to (-77.089997, 41.303441)
  Approx area: 7449.7 sq km

Download bounding box (WGS84):
  West:  -78.233233
  South: 40.612784
  East:  -77.089997
  North: 41.303441

Download USGS 3DEP Terrain

The tutorial downloads USGS terrain tiles through the RAS Mapper GUI. ras-commander provides the Python-native path through Usgs3depAws.

Set RUN_USGS_DOWNLOAD = True in the parameter cell to query and download 1m USGS 3DEP tiles directly from the public AWS bucket.

How it works: 1. Uses a cached USGS tile index when available, otherwise downloads it on first run. 2. Finds projects covering your area. 3. Downloads file lists from intersecting projects. 4. Checks which tiles intersect your area. 5. Downloads only intersecting tiles.

Current gap: RAS Mapper's product type/year filter UI is only partially mirrored. Usgs3depAws.download_tiles() supports project_name and min_year, but full data type/year filtering and non-1m download parity are tracked by CLB-271.

Python

# Find USGS 3DEP projects covering this area
if RUN_USGS_DOWNLOAD and DOWNLOAD_BBOX and GEOPANDAS_AVAILABLE:
    print("Finding USGS 3DEP 1m projects for this area...")
    print("This will download a cached tile index on first use if it is not already present.")
    print(f"Cache folder: {run_path(USGS_CACHE_FOLDER)}")

    try:
        projects = Usgs3depAws.find_tiles_for_bbox(
            bbox=DOWNLOAD_BBOX,
            resolution=TARGET_RESOLUTION,
            cache_folder=USGS_CACHE_FOLDER
        )

        print()
        print(f"Found {len(projects)} project(s) with {TARGET_RESOLUTION}m data:")
        for idx, proj in projects.iterrows():
            project_label = proj.get("project", proj.get("proj_name", proj.get("demname", "unknown")))
            published = proj.get("pub_date", "unknown")
            print(f"  - {project_label} (published: {published})")

        USGS_PROJECTS = projects
    except Exception as e:
        print(f"Could not find projects: {e}")
        USGS_PROJECTS = None
else:
    print("Skipping USGS project query. Set RUN_USGS_DOWNLOAD = True to enable it.")
    USGS_PROJECTS = None

Text Only

Skipping USGS project query. Set RUN_USGS_DOWNLOAD = True to enable it.

Python

# Download 1m DEM tiles from USGS 3DEP (AWS S3)
if RUN_USGS_DOWNLOAD and DOWNLOAD_BBOX and USGS_PROJECTS is not None and len(USGS_PROJECTS) > 0:
    print(f"Downloading {TARGET_RESOLUTION}m DEM tiles for project area...")
    print("Note: This checks each tile in intersecting project(s) for overlap with the full project extent.")
    print("      Existing files in the notebook project folder are reused when sizes match.")
    print()

    try:
        tile_folder = terrain_folder / "usgs_tiles"
        downloaded_tiles = Usgs3depAws.download_tiles(
            bbox=DOWNLOAD_BBOX,
            resolution=TARGET_RESOLUTION,
            output_folder=tile_folder,
            cache_folder=USGS_CACHE_FOLDER
        )

        print()
        print(f"Downloaded {len(downloaded_tiles)} tiles:")
        total_size_mb = 0
        for tile in downloaded_tiles:
            size_mb = tile.stat().st_size / 1024 / 1024
            total_size_mb += size_mb
            print(f"  {tile.name} ({size_mb:.2f} MB)")

        print()
        print(f"Total size: {total_size_mb:.2f} MB")

        DOWNLOADED_TILES = downloaded_tiles
    except Exception as e:
        print(f"Download failed: {e}")
        import traceback
        traceback.print_exc()
        DOWNLOADED_TILES = None
else:
    print("Skipping USGS tile download. Set RUN_USGS_DOWNLOAD = True after reviewing project extent and storage needs.")
    DOWNLOADED_TILES = None

Text Only

Skipping USGS tile download. Set RUN_USGS_DOWNLOAD = True after reviewing project extent and storage needs.

Python

# Create mosaic from downloaded tiles.
# Prefer the VRT workflow first. Usgs3depAws.create_vrt() uses HEC-RAS bundled
# GDAL tools, and rasterio merge is only the fallback if the VRT path fails.
if RUN_USGS_DOWNLOAD and DOWNLOADED_TILES and len(DOWNLOADED_TILES) > 0:
    print(f"Creating mosaic from {len(DOWNLOADED_TILES)} tiles...")
    output_tiff = terrain_folder / f"{NEW_TERRAIN_NAME}.tif"
    vrt_path = terrain_folder / f"{NEW_TERRAIN_NAME}.vrt"

    try:
        print("Attempting VRT-based mosaic with HEC-RAS GDAL tools...")

        vrt = Usgs3depAws.create_vrt(
            tile_files=DOWNLOADED_TILES,
            output_vrt=vrt_path
        )
        print(f"Created VRT mosaic: {vrt.name}")

        print()
        print("Converting VRT to TIFF...")
        print("Note: This may take several minutes for large mosaics")
        if CREATE_OVERVIEWS:
            print("Creating TIFF overviews is enabled")
        else:
            print("Creating TIFF overviews is disabled for faster execution")

        final_tiff = RasTerrain.vrt_to_tiff(
            vrt_path=vrt,
            output_path=output_tiff,
            compression="LZW",
            create_overviews=CREATE_OVERVIEWS,
            hecras_version=RAS_VERSION
        )
    except Exception as e:
        print(f"VRT path failed: {e}")
        print("Falling back to rasterio merge...")

        from rasterio.merge import merge

        src_files = [rasterio.open(tile) for tile in DOWNLOADED_TILES]
        mosaic, out_transform = merge(src_files)

        out_meta = src_files[0].meta.copy()
        out_meta.update({
            "driver": "GTiff",
            "height": mosaic.shape[1],
            "width": mosaic.shape[2],
            "transform": out_transform,
            "compress": "lzw"
        })

        for src in src_files:
            src.close()

        with rasterio.open(output_tiff, "w", **out_meta) as dest:
            dest.write(mosaic)

        final_tiff = output_tiff
        print(f"Merged {len(DOWNLOADED_TILES)} tiles using rasterio")

    print()
    print(f"Saved: {run_path(final_tiff)}")
    print(f"  Size: {final_tiff.stat().st_size / 1024 / 1024:.2f} MB")

    with rasterio.open(final_tiff) as src:
        print(f"  Shape: {src.height} x {src.width} pixels")
        print(f"  CRS: {src.crs}")
        print(f"  Resolution: {src.res[0]:.2f} x {src.res[1]:.2f}")
        print(f"  Bounds: {src.bounds}")

    DOWNLOADED_TIFF = final_tiff
else:
    print("No downloaded tiles to mosaic")
    DOWNLOADED_TIFF = None

# For single-raster terrain creation, use the downloaded mosaic when present.
# Otherwise fall back to an existing project TIFF if the example project includes one.
SOURCE_TERRAIN_TIFF = DOWNLOADED_TIFF
if SOURCE_TERRAIN_TIFF is None and existing_tiffs:
    SOURCE_TERRAIN_TIFF = existing_tiffs[0]
    print(f"Using existing project raster as source terrain: {SOURCE_TERRAIN_TIFF.name}")
elif SOURCE_TERRAIN_TIFF is None:
    print("No source terrain TIFF is available for terrain HDF creation")

Text Only

No downloaded tiles to mosaic
Using existing project raster as source terrain: Terrain50.baldeagledem.tif

Create Terrain HDF

Use RasTerrain.create_terrain_hdf() to create HEC-RAS terrain from the downloaded TIFF.

This uses RasProcess.exe CreateTerrain from HEC-RAS 6.6+ which: - Creates multi-resolution pyramid levels (7 levels, 0-6) - Enables TIN stitching for multi-source terrain - Optimizes tile-based storage for HEC-RAS rendering - May need a much longer timeout than small examples when using the full Bald Eagle extent

Python

# Find or create projection file
prj_files = list(terrain_folder.glob("*.prj"))
print(f"Found {len(prj_files)} projection files")

if prj_files:
    projection_prj = prj_files[0]
    print(f"Using existing: {projection_prj.name}")
elif SOURCE_TERRAIN_TIFF and SOURCE_TERRAIN_TIFF.exists():
    # Generate PRJ from the TIFF
    projection_prj = terrain_folder / "Projection.prj"
    print(f"Generating PRJ from {SOURCE_TERRAIN_TIFF.name}...")

    RasTerrain._generate_prj_from_raster(
        SOURCE_TERRAIN_TIFF,
        projection_prj
    )
    print(f"Created: {projection_prj.name}")
else:
    print("No projection file available - terrain creation will fail")
    projection_prj = None

Text Only

Found 1 projection files
Using existing: Projection.prj

Python

# Create terrain HDF using RasTerrain class
if RUN_TERRAIN_CREATION and SOURCE_TERRAIN_TIFF and SOURCE_TERRAIN_TIFF.exists() and projection_prj and projection_prj.exists():
    output_hdf = terrain_folder / f"{NEW_TERRAIN_NAME}.hdf"

    print("Creating terrain HDF:")
    print(f"  Input TIFF: {SOURCE_TERRAIN_TIFF.name}")
    print(f"  Output HDF: {output_hdf.name}")
    print(f"  Projection: {projection_prj.name}")
    print(f"  Units: {VERTICAL_UNITS}")
    print(f"  Timeout: {TERRAIN_HDF_TIMEOUT_SECONDS} seconds")
    print()

    try:
        result_hdf = RasTerrain.create_terrain_hdf(
            input_rasters=[SOURCE_TERRAIN_TIFF],
            output_hdf=output_hdf,
            projection_prj=projection_prj,
            units=VERTICAL_UNITS,
            stitch=True,
            hecras_version=RAS_VERSION,
            timeout_seconds=TERRAIN_HDF_TIMEOUT_SECONDS,
        )

        print()
        print("SUCCESS! Terrain HDF created")
        print(f"  File: {run_path(result_hdf)}")
        print(f"  Size: {result_hdf.stat().st_size / 1024 / 1024:.2f} MB")

        NEW_TERRAIN_HDF = result_hdf
    except Exception as e:
        print()
        print("ERROR: Terrain creation failed")
        print(f"  {e}")
        print()
        print("Troubleshooting:")
        print("  1. Ensure HEC-RAS 6.6+ is installed")
        print(f"  2. Check RasProcess.exe exists: {RAS_PROCESS_EXE.exists()}")
        print("  3. Verify TIFF file is valid (open in QGIS/ArcGIS)")
        print("  4. Increase TERRAIN_HDF_TIMEOUT_SECONDS for very large full-extent terrain builds")
        NEW_TERRAIN_HDF = None
else:
    print("Skipping terrain HDF creation. Set RUN_TERRAIN_CREATION = True to invoke RasProcess.exe.")
    if not SOURCE_TERRAIN_TIFF or not SOURCE_TERRAIN_TIFF.exists():
        print("  Missing: source terrain TIFF file")
    if not projection_prj or not projection_prj.exists():
        print("  Missing: Projection PRJ file")
    NEW_TERRAIN_HDF = None

Text Only

Skipping terrain HDF creation. Set RUN_TERRAIN_CREATION = True to invoke RasProcess.exe.

Multi-Source Terrain Priority Ordering

The official tutorial creates a terrain from channel bathymetry over the overbank/base raster. ras-commander supports that merge when both rasters already exist: pass the channel raster first so it has priority where rasters overlap.

The missing part is creating Channel.tif from cross-section channel geometry. That export is tracked with CLB-253 because it is the same gap used by the Terrain Modification tutorial.

Python

# Demonstrate the multi-source terrain merge API when a channel raster is available.
# Expected priority order: channel raster first, base/overbank raster second.
CHANNEL_TERRAIN_TIFF = terrain_folder / "Channel.tif"
BASE_TERRAIN_TIFF = SOURCE_TERRAIN_TIFF
WITH_CHANNEL_HDF = None

if RUN_MULTI_SOURCE_TERRAIN_CREATION and CHANNEL_TERRAIN_TIFF.exists() and BASE_TERRAIN_TIFF and BASE_TERRAIN_TIFF.exists():
    print("Creating multi-source terrain with channel priority:")
    print(f"  1. Channel raster: {CHANNEL_TERRAIN_TIFF.name}")
    print(f"  2. Base raster: {BASE_TERRAIN_TIFF.name}")

    if projection_prj and projection_prj.exists():
        output_hdf = terrain_folder / f"{NEW_TERRAIN_NAME}_WithChannel.hdf"
        WITH_CHANNEL_HDF = RasTerrain.create_terrain_hdf(
            input_rasters=[CHANNEL_TERRAIN_TIFF, BASE_TERRAIN_TIFF],
            output_hdf=output_hdf,
            projection_prj=projection_prj,
            units=VERTICAL_UNITS,
            stitch=True,
            hecras_version=RAS_VERSION,
            timeout_seconds=TERRAIN_HDF_TIMEOUT_SECONDS,
        )
    else:
        WITH_CHANNEL_HDF = RasTerrain.create_terrain_from_rasters(
            input_rasters=[CHANNEL_TERRAIN_TIFF, BASE_TERRAIN_TIFF],
            output_folder=terrain_folder,
            terrain_name=f"{NEW_TERRAIN_NAME}_WithChannel",
            units=VERTICAL_UNITS,
            stitch=True,
            hecras_version=RAS_VERSION,
            generate_prj=True,
        )

    print(f"Created multi-source terrain: {run_path(WITH_CHANNEL_HDF)}")
else:
    print("Skipping multi-source terrain creation.")
    print("  Set RUN_MULTI_SOURCE_TERRAIN_CREATION = True when channel/base rasters are available.")
    if not CHANNEL_TERRAIN_TIFF.exists():
        print(f"  Missing channel raster: {CHANNEL_TERRAIN_TIFF.name}")
        print("  Exporting channel bathymetry from XS data is tracked by CLB-253.")
    if not BASE_TERRAIN_TIFF or not BASE_TERRAIN_TIFF.exists():
        print("  Missing base terrain raster")

Text Only

Skipping multi-source terrain creation.
  Set RUN_MULTI_SOURCE_TERRAIN_CREATION = True when channel/base rasters are available.
  Missing channel raster: Channel.tif
  Exporting channel bathymetry from XS data is tracked by CLB-253.

Fallback: Use Existing Terrain

If USGS download or terrain creation failed, use the existing terrain from the example project.

Python

# Fallback to existing terrain if new terrain wasn't created
if NEW_TERRAIN_HDF is None or not NEW_TERRAIN_HDF.exists():
    existing_hdfs = list(terrain_folder.glob("*.hdf"))

    if existing_hdfs:
        NEW_TERRAIN_HDF = existing_hdfs[0]
        print(f"Using existing terrain: {NEW_TERRAIN_HDF.name}")
        print(f"  Size: {NEW_TERRAIN_HDF.stat().st_size / 1024 / 1024:.2f} MB")
        # Update name for RasMap registration
        NEW_TERRAIN_NAME = NEW_TERRAIN_HDF.stem
    else:
        print("No terrain HDF available")
else:
    print(f"Using newly created terrain: {NEW_TERRAIN_HDF.name}")

Text Only

Using existing terrain: Terrain50.hdf
  Size: 4.69 MB

Register Terrain in RAS Mapper

After creating a terrain HDF, register it in the project's .rasmap file so HEC-RAS can use it.

Passing projection_prj= also updates the .rasmap RASProjectionFilename reference, which covers the tutorial's projection assignment during terrain registration. A standalone set_project_projection() style API is tracked by CLB-270.

Python

# Find the .rasmap file
rasmap_files = list(project_path.glob("*.rasmap"))
print(f"Found {len(rasmap_files)} .rasmap files:")
for rm in rasmap_files:
    print(f"  {rm.name} ({rm.stat().st_size / 1024:.1f} KB)")

Text Only

Found 1 .rasmap files:
  BaldEagleDamBrk.rasmap (43.9 KB)

Python

# Inspect existing terrain configuration in .rasmap
import xml.etree.ElementTree as ET

if rasmap_files:
    rasmap_file = rasmap_files[0]
    tree = ET.parse(rasmap_file)
    root = tree.getroot()

    # Find Terrains section
    terrains = root.find(".//Terrains")
    if terrains is not None:
        print(f"Existing Terrains section in {rasmap_file.name}:")

        # List terrain layers
        layers = terrains.findall("Layer")
        print(f"\n  Terrain layers ({len(layers)}):")
        for layer in layers:
            print(f"    - Name: {layer.get('Name')}")
            print(f"      Filename: {layer.get('Filename')}")
    else:
        print("No Terrains section found in .rasmap file")

Text Only

Existing Terrains section in BaldEagleDamBrk.rasmap:

  Terrain layers (1):
    - Name: Terrain50
      Filename: .\Terrain\Terrain50.hdf

Python

# Add the new terrain layer to the .rasmap file
if rasmap_files and NEW_TERRAIN_HDF and NEW_TERRAIN_HDF.exists():
    rasmap_file = rasmap_files[0]

    print(f"Adding terrain layer to {rasmap_file.name}:")
    print(f"  Terrain HDF: {NEW_TERRAIN_HDF.name}")
    print(f"  Layer name: {NEW_TERRAIN_NAME}")

    try:
        RasMap.add_terrain_layer(
            terrain_hdf=NEW_TERRAIN_HDF,
            rasmap_path=rasmap_file,
            layer_name=NEW_TERRAIN_NAME,
            projection_prj=projection_prj if projection_prj and projection_prj.exists() else None
        )
        print(f"\nTerrain layer added successfully!")
    except Exception as e:
        print(f"\nError adding terrain layer: {e}")
elif not NEW_TERRAIN_HDF or not NEW_TERRAIN_HDF.exists():
    print(f"No terrain HDF to add")
else:
    print("No .rasmap files found in project")

Text Only

2026-08-29 20:34:25 - ras_commander.RasMap - INFO - Replaced terrain layer 'Terrain50' in .rasmap


Adding terrain layer to BaldEagleDamBrk.rasmap:
  Terrain HDF: Terrain50.hdf
  Layer name: Terrain50

Terrain layer added successfully!

Python

# Verify the terrain layer was added
if rasmap_files:
    rasmap_file = rasmap_files[0]
    tree = ET.parse(rasmap_file)
    root = tree.getroot()

    terrains = root.find(".//Terrains")
    if terrains is not None:
        layers = terrains.findall("Layer")
        print(f"Updated Terrains section ({len(layers)} layers):")
        for layer in layers:
            name = layer.get('Name')
            marker = " <-- NEW" if name == NEW_TERRAIN_NAME else ""
            print(f"  - {name}: {layer.get('Filename')}{marker}")

Text Only

Updated Terrains section (1 layers):
  - Terrain50: .\Terrain\Terrain50.hdf <-- NEW

List Registered Terrain Layers

RasMap.add_terrain_layer() updates the .rasmap catalog. RasMap.list_terrain_layers() reads that catalog back as a dataframe so automation can confirm the new layer name, source filename, resolved path, checked state, terrain type, resampling method, and surface-on flag without manual XML parsing.

Python

terrain_layers = RasMap.list_terrain_layers(project_path)
print(f"Registered terrain layers: {len(terrain_layers)}")
display(terrain_layers[[
    "name",
    "filename",
    "checked",
    "type",
    "resample_method",
    "surface_on",
]])

Text Only

Registered terrain layers: 1

	name	filename	checked	type	resample_method	surface_on
0	Terrain50	.\Terrain\Terrain50.hdf	True	TerrainLayer	near	True

Export an Already Registered Terrain

The creation APIs above and native registered-terrain export solve different problems:

• Usgs3depAws.create_vrt() and RasTerrain.vrt_to_tiff() consolidate loose raster files before terrain creation. RasTerrain.create_terrain_hdf() and RasTerrain.create_terrain_from_rasters() then build a new RAS Terrain HDF, which can be registered with RasMap.add_terrain_layer().
• RasTerrain.export_rasmapper_terrain() starts from an already registered RAS Mapper terrain. It delegates RAS Mapper's Export to Single Raster operation so the selected layer's source ordering, different source resolutions, priorities, stitches, masks, and optional vector modifications remain native RAS Mapper behavior. A requested 1x, 2x, 4x, or 8x factor defines the output resolution from the authoritative source grid.

The guarded example below selects an exact layer from RasMap.list_terrain_layers(), derives a small central window from the source bounds, and writes a validated derivative without registering it back into the project. The opt-in flag reads RAS_COMMANDER_RUN_REGISTERED_TERRAIN_EXPORT, so routine runs remain safe while reproducible qualification runs can set that environment variable to 1. It uses HEC-RAS 6.6 explicitly because that release is qualified for this export; RAS_VERSION above remains the independent version used by the rest of this terrain-creation tutorial. See 931_native_rasmapper_terrain_export.ipynb for the complete bounded workflow and receipt, and 316_terrain_modifications.ipynb for real modification-off/on pixel evidence.

This export does not replace the loose-raster creation APIs, and it does not create the cross-section channel bathymetry raster tracked by CLB-253.

Python

# Keep native export opt-in: it launches the qualified RAS Mapper helper.
RUN_REGISTERED_TERRAIN_EXPORT = (
    os.environ.get("RAS_COMMANDER_RUN_REGISTERED_TERRAIN_EXPORT", "0") == "1"
)
REGISTERED_TERRAIN_EXPORT_RAS_VERSION = "6.6"
REGISTERED_TERRAIN_EXPORT_DOWNSAMPLE_FACTOR = 2
REGISTERED_TERRAIN_EXPORT_NAME = None  # Set an exact registered layer name to override selection.
registered_export_result = None

if RUN_REGISTERED_TERRAIN_EXPORT:
    available_names = terrain_layers["name"].astype(str).tolist()
    if REGISTERED_TERRAIN_EXPORT_NAME is not None:
        selected_terrain_name = REGISTERED_TERRAIN_EXPORT_NAME
    elif NEW_TERRAIN_NAME in available_names:
        selected_terrain_name = NEW_TERRAIN_NAME
    elif PROJECT_NAME == "Muncie" and "TerrainWithChannel" in available_names:
        selected_terrain_name = "TerrainWithChannel"
    else:
        checked_names = terrain_layers.loc[terrain_layers["checked"], "name"].astype(str).tolist()
        if len(checked_names) == 1:
            selected_terrain_name = checked_names[0]
        elif len(available_names) == 1:
            selected_terrain_name = available_names[0]
        else:
            raise RuntimeError("Set REGISTERED_TERRAIN_EXPORT_NAME to one exact registered terrain.")

    selected_rows = terrain_layers.loc[terrain_layers["name"] == selected_terrain_name]
    if len(selected_rows) != 1:
        raise RuntimeError(
            f"Expected exactly one registered terrain named {selected_terrain_name!r}; "
            f"found {len(selected_rows)}"
        )
    if PROJECT_BOUNDS is None:
        raise RuntimeError("A source raster extent is required for the bounded export.")

    center_x = (PROJECT_BOUNDS.left + PROJECT_BOUNDS.right) / 2.0
    center_y = (PROJECT_BOUNDS.bottom + PROJECT_BOUNDS.top) / 2.0
    half_width = min((PROJECT_BOUNDS.right - PROJECT_BOUNDS.left) * 0.05, 500.0)
    half_height = min((PROJECT_BOUNDS.top - PROJECT_BOUNDS.bottom) * 0.05, 500.0)
    registered_export_extent = (
        center_x - half_width,
        center_y - half_height,
        center_x + half_width,
        center_y + half_height,
    )
    registered_export_folder = RUN_ROOT / "registered_terrain_export"
    registered_export_folder.mkdir(parents=True, exist_ok=True)
    registered_export_tif = registered_export_folder / (
        f"{selected_terrain_name}_{REGISTERED_TERRAIN_EXPORT_DOWNSAMPLE_FACTOR}x.tif"
    )

    registered_export_result = RasTerrain.export_rasmapper_terrain(
        ras_project_path=project_path,
        output_tif=registered_export_tif,
        terrain_name=selected_terrain_name,
        extent=registered_export_extent,
        downsample_factor=REGISTERED_TERRAIN_EXPORT_DOWNSAMPLE_FACTOR,
        rasterize_modifications=True,
        overwrite=False,
        timeout_seconds=300.0,
        hecras_version=REGISTERED_TERRAIN_EXPORT_RAS_VERSION,
    )
    if not registered_export_result:
        raise RuntimeError(registered_export_result.error)

    validation = registered_export_result.validation
    assert validation["driver"] == "GTiff"
    assert validation["band_count"] == 1
    assert validation["data_type"] == "Float32"
    assert validation["crs_present"] is True
    assert validation["sidecars"] == []
    assert registered_export_result.source_inventory["intersects_output"].all()
    assert int(registered_export_result.source_inventory["authoritative_grid"].sum()) == 1

    evidence_rows = []
    for _, source in registered_export_result.source_inventory.sort_values("index").iterrows():
        evidence_rows.append(
            {
                "source": Path(source["filename"]).name,
                "level_zero_cell_size": float(source["cell_sizes"][0]),
                "authoritative": bool(source["authoritative_grid"]),
                "intersects": bool(source["intersects_output"]),
            }
        )
    print(f"Native export: {selected_terrain_name} at {registered_export_result.output_cell_size:g} units")
    print(f"Output: {run_path(registered_export_result.output_path)}")
    print(f"Receipt: {run_path(registered_export_result.receipt_path)}")
    print(f"Elapsed: {registered_export_result.elapsed_seconds:.3f} seconds")
    display(evidence_rows)

    import matplotlib.pyplot as plt
    with rasterio.open(registered_export_result.output_path) as dataset:
        export_array = dataset.read(1, masked=True)
        export_bounds = dataset.bounds

    fig, (ax_raster, ax_evidence) = plt.subplots(
        1, 2, figsize=(11.5, 4.8), gridspec_kw={"width_ratios": [1.3, 1.0]}
    )
    image = ax_raster.imshow(
        export_array,
        extent=(export_bounds.left, export_bounds.right, export_bounds.bottom, export_bounds.top),
        cmap="terrain",
        aspect="auto",
    )
    ax_raster.set_title("Bounded RAS Mapper terrain")
    ax_raster.set_xlabel("Project X")
    ax_raster.set_ylabel("Project Y")
    fig.colorbar(image, ax=ax_raster, shrink=0.82, label="Elevation")

    source_lines = [
        f"{row['source']}: {row['level_zero_cell_size']:g}"
        + (" (authoritative grid)" if row["authoritative"] else "")
        for row in evidence_rows
    ]
    ax_evidence.axis("off")
    ax_evidence.set_title("Source and validation evidence")
    ax_evidence.text(
        0.0, 0.95,
        f"Terrain: {selected_terrain_name}\n"
        f"Sources: {len(evidence_rows)}\n"
        f"Output: {validation['columns']} x {validation['rows']} pixels\n"
        f"Cell size: {registered_export_result.output_cell_size:g}\n"
        f"Factor: {registered_export_result.downsample_factor}x\n"
        f"Band: {validation['data_type']}\n"
        f"CRS: {validation['crs_present']}\n"
        f"Checksum: {validation['checksum']}\n\n" + "\n".join(source_lines),
        va="top", family="monospace", fontsize=9.3, linespacing=1.3,
    )
    fig.suptitle(f"{selected_terrain_name}: native {REGISTERED_TERRAIN_EXPORT_DOWNSAMPLE_FACTOR}x export")
    fig.tight_layout()
    plt.show()
else:
    print("Registered-terrain export skipped. Set RAS_COMMANDER_RUN_REGISTERED_TERRAIN_EXPORT=1 to run it.")

Text Only

Native export: Terrain50 at 40 units
Output: run\registered_terrain_export\Terrain50_2x.tif
Receipt: run\registered_terrain_export\Terrain50_2x.tif.receipt.json
Elapsed: 2.561 seconds



[{'source': 'Terrain50.baldeagledem.tif',
  'level_zero_cell_size': 36.504512049933,
  'authoritative': False,
  'intersects': True},
 {'source': 'Terrain50.dtm_20ft.tif',
  'level_zero_cell_size': 20.0,
  'authoritative': True,
  'intersects': True}]

Associate Terrain With Geometry HDF

Registering a terrain in .rasmap makes it available to RAS Mapper. Associating a terrain with a compiled geometry HDF writes the HEC-RAS /Geometry attributes used by geometry preprocessing/property-table workflows. This is intentionally Python-native for normal automation; it requires an existing .g##.hdf and does not compile a text .g## file into HDF.

Python

geometry_hdf_path = None
if getattr(ras, "geom_df", None) is not None and not ras.geom_df.empty:
    mesh_rows = ras.geom_df
    if "has_2d_mesh" in mesh_rows.columns:
        mesh_rows = mesh_rows[mesh_rows["has_2d_mesh"]]
    if not mesh_rows.empty and "hdf_path" in mesh_rows.columns:
        candidate = Path(mesh_rows.iloc[0]["hdf_path"])
        if candidate.exists():
            geometry_hdf_path = candidate

if geometry_hdf_path and NEW_TERRAIN_HDF and NEW_TERRAIN_HDF.exists():
    print(f"Associating terrain with geometry HDF: {geometry_hdf_path.name}")
    RasMap.associate_geometry_layers(
        project_path,
        geometry_hdf_path,
        terrain_hdf_path=NEW_TERRAIN_HDF,
    )

    association = RasMap.get_hdf_geometry_association(geometry_hdf_path)
    print("Updated geometry association:")
    for key in ["terrain_hdf_path", "terrain_layer_name", "terrain_date"]:
        value = association.get(key)
        if key.endswith("_path") and value:
            value = Path(value).name
        print(f"  {key}: {value}")
else:
    print("Skipping geometry association because no existing geometry HDF or terrain HDF was found.")

Text Only

Associating terrain with geometry HDF: BaldEagleDamBrk.g06.hdf


2026-08-29 20:34:32 - ras_commander.geom.GeomMesh - INFO - Updated geometry associations on BaldEagleDamBrk.g06.hdf: infiltration_hdf_path, landcover_hdf_path, terrain_hdf_path


Updated geometry association:
  terrain_hdf_path: Terrain50.hdf
  terrain_layer_name: Terrain50
  terrain_date: None

Optional Native CLI Validator

RasProcess.validate_geometry_association_cli() is a reference/validation path around RasProcess.exe SetGeometryAssociation. It mutates the supplied HDF in place, so run it only on a disposable copy or a model you intentionally want the native executable to touch. The main public workflow remains RasMap.associate_geometry_layers() above.

Python

RUN_NATIVE_ASSOCIATION_VALIDATOR = False

if RUN_NATIVE_ASSOCIATION_VALIDATOR and geometry_hdf_path and NEW_TERRAIN_HDF:
    import pandas as pd

    cli_result = RasProcess.validate_geometry_association_cli(
        geometry_hdf_path,
        terrain_hdf_path=NEW_TERRAIN_HDF,
        ras_version=RAS_VERSION,
    )
    print(f"Native validator passed: {cli_result['passed']}")
    print(f"Return code: {cli_result['return_code']}")
    if cli_result["mismatches"]:
        display(pd.DataFrame(cli_result["mismatches"]))
else:
    print("Native RasProcess.exe validator skipped. Set RUN_NATIVE_ASSOCIATION_VALIDATOR = True to run it.")

Text Only

Native RasProcess.exe validator skipped. Set RUN_NATIVE_ASSOCIATION_VALIDATOR = True to run it.

Summary

This notebook maps the official HEC-RAS Creating a RAS Terrain tutorial to current ras-commander APIs:

1. Official Tutorial Reference - Linked the USACE HEC-RAS tutorial URL and mapped tutorial operations to API coverage.

2. Extract Project - Gets project extents from existing terrain using a real example project.

3. Download USGS 3DEP - Uses Usgs3depAws to identify and download intersecting 1m terrain tiles from AWS S3 when RUN_USGS_DOWNLOAD = True.

4. Create Mosaic - Uses Usgs3depAws.create_vrt() and RasTerrain.vrt_to_tiff():

5. Uses HEC-RAS bundled GDAL tools.
6. Falls back to rasterio merge only if the VRT path fails.

7. Optionally builds TIFF overviews.

8. Create Terrain HDF - Uses RasTerrain.create_terrain_hdf() when RUN_TERRAIN_CREATION = True:

9. Multi-resolution pyramid levels.
10. TIN stitching for seamless terrain.

11. Optimized for HEC-RAS rendering.

12. Multi-Source Terrain - Demonstrates RasTerrain.create_terrain_from_rasters() priority ordering when a channel raster is already available.

13. Register in RAS Mapper - Uses RasMap.add_terrain_layer() and optional projection_prj= to update terrain and project projection references.

14. List Terrain Layers - Uses RasMap.list_terrain_layers() for dataframe-based discovery.

15. Export Registered Terrain - Uses a bounded, opt-in RasTerrain.export_rasmapper_terrain() call to consolidate an already registered terrain through native RAS Mapper behavior.

16. Associate Geometry - Uses RasMap.associate_geometry_layers() to write terrain association attributes to an existing geometry HDF.

17. Native Validation Option - Includes guarded RasProcess.validate_geometry_association_cli() usage for disposable/reference checks.

Remaining Feature Gaps

• CLB-270: standalone project projection assignment API.
• CLB-271: RAS Mapper-style USGS product type/year filtering and non-1m download parity.
• XS channel bathymetry GeoTIFF export, linked to CLB-253.
• CLB-272: terrain hillshade, contour, and stitch TIN edge display settings.

Next Steps

• Open project in HEC-RAS and verify terrain displays correctly after an intentional terrain build.
• Use the Linear follow-up issues for missing API work instead of adding GUI-only workarounds to this notebook.

---

Appendix: Terrain HDF Structure Analysis

The following cells analyze the internal structure of HEC-RAS terrain HDF files for reference.

Python

def print_hdf_structure(hdf_file, prefix="", max_depth=3, current_depth=0):
    """Recursively print HDF5 file structure with dataset info."""
    for key in hdf_file.keys():
        item = hdf_file[key]
        if isinstance(item, h5py.Group):
            print(f"{prefix}{key}/")
            if current_depth < max_depth:
                print_hdf_structure(item, prefix + "  ", max_depth, current_depth + 1)
        else:
            # Dataset - show shape and dtype
            shape_str = f"shape={item.shape}" if hasattr(item, 'shape') else ""
            dtype_str = f"dtype={item.dtype}" if hasattr(item, 'dtype') else ""
            print(f"{prefix}{key} ({shape_str}, {dtype_str})")


# Find terrain HDF files to analyze
terrain_hdfs = list(terrain_folder.glob("*.hdf"))
print(f"Found {len(terrain_hdfs)} terrain HDF files:")
for hdf in terrain_hdfs:
    print(f"  {hdf.name} ({hdf.stat().st_size / 1024 / 1024:.2f} MB)")

Text Only

Found 1 terrain HDF files:
  Terrain50.hdf (4.69 MB)

Python

# Open and explore the first terrain HDF
if terrain_hdfs:
    terrain_hdf = terrain_hdfs[0]
    print(f"\nExploring: {terrain_hdf.name}")
    print("=" * 60)

    with h5py.File(terrain_hdf, 'r') as hdf:
        print_hdf_structure(hdf, max_depth=2)
else:
    print("No terrain HDF files found in Terrain folder")

Text Only

Exploring: Terrain50.hdf
============================================================
Terrain/
  Stitch TIN Points (shape=(72592, 4), dtype=float64)
  Stitch TIN Triangles (shape=(72301, 3), dtype=int32)
  Stitches (shape=(103000, 7), dtype=float64)
  Terrain50.baldeagledem/
    0/
    1/
    2/
    3/
    4/
    5/
    6/
  Terrain50.dtm_20ft/
    0/
    1/
    2/
    3/
    4/
    5/

Python

# Inspect the stitching data structure
if terrain_hdfs:
    with h5py.File(terrain_hdfs[0], 'r') as hdf:
        terrain_group = hdf.get('Terrain')

        if terrain_group:
            # Check for stitch data
            if 'Stitch TIN Points' in terrain_group:
                stitch_points = terrain_group['Stitch TIN Points'][:]
                print(f"Stitch TIN Points:")
                print(f"  Shape: {stitch_points.shape}")
                print(f"  Columns: X, Y, Z, M (material/priority)")
                print(f"  Sample (first 3 points):")
                for i, pt in enumerate(stitch_points[:3]):
                    print(f"    [{i}] X={pt[0]:.2f}, Y={pt[1]:.2f}, Z={pt[2]:.2f}, M={pt[3]:.0f}")

            if 'Stitch TIN Triangles' in terrain_group:
                stitch_tris = terrain_group['Stitch TIN Triangles'][:]
                print(f"\nStitch TIN Triangles:")
                print(f"  Shape: {stitch_tris.shape}")
                print(f"  Columns: p0, p1, p2 (vertex indices)")

            # List source terrain layers
            print(f"\nSource terrain layers:")
            for key in terrain_group.keys():
                if key not in ['Stitch TIN Points', 'Stitch TIN Triangles', 'Stitches']:
                    print(f"  {key}")
        else:
            print("No Terrain group found in HDF")

Text Only

Stitch TIN Points:
  Shape: (72592, 4)
  Columns: X, Y, Z, M (material/priority)
  Sample (first 3 points):
    [0] X=2075224.54, Y=371231.14, Z=796.59, M=20
    [1] X=2075224.54, Y=371251.14, Z=803.22, M=20
    [2] X=2075244.54, Y=371251.14, Z=803.09, M=20

Stitch TIN Triangles:
  Shape: (72301, 3)
  Columns: p0, p1, p2 (vertex indices)

Source terrain layers:
  Terrain50.baldeagledem
  Terrain50.dtm_20ft

Python

# Inspect pyramid levels in a source terrain layer
if terrain_hdfs:
    with h5py.File(terrain_hdfs[0], 'r') as hdf:
        terrain_group = hdf.get('Terrain')

        if terrain_group:
            # Find first source layer (not stitch data)
            source_layers = [k for k in terrain_group.keys() 
                            if k not in ['Stitch TIN Points', 'Stitch TIN Triangles', 'Stitches']]

            if source_layers:
                layer_name = source_layers[0]
                layer = terrain_group[layer_name]
                print(f"Source layer: {layer_name}")
                print(f"\nPyramid levels (resolution hierarchy):")

                # List pyramid levels (0-6 typically)
                for level_key in sorted([k for k in layer.keys() if k.isdigit()]):
                    level = layer[level_key]
                    print(f"\n  Level {level_key}:")
                    for dataset_key in level.keys():
                        ds = level[dataset_key]
                        print(f"    {dataset_key}: shape={ds.shape}, dtype={ds.dtype}")

Text Only

Source layer: Terrain50.baldeagledem

Pyramid levels (resolution hierarchy):

  Level 0:
    Mask: shape=(918, 65536), dtype=uint8
    Min-Max: shape=(918, 3), dtype=float32
    Perimeter: shape=(918, 513), dtype=float32

  Level 1:
    Mask: shape=(238, 65536), dtype=uint8
    Min-Max: shape=(238, 3), dtype=float32
    Perimeter: shape=(238, 513), dtype=float32

  Level 2:
    Mask: shape=(63, 65536), dtype=uint8
    Min-Max: shape=(63, 3), dtype=float32
    Perimeter: shape=(63, 513), dtype=float32

  Level 3:
    Mask: shape=(20, 65536), dtype=uint8
    Min-Max: shape=(20, 3), dtype=float32
    Perimeter: shape=(20, 513), dtype=float32

  Level 4:
    Mask: shape=(6, 65536), dtype=uint8
    Min-Max: shape=(6, 3), dtype=float32
    Perimeter: shape=(6, 513), dtype=float32

  Level 5:
    Mask: shape=(2, 65536), dtype=uint8
    Min-Max: shape=(2, 3), dtype=float32
    Perimeter: shape=(2, 513), dtype=float32

  Level 6:
    Mask: shape=(1, 65536), dtype=uint8
    Min-Max: shape=(1, 3), dtype=float32
    Perimeter: shape=(1, 513), dtype=float32

GUI Verification (Optional)

To verify the terrain displays correctly in HEC-RAS:

1. Open the project in HEC-RAS.
2. Launch RAS Mapper.
3. Check that the terrain appears under Terrains.
4. Check the new terrain layer to display it.
5. If using a multi-source terrain, compare the base terrain and channel-priority terrain.

The official tutorial also enables hillshade, contours, and stitch TIN edge display. Those RAS Mapper display settings are tracked by CLB-272.

Cleanup

Python

# Optional: Clean up extracted project
CLEANUP_PROJECT = False  # Set to True to remove extracted project

if CLEANUP_PROJECT and project_path.exists():
    print(f"Cleaning up: {project_path}")
    shutil.rmtree(project_path, ignore_errors=True)
    print("Project folder removed")
else:
    print(f"Project preserved at: {run_path(project_path)}")
    print("Set CLEANUP_PROJECT = True to remove on next run")

Text Only

Project preserved at: run\example_projects\BaldEagleCrkMulti2D_920
Set CLEANUP_PROJECT = True to remove on next run