diff --git a/examples/dflowfm_local/dflowfm/structures.ini b/examples/dflowfm_local/dflowfm/structures.ini index 4065f58a..448ee977 100644 --- a/examples/dflowfm_local/dflowfm/structures.ini +++ b/examples/dflowfm_local/dflowfm/structures.ini @@ -1,4 +1,4 @@ -# written by HYDROLIB-core 0.5.2 +# written by HYDROLIB-core 0.8.0 [General] fileVersion = 3.00 @@ -9,7 +9,7 @@ id = bridge_1 # Unique structure id (max. 256 c name = bridge_1 # Given name in the user interface. type = bridge # Structure type; must read bridge branchId = branch_12 # Branch on which the structure is located. -chainage = 14.101 # Chainage on the branch (m). +chainage = 14.101398079368284 # Chainage on the branch (m). allowedFlowdir = positive # Possible values: both, positive, negative, none. csDefId = rect_h1.550_w2.000_cno_point # Id of Cross-Section Definition. shift = 2.34 # Vertical shift of the cross section definition [m]. Defined positive upwards. @@ -24,7 +24,7 @@ id = bridge_2 # Unique structure id (max. 256 name = bridge_2 # Given name in the user interface. type = bridge # Structure type; must read bridge branchId = branch_3 # Branch on which the structure is located. -chainage = 485.794 # Chainage on the branch (m). +chainage = 485.7938981330206 # Chainage on the branch (m). allowedFlowdir = both # Possible values: both, positive, negative, none. csDefId = rect_h1.800_w5.000_cyes_point # Id of Cross-Section Definition. shift = 2.37 # Vertical shift of the cross section definition [m]. Defined positive upwards. @@ -39,7 +39,7 @@ id = culvert_1 # Unique structure id (max. 2 name = culvert_1 # Given name in the user interface. type = culvert branchId = branch_13 # Branch on which the structure is located. -chainage = 331.832 # Chainage on the branch (m). +chainage = 331.83172489372555 # Chainage on the branch (m). allowedFlowDir = positive leftLevel = 2.56 rightLevel = 2.5 diff --git a/examples/dflowfm_local/geoms/bridges.geojson b/examples/dflowfm_local/geoms/bridges.geojson index 5e08cf6e..ea9f9591 100644 --- a/examples/dflowfm_local/geoms/bridges.geojson +++ b/examples/dflowfm_local/geoms/bridges.geojson @@ -1,8 +1,63 @@ { -"type": "FeatureCollection", -"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32647" } }, -"features": [ -{ "type": "Feature", "properties": { "length": 10.0, "pillarwidth": null, "formfactor": null, "friction": 0.023, "chainage": 14.101, "branchid": "branch_12", "outletlosscoeff": 0.2, "id": "bridge_1", "allowedflowdir": "positive", "type": "bridge", "csdefid": "rect_h1.550_w2.000_cno_point", "inletlosscoeff": 0.2, "frictiontype": "Manning", "shift": 2.34 }, "geometry": { "type": "Point", "coordinates": [ 665021.849053248995915, 1525588.928306785644963 ] } }, -{ "type": "Feature", "properties": { "length": 5.0, "pillarwidth": null, "formfactor": null, "friction": 0.023, "chainage": 485.794, "branchid": "branch_3", "outletlosscoeff": 0.82, "id": "bridge_2", "allowedflowdir": "both", "type": "bridge", "csdefid": "rect_h1.800_w5.000_cyes_point", "inletlosscoeff": 0.8, "frictiontype": "Manning", "shift": 2.37 }, "geometry": { "type": "Point", "coordinates": [ 665599.689329552114941, 1525859.454655284294859 ] } } -] + "type": "FeatureCollection", + "crs": { + "type": "name", + "properties": { + "name": "urn:ogc:def:crs:EPSG::32647" + } + }, + "features": [ + { + "type": "Feature", + "properties": { + "allowedflowdir": "positive", + "csdefid": "rect_h1.550_w2.000_cno_point", + "branchid": "branch_12", + "formfactor": null, + "inletlosscoeff": 0.2, + "pillarwidth": null, + "length": 10.0, + "frictiontype": "Manning", + "friction": 0.023, + "chainage": 14.101398079368284, + "id": "bridge_1", + "outletlosscoeff": 0.2, + "shift": 2.34, + "type": "bridge" + }, + "geometry": { + "type": "Point", + "coordinates": [ + 665022.016103448579088, + 1525588.835311949253082 + ] + } + }, + { + "type": "Feature", + "properties": { + "allowedflowdir": "both", + "csdefid": "rect_h1.800_w5.000_cyes_point", + "branchid": "branch_3", + "formfactor": null, + "inletlosscoeff": 0.8, + "pillarwidth": null, + "length": 5.0, + "frictiontype": "Manning", + "friction": 0.023, + "chainage": 485.7938981330206, + "id": "bridge_2", + "outletlosscoeff": 0.82, + "shift": 2.37, + "type": "bridge" + }, + "geometry": { + "type": "Point", + "coordinates": [ + 665599.659389415639453, + 1525859.470726222032681 + ] + } + } + ] } diff --git a/examples/dflowfm_local/geoms/crosssections.geojson b/examples/dflowfm_local/geoms/crosssections.geojson index 24def10d..7fd01523 100644 --- a/examples/dflowfm_local/geoms/crosssections.geojson +++ b/examples/dflowfm_local/geoms/crosssections.geojson @@ -2,8 +2,8 @@ "type": "FeatureCollection", "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32647" } }, "features": [ -{ "type": "Feature", "properties": { "index": "1", "crsloc_id": "river_2_517.67", "crsloc_branchid": "river_2", "crsloc_chainage": 517.667, "crsloc_shift": -10.0, "crsloc_definitionid": "rect_h10.000_w50.000_c0_point", "crsdef_id": "rect_h10.000_w50.000_c0_point", "crsdef_type": "rectangle", "crsdef_branchid": "river_2", "crsdef_height": 10.0, "crsdef_width": 50.0, "crsdef_closed": "0", "crsdef_frictionid": "Manning_0.023", "frictiontype": "Manning", "frictionvalue": 0.023, "crsdef_thalweg": 0.0, "temp_id": "river_2_517.67", "crsdef_xyzcount": null, "crsdef_xcoordinates": null, "crsdef_ycoordinates": null, "crsdef_zcoordinates": null, "crsdef_frictionids": null, "crsdef_frictionpositions": null, "crsdef_diameter": null }, "geometry": { "type": "Point", "coordinates": [ 664224.816234820173122, 1526639.322860097046942 ] } }, -{ "type": "Feature", "properties": { "index": "0", "crsloc_id": "river_1_797.90", "crsloc_branchid": "river_1", "crsloc_chainage": 797.899, "crsloc_shift": 0.0, "crsloc_definitionid": "CRS_CP_test", "crsdef_id": "CRS_CP_test", "crsdef_type": "xyz", "crsdef_branchid": "river_1", "crsdef_height": null, "crsdef_width": null, "crsdef_closed": null, "crsdef_frictionid": null, "frictiontype": "Manning", "frictionvalue": 0.023, "crsdef_thalweg": 0.0, "temp_id": "river_1_797.90", "crsdef_xyzcount": 4.0, "crsdef_xcoordinates": "663467.053905815 663519.1541187783 663661.1271991032 663715.8324227147", "crsdef_ycoordinates": "1525610.846236772 1525571.7710770494 1525466.2681457987 1525421.98296478", "crsdef_zcoordinates": "1.5 -10.0 -10.0 1.5", "crsdef_frictionids": "Manning_0.023", "crsdef_frictionpositions": "0 312.3907275552155", "crsdef_diameter": null }, "geometry": { "type": "Point", "coordinates": [ 663592.067468464956619, 1525517.21589796221815 ] } }, +{ "type": "Feature", "properties": { "index": "1", "crsloc_id": "river_2_517.67", "crsloc_branchid": "river_2", "crsloc_chainage": 517.66731089307473, "crsloc_shift": -10.0, "crsloc_definitionid": "rect_h10.000_w50.000_c0_point", "crsdef_id": "rect_h10.000_w50.000_c0_point", "crsdef_type": "rectangle", "crsdef_branchid": "river_2", "crsdef_height": 10.0, "crsdef_width": 50.0, "crsdef_closed": "0", "crsdef_frictionid": "Manning_0.023", "frictiontype": "Manning", "frictionvalue": 0.023, "crsdef_thalweg": 0.0, "temp_id": "river_2_517.67", "crsdef_xyzcount": null, "crsdef_xcoordinates": null, "crsdef_ycoordinates": null, "crsdef_zcoordinates": null, "crsdef_frictionids": null, "crsdef_frictionpositions": null, "crsdef_diameter": null }, "geometry": { "type": "Point", "coordinates": [ 664218.310854589915834, 1526639.322860097046942 ] } }, +{ "type": "Feature", "properties": { "index": "0", "crsloc_id": "river_1_797.90", "crsloc_branchid": "river_1", "crsloc_chainage": 797.89912499925765, "crsloc_shift": 0.0, "crsloc_definitionid": "CRS_CP_test", "crsdef_id": "CRS_CP_test", "crsdef_type": "xyz", "crsdef_branchid": "river_1", "crsdef_height": null, "crsdef_width": null, "crsdef_closed": null, "crsdef_frictionid": null, "frictiontype": "Manning", "frictionvalue": 0.023, "crsdef_thalweg": 0.0, "temp_id": "river_1_797.90", "crsdef_xyzcount": 4.0, "crsdef_xcoordinates": "663467.053905815 663519.1541187783 663661.1271991032 663715.8324227147", "crsdef_ycoordinates": "1525610.846236772 1525571.7710770494 1525466.2681457987 1525421.98296478", "crsdef_zcoordinates": "1.5 -10.0 -10.0 1.5", "crsdef_frictionids": "Manning_0.023", "crsdef_frictionpositions": "0 312.3907275552155", "crsdef_diameter": null }, "geometry": { "type": "Point", "coordinates": [ 663683.65676968137268, 1525448.029921997338533 ] } }, { "type": "Feature", "properties": { "index": "branch_1", "crsloc_id": "branch_1_260.66", "crsloc_branchid": "branch_1", "crsloc_chainage": 260.66498236970841, "crsloc_shift": -10.0, "crsloc_definitionid": "rect_h10.000_w10.000_cno_branch", "crsdef_id": "rect_h10.000_w10.000_cno_branch", "crsdef_type": "rectangle", "crsdef_branchid": "branch_1", "crsdef_height": 10.0, "crsdef_width": 10.0, "crsdef_closed": "0", "crsdef_frictionid": "Manning_0.023", "frictiontype": "Manning", "frictionvalue": 0.023, "crsdef_thalweg": 0.0, "temp_id": "branch_1_260.66", "crsdef_xyzcount": null, "crsdef_xcoordinates": null, "crsdef_ycoordinates": null, "crsdef_zcoordinates": null, "crsdef_frictionids": null, "crsdef_frictionpositions": null, "crsdef_diameter": null }, "geometry": { "type": "Point", "coordinates": [ 665236.284100000746548, 1525175.886749999597669 ] } }, { "type": "Feature", "properties": { "index": "branch_2", "crsloc_id": "branch_2_13.86", "crsloc_branchid": "branch_2", "crsloc_chainage": 13.857973773609158, "crsloc_shift": -10.0, "crsloc_definitionid": "rect_h10.000_w10.000_cno_branch", "crsdef_id": "rect_h10.000_w10.000_cno_branch", "crsdef_type": "rectangle", "crsdef_branchid": "branch_2", "crsdef_height": 10.0, "crsdef_width": 10.0, "crsdef_closed": "0", "crsdef_frictionid": "Manning_0.023", "frictiontype": "Manning", "frictionvalue": 0.023, "crsdef_thalweg": 0.0, "temp_id": "branch_2_13.86", "crsdef_xyzcount": null, "crsdef_xcoordinates": null, "crsdef_ycoordinates": null, "crsdef_zcoordinates": null, "crsdef_frictionids": null, "crsdef_frictionpositions": null, "crsdef_diameter": null }, "geometry": { "type": "Point", "coordinates": [ 665363.387300000526011, 1525419.212349999696016 ] } }, { "type": "Feature", "properties": { "index": "branch_3", "crsloc_id": "branch_3_287.44", "crsloc_branchid": "branch_3", "crsloc_chainage": 287.43781358077678, "crsloc_shift": -10.0, "crsloc_definitionid": "rect_h10.000_w10.000_cno_branch", "crsdef_id": "rect_h10.000_w10.000_cno_branch", "crsdef_type": "rectangle", "crsdef_branchid": "branch_3", "crsdef_height": 10.0, "crsdef_width": 10.0, "crsdef_closed": "0", "crsdef_frictionid": "Manning_0.023", "frictiontype": "Manning", "frictionvalue": 0.023, "crsdef_thalweg": 0.0, "temp_id": "branch_3_287.44", "crsdef_xyzcount": null, "crsdef_xcoordinates": null, "crsdef_ycoordinates": null, "crsdef_zcoordinates": null, "crsdef_frictionids": null, "crsdef_frictionpositions": null, "crsdef_diameter": null }, "geometry": { "type": "Point", "coordinates": [ 665505.848199999891222, 1525684.700599999865517 ] } }, diff --git a/hydromt_delft3dfm/workflows/branches.py b/hydromt_delft3dfm/workflows/branches.py index 65ec6310..cf705f1d 100644 --- a/hydromt_delft3dfm/workflows/branches.py +++ b/hydromt_delft3dfm/workflows/branches.py @@ -9,10 +9,12 @@ import pyproj import shapely from hydromt import gis_utils +from hydromt.gis_utils import nearest_merge from scipy.spatial import distance -from shapely.geometry import LineString, MultiLineString, Point +from shapely.geometry import LineString, MultiLineString, MultiPoint, Point + +from hydromt_delft3dfm import graph_utils, mesh_utils -from .. import graph_utils, mesh_utils from ..gis_utils import cut_pieces, split_lines logger = logging.getLogger(__name__) @@ -1033,112 +1035,105 @@ def possibly_intersecting( return idx -# TODO copied from dhydamo geometry.py, update when available in main -# NOTE add option to write distance to nearest branch def find_nearest_branch( branches: gpd.GeoDataFrame, geometries: gpd.GeoDataFrame, method: str = "overal", - maxdist: int = 5, -): - """ - Determine nearest branch for each geometry. + maxdist: float = 5.0, +) -> gpd.GeoDataFrame: + """Determine the nearest branch for each geometry. - The nearest branch can be found by finding t from both ends (ends) or the - nearest branch from the geometry as a whole (overal), the centroid (centroid), - or intersecting (intersect). + The method of determination can vary. Parameters ---------- - branches : geopandas.GeoDataFrame - Geodataframe with branches. - geometries : geopandas.GeoDataFrame - Geodataframe with geometries to snap. - method : {'overal','intersecting','centroid','ends'} - Method for determine branch. Defaults to 'overal'. - maxdist: int or float - Maximum distance for finding nearest geometry. Defaults to 5. + branches : gpd.GeoDataFrame + Geodataframe containing branch geometries. + geometries : gpd.GeoDataFrame + Geodataframe containing geometries for which the nearest branch needs to + be found. + method : str, optional + Method to determine the nearest branch. Supports: + - 'overal': Find the nearest branch based on the geometry's location. + - 'intersecting': Convert the geometry to a centroid of its intersection + points with branches. + Default is 'overal'. + maxdist : float, optional + Maximum distance threshold for finding the nearest branch. Default is 5.0. + + Returns + ------- + gpd.GeoDataFrame + Geodataframe with additional columns: + - 'branch_id': ID of the nearest branch. + - 'branch_distance': Distance to the nearest branch. + - 'branch_offset': Offset along the branch. + + Raises + ------ + NotImplementedError + If the specified method is not among the allowed methods. """ # Check if method is in allowed methods - allowed_methods = ["intersecting", "overal", "centroid", "ends"] + allowed_methods = ["intersecting", "overal"] if method not in allowed_methods: raise NotImplementedError(f'Method "{method}" not implemented.') - # Add columns if not present - if "branch_id" not in geometries.columns: - geometries["branch_id"] = "" - if "branch_offset" not in geometries.columns: - geometries["branch_offset"] = np.nan - + # Depending on method, modify geometries if method == "intersecting": - # Determine intersection geometries per branch - geobounds = geometries.bounds.values.T - for branch in branches.itertuples(): - selectie = geometries.loc[ - possibly_intersecting(geobounds, branch.geometry) - ].copy() - intersecting = selectie.loc[selectie.intersects(branch.geometry).values] - - # For each geometrie, determine offset along branch - for geometry in intersecting.itertuples(): - # Determine distance of profile line along branch - geometries.at[geometry.Index, "branch_id"] = branch.Index - - # Calculate offset - branchgeo = branch.geometry - mindist = min(0.1, branchgeo.length / 2.0) - offset = round( - branchgeo.project( - branchgeo.intersection(geometry.geometry).centroid - ), - 3, - ) - offset = max(mindist, min(branchgeo.length - mindist, offset)) - geometries.at[geometry.Index, "branch_offset"] = offset + # Get the intersection points directly + for index, geom in geometries.iterrows(): + # Find branches that intersect with the current geometry + intersected_branches = branches[branches.intersects(geom["geometry"])] + + if not intersected_branches.empty: + # If there are multiple intersecting points, take the centroid + intersection_points = [ + geom["geometry"].intersection(branch["geometry"]) + for _, branch in intersected_branches.iterrows() + ] + centroid = MultiPoint(intersection_points).centroid + geometries.at[index, "geometry"] = centroid + + # Check for previous data and drop if exist + geometries.drop( + columns=["branch_id", "branch_distance", "branch_offset"], + inplace=True, + errors="ignore", + ) - else: - branch_bounds = branches.bounds.values.T - # In case of looking for the nearest, it is easier to iteratie over - # the geometries instead of the branches - for geometry in geometries.itertuples(): - # Find near branches - nearidx = possibly_intersecting( - branch_bounds, geometry.geometry, buffer=maxdist - ) - selectie = branches.loc[nearidx] - - if method == "overal": - # Determine distances to branches - dist = selectie.distance(geometry.geometry) - elif method == "centroid": - # Determine distances to branches - dist = selectie.distance(geometry.geometry.centroid) - elif method == "ends": - # Since a culvert can cross a channel, it is - crds = geometry.geometry.coords[:] - dist = ( - selectie.distance(Point(*crds[0])) - + selectie.distance(Point(*crds[-1])) - ) * 0.5 - - # Determine nearest - if dist.min() < maxdist: - branchidxmin = dist.idxmin() - geometries.at[geometry.Index, "branch_id"] = dist.idxmin() - geometries.at[geometry.Index, "branch_distance"] = dist.min() - if isinstance(geometry.geometry, Point): - geo = geometry.geometry - else: - geo = geometry.geometry.centroid - - # Calculate offset - branchgeo = branches.at[branchidxmin, "geometry"] - mindist = min(0.1, branchgeo.length / 2.0) - offset = max( - mindist, - min(branchgeo.length - mindist, round(branchgeo.project(geo), 3)), - ) - geometries.at[geometry.Index, "branch_offset"] = offset + # Use nearest_merge to get the nearest branches + result = nearest_merge(geometries, branches, max_dist=maxdist, columns=["geometry"]) + result.rename( + columns={"index_right": "branch_id", "distance_right": "branch_distance"}, + inplace=True, + ) + + # Select ones that are merged + valid_rows = result["branch_distance"] < maxdist + + # Interpolate the branch geometries based on index_right + branchgeo = branches.loc[result.loc[valid_rows, "branch_id"], "geometry"].values + maxdist = np.array([geo.length for geo in branchgeo]) + offset = np.array( + [ + geo.project(result.loc[idx, "geometry"]) + for idx, geo in zip(valid_rows[valid_rows].index, branchgeo) + ] + ) + result.loc[valid_rows, "branch_offset"] = np.where( + offset > maxdist, maxdist, offset + ) + snapped_geometries = [geo.interpolate(o) for geo, o in zip(branchgeo, offset)] + result.loc[valid_rows, "geometry"] = snapped_geometries + + # For rows where distance is greater than maxdist, set branch_id + # to empty and branch_offset to NaN + result.loc[~valid_rows, "branch_id"] = "" + result.loc[~valid_rows, "branch_offset"] = np.nan + result.loc[~valid_rows, "branch_distance"] = np.nan + + return result def snap_newbranches_to_branches_at_snappednodes( @@ -1207,9 +1202,8 @@ def snap_newbranches_to_branches_at_snappednodes( ].branch_chainage.to_list() snapped_line = MultiLineString(cut_pieces(branch.geometry, distances)) branches_snapped.at[branch_name, "geometry"] = snapped_line - branches_snapped.at[ - branch_name, "branchorder" - ] = branch_order # allow interpolation on the snapped branch + # allow interpolation on the snapped branch + branches_snapped.at[branch_name, "branchorder"] = branch_order # explode multilinestring after snapping branches_snapped = branches_snapped.explode(index_parts=False) @@ -1247,7 +1241,7 @@ def snap_geom_to_branches_and_drop_nonsnapped( Returns snapped geoms with branchid and chainage. Branches must have branchid. """ - find_nearest_branch( + geoms = find_nearest_branch( branches=branches, geometries=geoms, maxdist=snap_offset, diff --git a/hydromt_delft3dfm/workflows/crosssections.py b/hydromt_delft3dfm/workflows/crosssections.py index 9612303e..1647cb59 100644 --- a/hydromt_delft3dfm/workflows/crosssections.py +++ b/hydromt_delft3dfm/workflows/crosssections.py @@ -439,7 +439,7 @@ def set_xyz_crosssections( # snap to branch # setup branch_id - snap bridges to branch (inplace of bridges, # will add branch_id and branch_offset columns) - find_nearest_branch( + crosssections = find_nearest_branch( branches=branches, geometries=crosssections, method="intersecting" ) logger.warning( @@ -502,7 +502,7 @@ def set_xyz_crosssections( crosssections.branch_id.to_list(), "frictionvalue" ], "crsdef_frictionpositions": [ - f"0 {l}" for l in crosssections.geometry.length.to_list() + f"0 {l[-1]}" for l in crosssections.l.to_list() ], } ) @@ -582,7 +582,7 @@ def set_point_crosssections( # snap to branch # setup branch_id - snap bridges to branch # (inplace of bridges, will add branch_id and branch_offset columns) - find_nearest_branch( + crosssections = find_nearest_branch( branches=branches, geometries=crosssections, method="overal", maxdist=maxdist ) diff --git a/hydromt_delft3dfm/workflows/structures.py b/hydromt_delft3dfm/workflows/structures.py index 3d5c1b58..5d4ee240 100644 --- a/hydromt_delft3dfm/workflows/structures.py +++ b/hydromt_delft3dfm/workflows/structures.py @@ -105,7 +105,9 @@ def prepare_1dstructures( # 4. snap structures to branches # setup branch_id - snap structures to branch (inplace of structures, # will add branch_id and branch_offset columns) - find_nearest_branch(branches=branches, geometries=gdf_st, maxdist=snap_offset) + gdf_st = find_nearest_branch( + branches=branches, geometries=gdf_st, maxdist=snap_offset + ) # setup failed - drop based on branch_offset that are not snapped to branch _old_ids = gdf_st.index.to_list() gdf_st.dropna(axis=0, inplace=True, subset=["branch_offset"])