[5pt] Incorporate roads into FIM

[AliForghani-NOAA]

This PR addresses the issue #1385 and includes the following enhancements:

• Ingests OSM roads as a new input data for FIM.

• Derives the threshold discharge required for a road inundation. To achieve higher spatial resolution, each road segment is further split at HydroID boundaries, and the minimum non-zero HAND value within each segment is extracted as the inundation threshold stage. The corresponding threshold discharge values are then interpolated from the HydroTables using these threshold stages.

• Develops a new tool to identify inundated roads for a given flood event and report the maximum flood depth for each road segment.

Note that pulling OSM data requires a new Python package called overpy. The updated pre-clipped dataset with new osm roads data has been prepared here: inputs/pre_clip_huc8/20250606/.

In-Depth Workflow Explanation

The workflow required developing three new scripts and a major update to src/aggregate_by_huc.py.

1. data/roads/pull_osm_roads.py (new script)
   This script downloads OSM road data. OSM data is structured using key-value pairs; for roads, the key we use is highway. In this script, we focus on five major road types (i.e., osm values): motorway, trunk, primary, secondary, and tertiary. For more details, see the OSM road tag documentation.

   Because OSM roads can be very long (e.g., >10 km), the script splits them using NWM catchment boundaries (found in data/inputs/pre_clip_huc8) to generate more localized and accurate FIMpact results. Each resulting road segment is assigned a new ID in the format osmid_catchmentid and is used as the input road dataset for the FIM pipeline. The updated pre-clipped dataset with new osm roads data has been prepared here: inputs/pre_clip_huc8/20250606/

2. src/process_roads_fimpact.py (new script)
   This script is run for each branch of an HUC using three inputs:

   • Road segments (which are already splitted based on NWM catchments)
   • The HAND raster
   • The HAND-generated HydroIDs gpkg

   HAND-generated HydroIDs offer additional spatial resolution, as a single road segment may intersect multiple HydroIDs (each identified by a unique feature_id). To account for this, the script splits road segments at HydroID boundaries and calculates the minimum HAND value (excluding zeros) within each segment to serve as the inundation threshold.

   We use rasterstats.zonal_stats() to extract HAND values; however, this function can return incorrect results when the input lines are jagged—often caused by overlaying roads with catchment polygons. To mitigate this, we first explode all multi-part geometries to ensure each segment is a single LineString. After running the zonal statistics, we group by osmid_catchid and HydroID, retaining only the segment with the minimum threshold HAND value within each HydroID to eliminate redundant exploded segments.

   Three new columns are added to the road dataset: threshold_hand, HydroID, and feature_id. The results are saved as osm_roads_fimpact_***.csv for each branch, where *** represents the branch number. Each CSV file contains one record per osmid_catchid within each HydroID, providing the minimum HAND value for that combination.

3. src/aggregate_by_huc.py (updated script)
   For each branch, the script retrieves the discharge value corresponding to each threshold_hand from the branch’s HydroTable (per HydroID) and assigns it as threshold_discharge. Any record with a threshold_hand value greater than 25m (the maximum stage listed in the HydroTables) is removed entirely. The outputs from all branches are combined into a single file: osm_roads_fimpact.csv.

4. tools/road_inundation.py (new script)
   This tool takes two inputs:

   • A FIM run directory (which includes the osm_roads_fimpact.csv file), and
   • A flow file.

   The script identifies road segments where the given flow (referred to as evaluated_discharge) exceeds the threshold discharge and flags them as inundated. It also looks up the stage corresponding to the evaluated_discharge (and call it evaluated_stage) and subtracts the evaluated_stage from the threshold_hand value to calculate the flood_depth.

   Records with negative flood depth are currently removed, as these may result from non-monotonic synthetic rating curves—most commonly observed in branch zero.

   Note that a single road segment may have multiple inundation records, originating from different branches or intersecting multiple HydroIDs. The code retains only the record with the maximum flood depth for each road segment.

   The figure below displays the output of the new tool, with inundated roads (with their flood depth) and non-inundated roads (black) overlaid on a FIM raster. Both results were generated from a common 50-year recurrence interval flow file for HUC 11070103.

Additions

• data/roads/pull_osm_roads.py
• src/process_roads_fimpact.py
• tools/road_inundation.py

Changes

• src/aggregate_by_huc.py -> Aggregates road FIMpact results by HUC
• src/delineate_hydros_and_produce_HAND.sh -> Calls the new src/process_roads_fimpact.py script
• src/bash_variables.env -> Updated to reference the new pre-clipped dataset
• fim_post_processing.sh -> Enables aggregating road FIMpact results by HUC

Testing

The new Docker image with the new overpy Python package was tested successfully on data/roads/pull_osm_roads.py to retrieve OSM road data for both CONUS and Alaska.

The road FIMpact functionality has been successfully tested on three HUCs ( 05030104, 11070103, and 12030203).

Deployment Plan (For developer use)

How does the changes affect the product?

• Code only?
• If applicable, has a deployment plan be created with the deployment person/team?
• Require new or adjusted data inputs? Does it have start, end and duration code (in UTC)?
• If new or updated data sets, has the FIM code been updated and tested with the new/adjusted data (subset is fine, but must be a subset of the new data)?
• Require new pre-clip set?
• Has new or updated python packages?

Issuer Checklist (For developer use)

You may update this checklist before and/or after creating the PR. If you're unsure about any of them, please ask, we're here to help! These items are what we are going to look for before merging your code.

• Informative and human-readable title, using the format: [_pt] PR: <description>
• Links are provided if this PR resolves an issue, or depends on another other PR
• If submitting a PR to the dev branch (the default branch), you have a descriptive Feature Branch name using the format: dev-<description-of-change> (e.g. dev-revise-levee-masking)
• Changes are limited to a single goal (no scope creep)
• The feature branch you're submitting as a PR is up to date (merged) with the latest dev branch
• pre-commit hooks were run locally
• Any change in functionality is tested
• New functions are documented (with a description, list of inputs, and expected output)
• Placeholder code is flagged / future todos are captured in comments
• CHANGELOG updated with template version number, e.g. 4.x.x.x
• Add yourself as an assignee in the PR as well as the FIM Technical Lead

Merge Checklist (For Technical Lead use only)

• Update CHANGELOG with latest version number and merge date
• Update the Citation.cff file to reflect the latest version number in the CHANGELOG
• If applicable, update README with major alterations

[ZahraGhahremani]

Ran fim_pipline and looked at the inundation results. I also tested road_inundation.py to look at the inundated roads.
HUC 11010003, 25 yr flow.

[RobHanna-NOAA]

When I look at the alaska_osm_roads.gpkg and conus_osm_roads.gpkg, the layer names could be changed. Weird layer name but maybe that was desired. Just checking. :)

[RobHanna-NOAA]

I only tested the pull_osm_roads early on. Sounds like we are all good on all fronts.