Skip to content

Headless Stored Map Generation Using RasMapper

This notebook demonstrates headless stored map generation using RasMapper through RasProcess.store_maps() to create: - Maximum Water Surface Elevation (WSE) rasters - Maximum Depth rasters - Inundation boundary polygons

from computed HEC-RAS plans using the BaldEagleCrkMulti2D example project.

Comparison with Other Methods

Method Speed Reliability Cloud-Compatible GUI Required Recommendation
Headless RasMapper stored maps (this notebook) Fastest (8-10 sec) Excellent No No Recommended (Windows)
GUI Automation (notebook 600) Slow (60+ sec) Fragile No Yes Last Resort

Prerequisites

  • HEC-RAS 6.x or 7.x installed with RAS Mapper support
  • Windows operating system
  • Required packages: ras-commander, rasterio, geopandas, shapely
Python
# =============================================================================
# DEVELOPMENT MODE TOGGLE
# =============================================================================
import logging
import warnings

# Keep example output focused on the stored-map workflow.
logging.getLogger().setLevel(logging.WARNING)
logging.getLogger("ras_commander").setLevel(logging.WARNING)
logging.getLogger("rasterio").setLevel(logging.ERROR)
logging.getLogger("rasterio._env").setLevel(logging.ERROR)
logging.getLogger("pyogrio").setLevel(logging.WARNING)
warnings.filterwarnings(
    "ignore",
    message="Setting the shape on a NumPy array has been deprecated.*",
    category=DeprecationWarning,
)

USE_LOCAL_SOURCE = True  # <-- Set to True to use local ras-commander source

if USE_LOCAL_SOURCE:
    import sys
    from pathlib import Path
    local_path = str(Path.cwd().parent)
    if local_path not in sys.path:
        sys.path.insert(0, local_path)
    print("LOCAL SOURCE MODE: Loading repository source")
else:
    print("PIP PACKAGE MODE: Loading installed ras-commander")

# Import ras-commander
from ras_commander import (
    RasProcess, RasMap, RasCmdr, RasExamples,
    init_ras_project, ras
)

# Additional imports
import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio
from rasterio.features import shapes
from rasterio.plot import show
from shapely.geometry import shape
from shapely.ops import unary_union
from pathlib import Path
import matplotlib.pyplot as plt


def notebook_display_path(path):
    """Return a stable relative path for notebook output messages."""
    path = Path(path)
    try:
        return path.resolve().relative_to(Path.cwd().resolve()).as_posix()
    except ValueError:
        return path.name


# Verify which version loaded
import ras_commander
logging.getLogger("ras_commander").setLevel(logging.WARNING)
print(f"Loaded ras-commander {ras_commander.__version__}")
Text Only
LOCAL SOURCE MODE: Loading repository source


Loaded ras-commander 0.98.2

Parameters

Python
# =============================================================================
# PARAMETERS
# =============================================================================

# Example project configuration
PROJECT_NAME = "BaldEagleCrkMulti2D"
RAS_VERSION = "7.0"
PLAN = "06"  # Dam break plan

# Profile to map ("Max", "Min", or specific timestamp)
PROFILE = "Max"

# Depth threshold for inundation boundary (feet)
DEPTH_THRESHOLD = 0.1

print(f"Project: {PROJECT_NAME}")
print(f"Plan: {PLAN}, Profile: {PROFILE}")
print(f"Depth threshold: {DEPTH_THRESHOLD} ft")
Text Only
Project: BaldEagleCrkMulti2D
Plan: 06, Profile: Max
Depth threshold: 0.1 ft

Step 1: Extract Example Project and Initialize

Python
# Extract example project (uses suffix to avoid conflicts with other notebooks)
project_path = RasExamples.extract_project(PROJECT_NAME, suffix="601")

# Initialize project
init_ras_project(project_path, RAS_VERSION)

# Output folder for floodplain maps
OUTPUT_BASE = project_path / "FloodplainMaps"

print(f"Project Name: {ras.project_name}")
print(f"Project Folder: {notebook_display_path(ras.project_folder)}")
print(f"Output folder: {notebook_display_path(OUTPUT_BASE)}")
print(f"\nPlans in project:")
print(ras.plan_df[['plan_number', 'Plan Title']].to_string())
Text Only
Project Name: BaldEagleDamBrk
Project Folder: example_projects/BaldEagleCrkMulti2D_601
Output folder: example_projects/BaldEagleCrkMulti2D_601/FloodplainMaps

