[13pt] PR: Modifying OSM bridge scripts

data/bridges/pull_osm_bridges.py

@@ -3,14 +3,16 @@
 import logging
 import os
 import traceback
+import warnings
 from concurrent.futures import ProcessPoolExecutor, as_completed
 from pathlib import Path
 
 import geopandas as gpd
 import osmnx as ox
 import pandas as pd
 import pyproj
-from shapely.geometry import shape
+from networkx import Graph, connected_components
+from shapely.geometry import LineString, shape
 
 
 CRS = "epsg:5070"
@@ -19,6 +21,32 @@
 # Save all OSM bridge features by HUC8 to a specified folder location.
 # Bridges will have point geometry converted to linestrings if needed
 #
+# Dissolve touching lines
+def find_touching_groups(gdf):
+    # Create a graph
+    graph = Graph()
+
+    # Add nodes for each geometry
+    graph.add_nodes_from(gdf.index)
+
+    # Create spatial index for efficient querying
+    spatial_index = gdf.sindex
+
+    # For each geometry, find touching geometries and add edges to the graph
+    for idx, geometry in gdf.iterrows():
+        possible_matches_index = list(spatial_index.intersection(geometry['geometry'].bounds))
+        possible_matches = gdf.iloc[possible_matches_index]
+        precise_matches = possible_matches[possible_matches.intersects(geometry['geometry'])]
+
+        for match_idx in precise_matches.index:
+            if match_idx != idx:
+                graph.add_edge(idx, match_idx)
+
+    # Find connected components
+    groups = list(connected_components(graph))
+    return groups
+
+
 def pull_osm_features_by_huc(huc_bridge_file, huc_num, huc_geom):
     """
     Returns: The huc number but only if it failed, so we can make a master list of failed HUCs.
@@ -40,13 +68,47 @@ def pull_osm_features_by_huc(huc_bridge_file, huc_num, huc_geom):
             logging.info(f"osmnx pull for {huc_num} came back with no records")
             return huc_num
 
+        # Create bridge_type column
+        # Check if 'highway' column exists
+        if 'highway' not in gdf.columns:
+            gdf['highway'] = None
+
+        # Check if 'railway' column exists
+        if 'railway' not in gdf.columns:
+            gdf['railway'] = None
+
+        # Create the bridge_type column by combining above information
+        gdf['bridge_type'] = gdf.apply(
+            lambda row: (
+                f"highway-{row['highway']}" if pd.notna(row['highway']) else f"railway-{row['railway']}"
+            ),
+            axis=1,
+        )
+        gdf.reset_index(inplace=True)
+        # Remove abandoned bridges
+        gdf = gdf[gdf['bridge'] != 'abandoned']
+
         cols_to_drop = []
         for col in gdf.columns:
             if any(isinstance(val, list) for val in gdf[col]):
                 cols_to_drop.append(col)
 
         # This a common and know duplicate column name (and others)
-        bad_column_names = ["atv", "fixme", "FIXME"]
+        bad_column_names = [
+            "atv",
+            "fixme",
+            "FIXME",
+            "NYSDOT_ref",
+            "REF",
+            "fid",
+            "fixme:maxspeed",
+            "LAYER",
+            "unsigned_ref",
+            "Fut_Ref",
+            "Ref",
+            "FIXME:ref",
+            "id",
+        ]
         for bad_cn in bad_column_names:
             if bad_cn in gdf.columns:
                 cols_to_drop.append(bad_cn)
@@ -57,7 +119,37 @@ def pull_osm_features_by_huc(huc_bridge_file, huc_num, huc_geom):
         gdf1 = gdf[gdf.geometry.apply(lambda x: x.geom_type == 'LineString')]
 
         gdf1 = gdf1.to_crs(CRS)
-        gdf1.to_file(huc_bridge_file, driver="GPKG")
+
+        # Perform dissolve touching lines
+        buffered = gdf1.copy()
+        buffered['geometry'] = buffered['geometry'].buffer(0.0001)
+        # Find groups of touching geometries
+        touching_groups = find_touching_groups(buffered)
+
+        # Dissolve each group separately
+        warnings.filterwarnings('ignore')
+        dissolved_groups = []
+        for group in touching_groups:
+            group_gdf = buffered.loc[list(group)]
+            if not group_gdf.empty:
+                dissolved_group = group_gdf.dissolve()
+                single_part_group = dissolved_group.explode(index_parts=False)
+                dissolved_groups.append(single_part_group)
+
+        # Combine dissolved groups and reconstruct GeoDataFrame
+        if dissolved_groups:
+            dissolved_gdf = pd.concat(dissolved_groups, ignore_index=True)
+            final_gdf = gpd.GeoDataFrame(dissolved_gdf, crs=buffered.crs)
+        else:
+            final_gdf = buffered.copy()
+
+        # Polygon to linestring
+        final_gdf['geometry'] = final_gdf['geometry'].apply(
+            lambda geom: LineString(geom.exterior.coords) if geom.geom_type == 'Polygon' else geom
+        )
+        # Reconstruct the GeoDataFrame to remove fragmentation
+        final_gdf = final_gdf.copy()
+        final_gdf.to_file(huc_bridge_file, driver="GPKG")
 
         # returns the HUC but only if it failed so we can keep a list of failed HUCs
         return ""
@@ -95,7 +187,7 @@ def combine_huc_features(output_dir):
     section_time = dt.datetime.now(dt.timezone.utc)
     logging.info(f"  .. started: {section_time.strftime('%m/%d/%Y %H:%M:%S')}")
 
-    all_bridges_gdf = all_bridges_gdf_raw[['osmid', 'name', 'geometry']]
+    all_bridges_gdf = all_bridges_gdf_raw[['osmid', 'name', 'bridge_type', 'geometry']]
     all_bridges_gdf.to_file(osm_bridge_file, driver="GPKG")
 
     return
@@ -185,6 +277,15 @@ def process_osm_bridges(wbd_file, output_folder, number_of_jobs):
     # all huc8 processing must be completed before this function call
     combine_huc_features(output_folder)
 
+    # Clean up individual HUC8 files
+    logging.info('Deleting individual HUC8 files as a final cleanup step')
+    huc_files = Path(output_folder).glob('huc_*_osm_bridges.gpkg')
+    for huc_file in huc_files:
+        try:
+            os.remove(huc_file)
+        except Exception as e:
+            logging.info(f"Error deleting {huc_file}: {str(e)}")
+
     if len(failed_HUCs_list) > 0:
         logging.info("\n+++++++++++++++++++")
         logging.info("HUCs that failed to download from OSM correctly are:")
@@ -236,6 +337,7 @@ def __setup_logger(outputs_dir):
     logging.basicConfig(
         filename=log_file_path, level=logging.INFO, format='[%(asctime)s] %(message)s', datefmt='%H:%M:%S'
     )
+    logging.captureWarnings(True)
 
     # set up logging to console
     console = logging.StreamHandler()
@@ -260,7 +362,7 @@ def __setup_logger(outputs_dir):
         It should be run only as often as the user thinks OSM has had any important updates.
         - As written, the code will skip any HUC that there's already a file for.
         - Each HUC8's worth of OSM bridge features is saved out individually, then merged together
-        into one. The HUC8 files can be deleted if desired, as an added final cleanup step.
+        into one.
     '''
 
     parser = argparse.ArgumentParser(description='Acquires and saves Open Street Map bridge features')

docs/CHANGELOG.md

@@ -1,6 +1,18 @@
 All notable changes to this project will be documented in this file.
 We follow the [Semantic Versioning 2.0.0](http://semver.org/) format.
 
+## v4.5.11.0 - 2024-10-11 - [PR#1298](https://github.com/NOAA-OWP/inundation-mapping/pull/1298)
+
+This PR addresses four issues regarding OSM bridges. It dissolves touching bridge lines so each bridge has a single linestring. It also removes abandoned bridge from the dataset and it adds bridge type field to bridge centroids. As part of this PR, `max_hand_ft` and `max_discharge_cfs` columns are added to `osm_bridge_centroids.gkpg`.
+
+### Changes
+
+- `data/bridges/pull_osm_bridges.py`
+- `src/heal_bridges_osm.py`
+
+<br/><br/>
+
+
 ## v4.5.10.3 - 2024-10-11 - [PR#1306](https://github.com/NOAA-OWP/inundation-mapping/pull/1306)
 
 Extends outlet levelpath(s) outside HUC.

src/bash_variables.env

@@ -7,7 +7,7 @@ export                             ALASKA_CRS=EPSG:3338
 
 # NOTE: $inputsDir is defined in Dockerfile
 
-export                       pre_clip_huc_dir=${inputsDir}/pre_clip_huc8/20240917
+export                       pre_clip_huc_dir=${inputsDir}/pre_clip_huc8/20241002
 
 export                              input_DEM=${inputsDir}/dems/3dep_dems/10m_5070/20240916/hand_seamless_3dep_dem_10m_5070.vrt
 export                       input_DEM_Alaska=${inputsDir}/dems/3dep_dems/10m_South_Alaska/20240912/FIM_3dep_dem_South_Alaska_10m.vrt
@@ -32,7 +32,7 @@ export                          input_WBD_gdb=${inputsDir}/wbd/WBD_National_EPSG
 export                   input_WBD_gdb_Alaska=${inputsDir}/wbd/WBD_National_South_Alaska.gpkg
 export                 input_calib_points_dir=${inputsDir}/rating_curve/water_edge_database/calibration_points/
 export                        bathymetry_file=${inputsDir}/bathymetry/bathymetry_adjustment_data.gpkg
-export                            osm_bridges=${inputsDir}/osm/bridges/240426/osm_all_bridges.gpkg
+export                            osm_bridges=${inputsDir}/osm/bridges/241002/osm_all_bridges.gpkg
 
 # input file location with nwm feature_id and recurrence flow values
 export               bankfull_flows_file=${inputsDir}/rating_curve/bankfull_flows/nwm3_high_water_threshold_cms.csv

src/heal_bridges_osm.py

@@ -94,6 +94,12 @@ def flows_from_hydrotable(bridge_pnts, hydroTable):
         axis=1,
         result_type='expand',
     )
+    # Convert stages and dischrages to ft and cfs respectively
+    bridge_pnts['max_hand_ft'] = bridge_pnts['max_hand'] * 3.28084
+    bridge_pnts['max_hand_75_ft'] = bridge_pnts['max_hand_75'] * 3.28084
+    bridge_pnts['max_discharge_cfs'] = bridge_pnts['max_discharge'] * 35.3147
+    bridge_pnts['max_discharge_75_cfs'] = bridge_pnts['max_discharge75'] * 35.3147
+
     return bridge_pnts