2D Spatial Result Queries with HdfResultsQuery

Query water surface elevation, depth, and velocity at arbitrary (x,y) coordinates, extract profiles along transects, compute flood extent with engineering filters, and generate domain-wide statistics -- all using scipy KDTree spatial indexing.

Text Only

from ras_commander import *
from ras_commander.hdf import HdfResultsQuery
from pathlib import Path

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

plt.style.use("seaborn-v0_8-whitegrid")
pd.set_option("display.max_columns", 20)
pd.set_option("display.width", 120)

Setup and Project Extraction

Use the reproducible BaldEagleCrkMulti2D example project, initialize it with ras-commander, and compute plan 06 before issuing any spatial queries. The notebook uses plan numbers directly so the decorators resolve the plan HDF and companion geometry HDF for us.

Text Only

PROJECT_NAME = "BaldEagleCrkMulti2D"
PLAN_NUMBER = "06"
PROJECT_SUFFIX = "415"
RAS_VERSION_OR_PATH = "7.0"  # Replace with a full Ras.exe path if needed.

try:
    from _notebook_prereqs import ensure_plan_result_hdf, get_or_extract_example_project
except ModuleNotFoundError:
    import sys

    examples_dir = Path.cwd() if Path.cwd().name == "examples" else Path.cwd() / "examples"
    if examples_dir.exists() and str(examples_dir) not in sys.path:
        sys.path.insert(0, str(examples_dir))
    from _notebook_prereqs import ensure_plan_result_hdf, get_or_extract_example_project

project_path = get_or_extract_example_project(PROJECT_NAME, suffix=PROJECT_SUFFIX)
init_ras_project(project_path, RAS_VERSION_OR_PATH)

# Ensure the required results HDF is present before analysis cells run.
_plan_hdf_check = ensure_plan_result_hdf(
    PLAN_NUMBER,
    ras_object=ras,
    num_cores=2,
    clear_geompre=False,
)

plan_columns = [
    col for col in ["plan_number", "plan_title", "Geom File", "HDF_Results_Path"]
    if col in ras.plan_df.columns
]
plan_overview = ras.plan_df.loc[:, plan_columns].copy()
plan_overview["plan_number"] = plan_overview["plan_number"].astype(str).str.zfill(2)

plan_hdf = Path(_plan_hdf_check)

print(f"Project folder: {project_path}")
print(f"Plan HDF: {plan_hdf.name}")
plan_overview.head()

Build Reproducible Query Locations

Rather than hard-coding coordinates, pull maximum-depth and maximum-WSE envelopes, keep wet cells, and reuse those cell-center coordinates for point, batch, and profile queries. This keeps the notebook reproducible even if the extracted project folder changes.

Text Only

# get_mesh_max_ws returns cell-center GeoDataFrame with max WSE values
max_wse_gdf = HdfResultsMesh.get_mesh_max_ws(PLAN_NUMBER)

# Compute max depth from WSE and cell minimum elevation (from geometry HDF)
import h5py

geom_number = ras.plan_df.loc[
    ras.plan_df["plan_number"].astype(str).str.zfill(2) == PLAN_NUMBER,
    "Geom File",
].iloc[0]
geom_hdf_path = ras.project_folder / f"{ras.project_name}.g{str(geom_number).zfill(2)}.hdf"

cell_min_elev = {}
with h5py.File(geom_hdf_path, "r") as gf:
    base = "Geometry/2D Flow Areas"
    for mesh_name in gf[base]:
        elev_path = f"{base}/{mesh_name}/Cells Minimum Elevation"
        if elev_path in gf:
            cell_min_elev[mesh_name] = gf[elev_path][()]

# Add minimum elevation to GeoDataFrame and compute depth
elev_series = []
for _, row in max_wse_gdf.iterrows():
    mesh = row["mesh_name"]
    cid = int(row["cell_id"])
    if mesh in cell_min_elev and cid < len(cell_min_elev[mesh]):
        elev_series.append(float(cell_min_elev[mesh][cid]))
    else:
        elev_series.append(float("nan"))

max_wse_gdf["cell_min_elevation"] = elev_series
max_wse_gdf["maximum_depth"] = (
    max_wse_gdf["maximum_water_surface"] - max_wse_gdf["cell_min_elevation"]
).clip(lower=0)

analysis_cells = max_wse_gdf.assign(
    x=lambda df: df.geometry.x,
    y=lambda df: df.geometry.y,
)

wet_cells = analysis_cells.loc[analysis_cells["maximum_depth"] > 0.5].copy()
if wet_cells.empty:
    wet_cells = analysis_cells.nlargest(12, "maximum_depth").copy()

wet_cells = wet_cells.sort_values(["mesh_name", "x"]).reset_index(drop=True)
sample_count = min(6, len(wet_cells))
sample_index = np.linspace(0, len(wet_cells) - 1, sample_count, dtype=int)
sample_points = wet_cells.iloc[np.unique(sample_index)].copy().reset_index(drop=True)
sample_points["point_id"] = [f"P{i+1}" for i in range(len(sample_points))]

