Cloud-Native Results Export with ras2cng¶
Export HEC-RAS simulation results to cloud-native GeoParquet, build interactive flood depth maps, and generate PMTiles for web deployment using ras2cng.
USE_LOCAL_SOURCE = False
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(f"LOCAL SOURCE MODE: Loading from {local_path}/ras_commander")
else:
print("PIP PACKAGE MODE: Loading installed ras-commander")
from ras_commander import RasExamples, init_ras_project, RasCmdr, ras
from pathlib import Path
import geopandas as gpd
import ras_commander
print(f"ras-commander: {ras_commander.__file__}")
import ras2cng
print(f"ras2cng: {ras2cng.__file__}")
Prerequisites¶
- ras-commander installed:
pip install ras-commander - ras2cng installed:
pip install "ras2cng[duckdb]" - HEC-RAS installed: Required to run simulation (results export needs
.p##.hdffiles) - Optional:
tippecanoe+pmtilesCLIs for PMTiles generation
What You'll Learn¶
- Run a HEC-RAS plan and export results to GeoParquet
- Join results to mesh cell polygons for thematic mapping
- Query flood depth statistics with DuckDB SQL
- Display interactive flood depth maps in the notebook
- Archive a complete project (geometry + results + terrain)
- Generate PMTiles for serverless flood map deployment
Related Notebooks¶
- 960_cloud_native_geometry_export.ipynb — Export geometry to GeoParquet
- 410_2d_hdf_data_extraction.ipynb — Extract 2D results with ras-commander
- 110_single_plan_execution.ipynb — Execute plans with ras-commander
Parameters¶
PROJECT_NAME = "BaldEagleCrkMulti2D"
RAS_VERSION = "6.6"
OUTPUT_DIR = Path("out/961_cloud_native_results")
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
Example 1: Results Export with Polygon Join¶
Extract a project, run a simulation, then export the results joined to mesh cell polygons for thematic flood mapping.
Extract and Initialize Project¶
project_path = RasExamples.extract_project(PROJECT_NAME)
ras_obj = init_ras_project(project_path, RAS_VERSION)
print(f"Project: {project_path}")
print(f"Plans: {len(ras.plan_df)}")
print(f"Geometries: {len(ras.geom_df)}")
ras.plan_df[['plan_number', 'Plan Title', 'Geom File']].head()
Run Simulation¶
Execute the first plan to generate results HDF files. Smart skip will avoid re-running if results are already current.
plan_number = ras.plan_df['plan_number'].iloc[0]
print(f"Running plan {plan_number}...")
try:
RasCmdr.compute_plan(plan_number)
print(f"Plan {plan_number} completed.")
except Exception as e:
print(f"Execution failed: {e}")
print("Continuing with any existing HDF files...")
plan_hdfs = sorted(project_path.glob("*.p??.hdf"))
geom_hdfs = sorted(project_path.glob("*.g??.hdf"))
print(f"Plan HDF files: {[f.name for f in plan_hdfs]}")
print(f"Geometry HDF files: {[f.name for f in geom_hdfs]}")
if not plan_hdfs:
print("\nNo plan HDF files found. Run HEC-RAS first to generate results.")
plan_hdf = None
else:
plan_hdf = plan_hdfs[0]
print(f"\nUsing plan HDF: {plan_hdf.name} ({plan_hdf.stat().st_size / 1e6:.1f} MB)")
# Match geometry HDF to the plan being used
plan_row = ras.plan_df[ras.plan_df['plan_number'] == plan_number].iloc[0]
geom_num = plan_row.get('Geom File', '01')
geom_hdf_match = project_path / f"{ras.project_name}.g{geom_num}.hdf"
if geom_hdf_match.exists():
geom_hdf = geom_hdf_match
else:
geom_hdf = geom_hdfs[0] if geom_hdfs else None
if geom_hdf:
print(f"Using geom HDF: {geom_hdf.name} (matched to plan {plan_number})")
Export Mesh Geometry¶
First, export the mesh cell polygons. These serve as the base layer for joining results.
from ras2cng import export_geometry_layers
mesh_out = OUTPUT_DIR / "mesh_cells.parquet"
if geom_hdf:
export_geometry_layers(geom_hdf, mesh_out, layer="mesh_cells")
gdf_mesh = gpd.read_parquet(mesh_out)
print(f"Exported {len(gdf_mesh)} mesh cells")
elif plan_hdf:
export_geometry_layers(plan_hdf, mesh_out, layer="mesh_cells")
gdf_mesh = gpd.read_parquet(mesh_out)
print(f"Exported {len(gdf_mesh)} mesh cells (from plan HDF)")
else:
print("No HDF files available for geometry export")
gdf_mesh = None
Export Results Joined to Polygons¶
Export Maximum Water Surface results joined to mesh cell polygons. Each row becomes a polygon with the result value as an attribute — ready for thematic mapping.
from ras2cng import export_results_layer, export_all_variables
gdf_wse = None
if plan_hdf and mesh_out.exists():
wse_out = OUTPUT_DIR / "maximum_wse.parquet"
try:
export_results_layer(
plan_hdf,
wse_out,
variable="Maximum Water Surface",
geom_file=mesh_out,
)
gdf_wse = gpd.read_parquet(wse_out)
print(f"Exported {len(gdf_wse)} cells with Maximum Water Surface")
print(f"WSE range: {gdf_wse['maximum_water_surface'].min():.2f} - {gdf_wse['maximum_water_surface'].max():.2f}")
gdf_wse.head()
except Exception as e:
print(f"Results export failed: {e}")
print("This may occur if the plan has no 2D summary results.")
else:
print("No plan HDF or mesh geometry available for results export")
DuckDB Analytics on Results¶
Query water surface elevation statistics by mesh area using DuckDB.
from ras2cng import query_parquet
if gdf_wse is not None:
df_flood_stats = query_parquet(
wse_out,
"""
SELECT
mesh_name,
COUNT(*) AS total_cells,
ROUND(MIN(maximum_water_surface), 2) AS min_wse,
ROUND(MAX(maximum_water_surface), 2) AS max_wse,
ROUND(AVG(maximum_water_surface), 2) AS avg_wse
FROM _
GROUP BY mesh_name
ORDER BY max_wse DESC
"""
)
print("Water surface elevation statistics by mesh area:")
df_flood_stats
else:
print("No results available for DuckDB analytics.")
if gdf_wse is not None:
df_high_wse = query_parquet(
wse_out,
"""
SELECT mesh_name, cell_id, ROUND(maximum_water_surface, 2) AS wse_ft
FROM _
ORDER BY maximum_water_surface DESC
LIMIT 20
"""
)
print("Top 20 highest WSE cells:")
df_high_wse
else:
print("No results available.")
Interactive Water Surface Elevation Map¶
Display the WSE results as a choropleth map. Cells are colored by maximum water surface elevation.
if gdf_wse is not None:
try:
gdf_wse_wgs84 = gdf_wse.to_crs("EPSG:4326")
print(f"Displaying {len(gdf_wse)} cells with WSE results")
gdf_wse_wgs84.explore(
column="maximum_water_surface",
cmap="Blues",
tooltip=["mesh_name", "cell_id", "maximum_water_surface"],
style_kwds={"weight": 0.3, "fillOpacity": 0.6},
legend_kwds={"caption": "Maximum Water Surface (ft)"},
tiles="CartoDB positron",
)
except ImportError as e:
print(f"Interactive map requires folium and mapclassify: {e}")
print("Install with: pip install folium mapclassify")
print(f"\nWSE range: {gdf_wse['maximum_water_surface'].min():.2f} - {gdf_wse['maximum_water_surface'].max():.2f}")
else:
print("No results available for mapping.")
Example 2: Full Project Archive¶
Archive the entire project — all geometry files, plan results, and terrain — to a consolidated GeoParquet archive with a manifest.json catalog.
from ras2cng import archive_project, inspect_project
info = inspect_project(project_path)
print(f"Project: {info.name}")
print(f"Geometry files: {len(info.geom_files)}")
print(f"Plan files: {len(info.plan_files)}")
print(f"Has terrain: {bool(info.terrain_files)}")
archive_dir = OUTPUT_DIR / "archive"
archive_dir.mkdir(parents=True, exist_ok=True)
manifest = archive_project(
project_path,
archive_dir,
include_results=bool(plan_hdfs),
include_terrain=False,
)
print(f"\nArchive complete:")
print(f" Geometry configs: {len(manifest.geometry)}")
if hasattr(manifest, 'plans'):
print(f" Plan results: {len(manifest.plans)}")
print(f"\nArchive files:")
for f in sorted(archive_dir.glob("*")):
print(f" {f.name:40s} {f.stat().st_size / 1e3:>8.1f} KB")
archive_parquets = sorted(archive_dir.glob("*.parquet"))
if archive_parquets:
for pq in archive_parquets:
try:
gdf = gpd.read_parquet(pq)
if 'layer' in gdf.columns:
layers = gdf['layer'].unique()
print(f"{pq.name}: {len(gdf)} features, layers: {list(layers)}")
elif '_table' in gdf.columns:
tables = gdf['_table'].unique()
print(f"{pq.name}: {len(gdf)} rows, tables: {list(tables)}")
else:
print(f"{pq.name}: {len(gdf)} features")
except Exception as e:
import pandas as pd
try:
df = pd.read_parquet(pq)
print(f"{pq.name}: {len(df)} rows (non-spatial)")
except Exception:
print(f"{pq.name}: could not read ({e})")
PMTiles for Flood Results (Optional)¶
Generate vector PMTiles from the flood depth results for web map deployment.
import shutil
has_tools = bool(shutil.which("tippecanoe") and shutil.which("pmtiles"))
if has_tools and gdf_wse is not None:
from ras2cng import generate_pmtiles_from_input
pmtiles_out = OUTPUT_DIR / "flood_wse.pmtiles"
generate_pmtiles_from_input(
wse_out,
pmtiles_out,
layer_name="flood_wse",
min_zoom=8,
max_zoom=14,
)
print(f"PMTiles: {pmtiles_out.name} ({pmtiles_out.stat().st_size / 1e6:.1f} MB)")
print("\nDeploy to S3/R2/GitHub Pages for serverless flood maps.")
elif not has_tools:
print("PMTiles CLIs not found. Install: conda install -c conda-forge tippecanoe pmtiles")
print("\nCLI equivalent:")
print(f" ras2cng pmtiles {wse_out if gdf_wse is not None else 'results.parquet'} flood_wse.pmtiles --layer flood --min-zoom 8 --max-zoom 14")
else:
print("No results available for PMTiles generation.")
Flood Depth Web Map with PMTiles¶
Style the flood depth PMTiles with a color ramp in MapLibre GL JS:
<script src="https://unpkg.com/maplibre-gl/dist/maplibre-gl.js"></script>
<script src="https://unpkg.com/pmtiles/dist/pmtiles.js"></script>
<script>
let protocol = new pmtiles.Protocol();
maplibregl.addProtocol("pmtiles", protocol.tile);
const map = new maplibregl.Map({
container: "map",
style: {
version: 8,
sources: {
flood: {
type: "vector",
url: "pmtiles://./flood_depth.pmtiles",
},
},
layers: [
{
id: "flood-fill",
source: "flood",
"source-layer": "flood_depth",
type: "fill",
paint: {
"fill-color": [
"interpolate", ["linear"], ["get", "maximum_depth"],
0, "#ffffcc",
1, "#feb24c",
3, "#f03b20",
10, "#bd0026"
],
"fill-opacity": 0.7,
},
},
],
},
});
</script>
Output Summary¶
print("Output files:")
for f in sorted(OUTPUT_DIR.rglob("*")):
if f.is_file():
rel = f.relative_to(OUTPUT_DIR)
print(f" {str(rel):45s} {f.stat().st_size / 1e3:>8.1f} KB")
Key Takeaways¶
- Results join to geometry polygons — each mesh cell becomes a polygon with flood depth as an attribute
- DuckDB enables SQL analytics — aggregate by mesh area, filter deep cells, compute statistics
- GeoPandas
.explore()gives instant flood maps — choropleth visualization in the notebook - Full project archive creates a manifest-tracked GeoParquet bundle — one file per geometry, one per plan
- PMTiles enable serverless flood maps — deploy to static hosting, no tile server needed
CLI Quick Reference¶
# Full project archive (geometry + results + terrain)
ras2cng archive path/to/project/ ./archive/ --results --terrain
# Export results joined to polygon geometry
ras2cng results model.p01.hdf max_depth.parquet \
--geometry mesh_cells.parquet --var "Maximum Depth"
# Query flood statistics
ras2cng query max_depth.parquet \
"SELECT mesh_name, MAX(maximum_depth) FROM _ GROUP BY mesh_name"
# Generate flood depth PMTiles
ras2cng pmtiles max_depth.parquet flood.pmtiles --layer flood
Next Steps¶
- ras2cng documentation — Complete API and CLI reference
- ras2cng GitHub — Source code and issues
- 960_cloud_native_geometry_export.ipynb — Geometry-only export (no HEC-RAS required)