Plans in project:
   plan_number                               Plan Title
0           13                  PMF with Multi 2D Areas
1           15              1d-2D Dambreak Refined Grid
2           17                          2D to 1D No Dam
3           18                             2D to 2D Run
4           19                   SA to 2D Dam Break Run
5           03  Single 2D Area - Internal Dam Structure
6           04  SA to 2D Area Conn - 2D Levee Structure
7           02                 SA to Detailed 2D Breach
8           01             SA to Detailed 2D Breach FEQ
9           05          Single 2D area with Bridges FEQ
10          06            Gridded Precip - Infiltration

Step 2: Compute Plan (if needed)

Stored map generation needs computed HDF results. Compute the plan if it hasn't been run yet.

Python
# Check if plan has HDF results, compute if needed
hdf_path = ras.project_folder / f"{ras.project_name}.p{PLAN}.hdf"

if not hdf_path.exists():
    print(f"Computing plan {PLAN}...")
    RasCmdr.compute_plan(PLAN, num_cores=2)
    print(f"Plan {PLAN} complete.")
else:
    print(f"Plan {PLAN} already computed ({hdf_path.stat().st_size / (1024*1024):.1f} MB)")
Text Only
Computing plan 06...


Plan 06 complete.

Step 3: Check .rasmap Compatibility

Python
# Check and upgrade .rasmap if needed for RasProcess compatibility
result = RasMap.ensure_rasmap_compatible(auto_upgrade=True)

print(f"Status: {result['status']}")
print(f"Message: {result['message']}")
print(f"Version: {result['version']}")

if result['status'] == 'manual_needed':
    print("\nManual intervention required:")
    print("1. Open project in HEC-RAS")
    print("2. Click 'GIS Tools' > 'RAS Mapper'")
    print("3. Wait for RASMapper to open (this upgrades .rasmap)")
    print("4. Close RASMapper and HEC-RAS")
    print("5. Re-run this notebook")
else:
    print("\n.rasmap file is ready for stored map generation")
Text Only
Status: ready
Message: Already compatible (version 2.0.0)
Version: 2.0.0

.rasmap file is ready for stored map generation

Step 4: Identify Plans with HDF Results

Python
# Find all plans that have computed HDF results
plans_with_hdf = []

for _, row in ras.plan_df.iterrows():
    plan_num = row['plan_number']
    plan_hdf = ras.project_folder / f"{ras.project_name}.p{plan_num}.hdf"

    if plan_hdf.exists():
        plans_with_hdf.append({
            'plan_number': plan_num,
            'Plan Title': row.get('Plan Title', f'Plan {plan_num}'),
            'hdf_path': plan_hdf,
            'hdf_size_mb': plan_hdf.stat().st_size / (1024 * 1024)
        })

print(f"Found {len(plans_with_hdf)} plans with HDF results:\n")
for p in plans_with_hdf:
    print(f"  Plan {p['plan_number']}: {p['Plan Title']} ({p['hdf_size_mb']:.1f} MB)")
Text Only
Found 1 plans with HDF results:

  Plan 06: Gridded Precip - Infiltration (577.2 MB)

Step 5: Generate Max WSE and Depth Rasters

This uses RasProcess.store_maps() to generate stored maps headlessly through the HEC-RAS RasMapper stack. ras-commander returns the generated raster paths and applies georeferencing fixes when fix_georef=True.

Python
# Create output folder
OUTPUT_BASE.mkdir(parents=True, exist_ok=True)

# Track results
all_results = {}
failed_plans = []

print(f"Generating Max WSE and Depth rasters for {len(plans_with_hdf)} plans...")
print("=" * 70)

for i, plan_info in enumerate(plans_with_hdf):
    plan_num = plan_info['plan_number']
    plan_title = plan_info['Plan Title']

    print(f"\n[{i+1}/{len(plans_with_hdf)}] Processing Plan {plan_num}: {plan_title}")

    try:
        results = RasProcess.store_maps(
            plan_number=plan_num,
            profile=PROFILE,
            wse=True,
            depth=True,
            velocity=False,
            fix_georef=True,
            ras_version=RAS_VERSION,
            timeout=1800
        )

        all_results[plan_num] = results

        wse_count = len(results.get('wse', []))
        depth_count = len(results.get('depth', []))
        print(f"    Generated: {wse_count} WSE, {depth_count} Depth rasters")

    except Exception as e:
        print(f"    ERROR: {e}")
        failed_plans.append({'plan': plan_num, 'error': str(e)})