mesh_for_profile = sample_points["mesh_name"].mode().iloc[0]
profile_pool = wet_cells.loc[wet_cells["mesh_name"] == mesh_for_profile].reset_index(drop=True)
if len(profile_pool) < 2:
    profile_pool = analysis_cells.reset_index(drop=True)
if len(profile_pool) < 2:
    raise ValueError("Need at least two cell centers to build a profile transect")

# Use a local transect between neighboring wet cell centers to avoid sampling outside the mesh.
anchor_idx = len(profile_pool) // 2
anchor = profile_pool.iloc[anchor_idx]
neighbor_distances = (profile_pool["x"] - anchor["x"]) ** 2 + (profile_pool["y"] - anchor["y"]) ** 2
neighbor_distances.iloc[anchor_idx] = np.inf
neighbor_idx = int(neighbor_distances.idxmin())
profile_start = anchor
profile_end = profile_pool.loc[neighbor_idx]

sample_points[
    ["point_id", "mesh_name", "cell_id", "x", "y", "maximum_depth", "maximum_water_surface"]
].head()

Point Query

query_point() snaps an arbitrary coordinate to the nearest 2D cell center and returns a compact dictionary with the value, mesh name, cell id, and snap distance. Using time_index="max" makes the first example land on the peak envelope instead of the final saved timestep.

Text Only

