Creating a RAS Terrain with ras-commander¶

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

# =============================================================================
# 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(f"LOCAL SOURCE MODE: Loading from {local_path}/ras_commander")
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: {ras_commander.__file__}")

Python

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

NOTEBOOK_DIR = Path.cwd()
REPO_ROOT = NOTEBOOK_DIR.parent
EXAMPLE_PROJECTS_ROOT = REPO_ROOT / "example_projects"
USGS_CACHE_FOLDER = REPO_ROOT / "cache" / "usgs3dep"

# 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"Example projects root: {EXAMPLE_PROJECTS_ROOT}")
print(f"USGS cache folder: {USGS_CACHE_FOLDER}")

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()}")

Extract Example Project¶

We'll use a stable extracted project under the repo-root example_projects/ folder. 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: {project_path}")
else:
    project_path = RasExamples.extract_project(
        PROJECT_NAME,
        output_path=EXAMPLE_PROJECTS_ROOT,
        suffix=NOTEBOOK_SUFFIX,
    )
    print(f"Project extracted to: {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: {terrain_folder}")

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

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

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: {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

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

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: {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")

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

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: {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

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: {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}")
        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")

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}")

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)")

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")

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")

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}")

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",
    "resolved_path",
    "checked",
    "type",
    "resample_method",
    "surface_on",
]])

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"]:
        print(f"  {key}: {association.get(key)}")
else:
    print("Skipping geometry association because no existing geometry HDF or terrain HDF was found.")

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.")

Summary¶

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

Official Tutorial Reference - Linked the USACE HEC-RAS tutorial URL and mapped tutorial operations to API coverage.
Extract Project - Gets project extents from existing terrain using a real example project.
Download USGS 3DEP - Uses Usgs3depAws to identify and download intersecting 1m terrain tiles from AWS S3 when RUN_USGS_DOWNLOAD = True.
Create Mosaic - Uses Usgs3depAws.create_vrt() and RasTerrain.vrt_to_tiff():
Uses HEC-RAS bundled GDAL tools.
Falls back to rasterio merge only if the VRT path fails.
Optionally builds TIFF overviews.
Create Terrain HDF - Uses RasTerrain.create_terrain_hdf() when RUN_TERRAIN_CREATION = True:
Multi-resolution pyramid levels.
TIN stitching for seamless terrain.
Optimized for HEC-RAS rendering.
Multi-Source Terrain - Demonstrates RasTerrain.create_terrain_from_rasters() priority ordering when a channel raster is already available.
Register in RAS Mapper - Uses RasMap.add_terrain_layer() and optional projection_prj= to update terrain and project projection references.
List Terrain Layers - Uses RasMap.list_terrain_layers() for dataframe-based discovery.
Associate Geometry - Uses RasMap.associate_geometry_layers() to write terrain association attributes to an existing geometry HDF.
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)")

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")

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")

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}")

GUI Verification (Optional)¶

To verify the terrain displays correctly in HEC-RAS:

Open the project in HEC-RAS.
Launch RAS Mapper.
Check that the terrain appears under Terrains.
Check the new terrain layer to display it.
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: {project_path}")
    print("Set CLEANUP_PROJECT = True to remove on next run")