print("\n" + "=" * 70)
print(f"Completed: {len(all_results)} plans processed successfully")
if failed_plans:
    print(f"Failed: {len(failed_plans)} plans")
    for fp in failed_plans:
        print(f"  - Plan {fp['plan']}: {fp['error']}")
Text Only
Generating Max WSE and Depth rasters for 1 plans...
======================================================================

[1/1] Processing Plan 06: Gridded Precip - Infiltration


    Generated: 1 WSE, 1 Depth rasters

======================================================================
Completed: 1 plans processed successfully

Step 6: Generate Inundation Boundary Polygons

Convert depth rasters to polygons representing the inundation boundary.

Python
def depth_raster_to_polygon(depth_tif_path: Path, depth_threshold: float = 0.1) -> gpd.GeoDataFrame:
    """
    Convert a depth raster to inundation boundary polygon(s).

    Args:
        depth_tif_path: Path to depth raster TIF file
        depth_threshold: Minimum depth to consider as inundated (feet)

    Returns:
        GeoDataFrame with inundation boundary polygon(s)
    """
    with rasterio.open(depth_tif_path) as src:
        depth_data = src.read(1)
        transform = src.transform
        crs = src.crs
        nodata = src.nodata

        # Create binary mask: 1 = inundated (depth > threshold), 0 = dry
        if nodata is not None:
            inundated_mask = (depth_data > depth_threshold) & (depth_data != nodata)
        else:
            inundated_mask = depth_data > depth_threshold

        inundated_mask = inundated_mask.astype(np.uint8)

        # Extract polygon shapes from binary mask
        polygon_shapes = list(shapes(
            inundated_mask,
            mask=inundated_mask == 1,
            transform=transform
        ))

        if not polygon_shapes:
            return gpd.GeoDataFrame(columns=['geometry'], crs=crs)

        # Convert to shapely geometries
        geometries = [shape(geom) for geom, value in polygon_shapes if value == 1]

        if not geometries:
            return gpd.GeoDataFrame(columns=['geometry'], crs=crs)

        # Merge all polygons into a single multipolygon
        merged = unary_union(geometries)

        return gpd.GeoDataFrame({'geometry': [merged]}, crs=crs)


print("Function defined: depth_raster_to_polygon()")
Text Only
Function defined: depth_raster_to_polygon()
Python
# Generate inundation boundary polygons for all plans with depth rasters
inundation_polygons = {}

print(f"Generating inundation boundary polygons (depth > {DEPTH_THRESHOLD} ft)...")
print("=" * 70)

for plan_num, results in all_results.items():
    depth_files = results.get('depth', [])

    if not depth_files:
        print(f"Plan {plan_num}: No depth rasters found, skipping")
        continue

    for depth_tif in depth_files:
        print(f"\nProcessing Plan {plan_num}: {depth_tif.name}")

        try:
            gdf = depth_raster_to_polygon(depth_tif, DEPTH_THRESHOLD)

            if gdf.empty:
                print(f"    No inundation areas found (all depths < {DEPTH_THRESHOLD} ft)")
                continue

            # Calculate area
            area_sq_ft = gdf.geometry.area.sum()
            area_acres = area_sq_ft / 43560

            # Save to shapefile
            output_shp = OUTPUT_BASE / f"Inundation_Boundary_Plan_{plan_num}.shp"
            gdf['plan'] = plan_num
            gdf['area_sqft'] = area_sq_ft
            gdf['area_acres'] = area_acres
            gdf.to_file(output_shp)

            inundation_polygons[plan_num] = {
                'shapefile': output_shp,
                'gdf': gdf,
                'area_acres': area_acres
            }

            print(f"    Inundation area: {area_acres:,.1f} acres")
            print(f"    Saved to: {output_shp.name}")

        except Exception as e:
            print(f"    ERROR: {e}")

print("\n" + "=" * 70)
print(f"Generated {len(inundation_polygons)} inundation boundary polygons")
Text Only
Generating inundation boundary polygons (depth > 0.1 ft)...
======================================================================

Processing Plan 06: Depth (Max).Terrain50.dtm_20ft.tif


    Inundation area: 31,558.5 acres
    Saved to: Inundation_Boundary_Plan_06.shp

======================================================================
Generated 1 inundation boundary polygons

Step 7: Summary Report