representative_point = sample_points.iloc[len(sample_points) // 2]
print(
    f"Using {representative_point['point_id']} on {representative_point['mesh_name']} "
    f"at ({representative_point['x']:.2f}, {representative_point['y']:.2f})"
)

point_result = HdfResultsQuery.query_point(
    PLAN_NUMBER,
    x=float(representative_point["x"]),
    y=float(representative_point["y"]),
    variable="wse",
    time_index="max",
)

point_result

Batch Point Query

query_points() handles many coordinates in one call and returns a tidy DataFrame. Here we query several wet-cell centers and use the maximum water-surface envelope so every sample is on the same peak-result basis.

Text Only

batch_results = (
    HdfResultsQuery.query_points(
        PLAN_NUMBER,
        sample_points[["x", "y"]].to_numpy(),
        variable="wse",
        time_index="max",
    )
    .rename(columns={"value": "wse_max"})
    .assign(point_id=sample_points["point_id"].values)
)

batch_results = batch_results[
    ["point_id", "x", "y", "mesh_name", "cell_id", "distance", "wse_max"]
]
batch_results.head()

Text Only

fig, ax = plt.subplots(figsize=(9, 6))
ax.scatter(
    wet_cells["x"],
    wet_cells["y"],
    s=8,
    color="#d9d9d9",
    alpha=0.35,
    label="Wet cell centers (> 0.5 depth)",
)
scatter = ax.scatter(
    batch_results["x"],
    batch_results["y"],
    c=batch_results["wse_max"],
    s=120,
    cmap="viridis",
    edgecolor="black",
    linewidth=0.8,
    label="Query points",
)

for row in batch_results.itertuples():
    ax.annotate(row.point_id, (row.x, row.y), xytext=(5, 5), textcoords="offset points")

ax.set_title("Batch spatial queries colored by maximum WSE")
ax.set_xlabel("X coordinate")
ax.set_ylabel("Y coordinate")
ax.legend(loc="best")
fig.colorbar(scatter, ax=ax, label="Maximum WSE")
plt.show()

Profile Extraction

query_profile() samples evenly spaced points along a transect. The start and end points below are taken from wet cells in the dominant mesh so the section crosses the active flow area without hard-coded coordinates.

Text Only

profile_df = HdfResultsQuery.query_profile(
    PLAN_NUMBER,
    x1=float(profile_start["x"]),
    y1=float(profile_start["y"]),
    x2=float(profile_end["x"]),
    y2=float(profile_end["y"]),
    variable="wse",
    n_points=75,
    time_index="max",
)

display(profile_df.head())

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(profile_df["station"], profile_df["value"], color="#005f73", linewidth=2.5)
ax.set_title(f"Maximum WSE along transect in {mesh_for_profile}")
ax.set_xlabel("Station along transect")
ax.set_ylabel("Water surface elevation")
ax.grid(True, alpha=0.3)
plt.show()

Flood Extent Analysis

flood_extent() can do more than a simple depth screen. The three cases below show raw depth, a shallow ponding exclusion intended for rain-on-grid screening, and a hazard-oriented filter using depth x velocity >= 3.0. For this example project the ponding filter is illustrative rather than a calibrated rainfall-depth threshold.

Text Only

raw_extent = HdfResultsQuery.flood_extent(
    PLAN_NUMBER,
    depth_threshold=0.10,
    time_index="max",
)

filtered_extent = HdfResultsQuery.flood_extent(
    PLAN_NUMBER,
    depth_threshold=0.10,
    precip_depth_threshold=0.25,
    ponding_velocity_max=0.10,
    time_index="max",
)

hazard_extent = HdfResultsQuery.flood_extent(
    PLAN_NUMBER,
    depth_threshold=0.50,
    dv_threshold=3.0,
    precip_depth_threshold=0.25,
    ponding_velocity_max=0.10,
    time_index="max",
)

extent_comparison = pd.DataFrame(
    [
        {"scenario": "Raw depth > 0.10", **raw_extent},
        {"scenario": "Ponding excluded", **filtered_extent},
        {"scenario": "Hazard DV >= 3.0", **hazard_extent},
    ]
)

# Display only columns that exist (some filters add extra keys)
display_cols = ["scenario", "wet_cells", "wet_fraction", "wet_area_acres",
                "max_depth", "mean_wet_depth"]
for optional_col in ["ponding_excluded_cells", "max_dv", "filters_applied"]:
    if optional_col in extent_comparison.columns:
        display_cols.append(optional_col)

display(extent_comparison[display_cols])

fig, ax = plt.subplots(figsize=(9, 4))
bars = ax.bar(
    extent_comparison["scenario"],
    extent_comparison["wet_area_acres"],
    color=["#94d2bd", "#ee9b00", "#bb3e03"],
)
ax.set_title("Flood extent shrinks as filters become more selective")
ax.set_ylabel("Wet area (acres)")
for bar in bars:
    height = bar.get_height()
    ax.text(
        bar.get_x() + bar.get_width() / 2,
        height,
        f"{height:,.1f}",
        ha="center",
        va="bottom",
    )
plt.xticks(rotation=10)
plt.show()

Domain Statistics

result_statistics() collapses an entire field into a compact summary. Comparing wet_only=True and wet_only=False makes it obvious how dry cells pull the distribution toward zero.

Text Only

stats_all = HdfResultsQuery.result_statistics(
    PLAN_NUMBER,
    variable="depth",
    time_index="max",
    wet_only=False,
)

stats_wet = HdfResultsQuery.result_statistics(
    PLAN_NUMBER,
    variable="depth",
    time_index="max",
    wet_only=True,
    depth_threshold=0.10,
)

stats_df = pd.DataFrame(
    [stats_all, stats_wet],
    index=["All cells", "Wet cells only"],
).T
display(stats_df.loc[["count", "min", "p10", "p25", "mean", "median", "p75", "p90", "p99", "max", "std"]])

# Clip the extreme tail so the histogram stays readable.
plot_depths = analysis_cells["maximum_depth"].clip(
    upper=analysis_cells["maximum_depth"].quantile(0.99)
)

fig, ax = plt.subplots(figsize=(9, 4))
ax.hist(plot_depths, bins=40, color="#0a9396", alpha=0.85)
ax.axvline(stats_wet["mean"], color="#ae2012", linestyle="--", linewidth=2, label="Wet-cell mean")
ax.set_title("Maximum depth distribution across 2D cells")
ax.set_xlabel("Maximum depth")
ax.set_ylabel("Cell count")
ax.legend()
plt.show()

Max Envelope vs Time Step

The same coordinate can answer two different questions. time_index=-1 asks for the model state at the final saved timestep, while time_index="max" asks for the peak envelope at that location over the full run. Comparing both at the same sample points highlights recession between the peak and the last output.

Text Only

last_step_results = HdfResultsQuery.query_points(
    PLAN_NUMBER,
    sample_points[["x", "y"]].to_numpy(),
    variable="wse",
    time_index=-1,
).rename(columns={"value": "wse_last"})

max_envelope_results = HdfResultsQuery.query_points(
    PLAN_NUMBER,
    sample_points[["x", "y"]].to_numpy(),
    variable="wse",
    time_index="max",
).rename(columns={"value": "wse_max"})

envelope_comparison = sample_points[["point_id", "mesh_name", "cell_id", "x", "y"]].copy()
envelope_comparison["wse_last"] = last_step_results["wse_last"].values
envelope_comparison["wse_max"] = max_envelope_results["wse_max"].values
envelope_comparison["difference"] = envelope_comparison["wse_max"] - envelope_comparison["wse_last"]

display(envelope_comparison.head())

fig, ax = plt.subplots(figsize=(9, 4))
positions = np.arange(len(envelope_comparison))
width = 0.38

ax.bar(
    positions - width / 2,
    envelope_comparison["wse_last"],
    width=width,
    label="Last timestep",
    color="#6c757d",
)
ax.bar(
    positions + width / 2,
    envelope_comparison["wse_max"],
    width=width,
    label="Max envelope",
    color="#1d3557",
)
ax.set_xticks(positions)
ax.set_xticklabels(envelope_comparison["point_id"])
ax.set_ylabel("Water surface elevation")
ax.set_title("Peak envelope vs final-state WSE at the same coordinates")
ax.legend()
plt.show()

Key Takeaways

• HdfResultsQuery.query_point() and query_points() turn arbitrary coordinates into model values without pre-built extraction features.
• query_profile() creates quick transect plots directly from the 2D mesh.
• flood_extent() supports depth-only screening, hazard-style depth x velocity filters, and shallow-ponding exclusions.
• result_statistics() is useful for quick QA/QC, especially when comparing wet-only statistics to the full domain.
• time_index="max" returns peak envelopes, while time_index=-1 returns the final saved model state.