Python
# Summary of all generated outputs
print("\n" + "=" * 70)
print("FLOODPLAIN MAPPING SUMMARY")
print("=" * 70)
print(f"\nProject: {ras.project_name}")
print(f"Profile: {PROFILE}")
print(f"Total plans processed: {len(all_results)}")
print(f"\n{'Plan':<6} {'WSE':<10} {'Depth':<10} {'Inundation Area (acres)':<25}")
print("-" * 55)

for plan_num, results in sorted(all_results.items()):
    wse_count = len(results.get('wse', []))
    depth_count = len(results.get('depth', []))

    if plan_num in inundation_polygons:
        area = f"{inundation_polygons[plan_num]['area_acres']:,.1f}"
    else:
        area = "N/A"

    print(f"{plan_num:<6} {wse_count:<10} {depth_count:<10} {area:<25}")

print("\n" + "=" * 70)
print(f"\nOutput files located in: {notebook_display_path(OUTPUT_BASE)}")
Text Only
======================================================================
FLOODPLAIN MAPPING SUMMARY
======================================================================

Project: BaldEagleDamBrk
Profile: Max
Total plans processed: 1

Plan   WSE        Depth      Inundation Area (acres)  
-------------------------------------------------------
06     1          1          31,558.5

======================================================================

Output files located in: example_projects/BaldEagleCrkMulti2D_601/FloodplainMaps

Step 8: Visualize Results

Python
# Visualize a sample depth raster and inundation boundary
sample_plan = list(all_results.keys())[0] if all_results else None

if sample_plan and 'depth' in all_results[sample_plan] and all_results[sample_plan]['depth']:
    depth_tif = all_results[sample_plan]['depth'][0]

    fig, axes = plt.subplots(1, 2, figsize=(16, 8))

    # Plot depth raster
    with rasterio.open(depth_tif) as src:
        show(src, ax=axes[0], cmap='Blues', title=f'Max Depth - Plan {sample_plan}')
    axes[0].set_xlabel('Easting')
    axes[0].set_ylabel('Northing')

    # Plot inundation boundary
    if sample_plan in inundation_polygons:
        gdf = inundation_polygons[sample_plan]['gdf']
        gdf.plot(ax=axes[1], facecolor='lightblue', edgecolor='darkblue', linewidth=1)
        axes[1].set_title(f'Inundation Boundary - Plan {sample_plan}')
        axes[1].set_xlabel('Easting')
        axes[1].set_ylabel('Northing')
        axes[1].set_aspect('equal')

    plt.tight_layout()
    plt.show()
else:
    print("No results available for visualization")

png

Optional: Export All Inundation Boundaries to GeoPackage

Python
# Combine all inundation boundaries into a single GeoPackage
if inundation_polygons:
    all_gdfs = []

    for plan_num, data in inundation_polygons.items():
        gdf = data['gdf'].copy()
        gdf['plan_number'] = plan_num
        all_gdfs.append(gdf)

    combined_gdf = gpd.GeoDataFrame(pd.concat(all_gdfs, ignore_index=True))

    output_gpkg = OUTPUT_BASE / "all_inundation_boundaries.gpkg"
    combined_gdf.to_file(output_gpkg, driver='GPKG')

    print(f"Combined inundation boundaries saved to: {notebook_display_path(output_gpkg)}")
    print(f"Total plans: {len(inundation_polygons)}")
else:
    print("No inundation polygons to export")
Text Only
Combined inundation boundaries saved to: example_projects/BaldEagleCrkMulti2D_601/FloodplainMaps/all_inundation_boundaries.gpkg
Total plans: 1

Technical Notes

RasMapper stored maps

RasProcess.store_maps() is the ras-commander entry point for headless stored map generation. It prepares the plan's RAS Mapper configuration, invokes the native stored-map workflow, and returns the generated WSE, depth, and velocity raster paths.

Inundation Boundary Generation

The inundation boundary is created by: 1. Reading the depth raster 2. Creating a binary mask where depth > threshold 3. Converting the mask to polygon features using rasterio.features.shapes 4. Merging all polygons into a single boundary using shapely.ops.unary_union

Performance

  • Headless RasMapper stored map generation is the fastest method for generating maps (8-10 seconds per plan)
  • Polygon conversion adds ~1-2 seconds per depth raster
  • For large projects with many plans, consider running in batches
CLB Engineering Corporation  ·  LLM Forward Engineering
RAS Commander Documentation is a free and open-source project maintained by CLB Engineering Corporation. For agencies and firms modernizing H&H workflows with LLM Forward approaches, contact CLB to partner with the engineers who wrote the automation.