Skip to content

Commit b495fa6

Browse files
pjMitchell490pjMitchell490
authored
Merge pull request #5 from pjMitchell490/main
Improved Performance, simplified code using networkx library
2 parents 222c2a2 + e7dd8d1 commit b495fa6

File tree

3 files changed

+152
-81
lines changed

3 files changed

+152
-81
lines changed

‎README.md‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,9 @@
11
# valley-centerline
22
Tool to pick the centerline of a valley (or other feature) based on the valley walls (or other edges)
33

4+
# valley-centerline-interpolation
5+
Much faster than valley-centerline, but only suitable for valleys with little meandering and mostly parallel walls
6+
47
## Ways to pick valley margins
58

69
In a GIS, create polylines on either side of the valley. The valley width that is chosen depends on the problem to be solved. For mapping terraces within the valley, it would be best to map at the break between the bedrock valley walls and the flat valley bottom and terraces, such that the defined valley encompasses all of the features of concern. Dealing with tributaries might be more difficult, but the idea that I would suggest is just to start with the simplest: snap a straight line across each tributary-valley mouth between the closest points before the valley starts to curve in towards the tributary. You might also start with a shorter portion of the river to build out a smaller valley-wall data set and see how it works with the algorithm that you come up with (which will be stored here), and then modify and/or extend it depending on how the results look.

‎valley-centerline-interpolation.py‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
import argparse
2+
import shapely
3+
import geopandas as gpd
4+
5+
SUBDIVISION_AMOUNT = 200
6+
7+
parser = argparse.ArgumentParser(description='Find the centerline between two lines.')
8+
parser.add_argument("input_wall_1", help="the first valley wall shapefile")
9+
parser.add_argument("input_wall_2", help="the second valley wall shapefile")
10+
parser.add_argument("output", help="the filepath for the output shapefile")
11+
12+
13+
args = parser.parse_args()
14+
input_1 = gpd.read_file(args.input_wall_1)
15+
input_2 = gpd.read_file(args.input_wall_2)
16+
output = args.output
17+
18+
19+
# Import line shapefiles
20+
'''
21+
22+
input_1 = gpd.read_file("C:\LSDTopoTools\Github\Valley-Centerline\WW\Input\WW_ValleyWall_N_10mpoints.shp")
23+
input_2 = gpd.read_file("C:\LSDTopoTools\Github\Valley-Centerline\WW\Input\WW_ValleyWall_S_10mpoints.shp")
24+
'''
25+
26+
27+
wall_1 = shapely.geometry.shape(input_1['geometry'][0])
28+
wall_2 = shapely.geometry.shape(input_2['geometry'][0])
29+
30+
wall_1_points = []
31+
wall_2_points = []
32+
centerline_points = []
33+
34+
for i in range(SUBDIVISION_AMOUNT + 1):
35+
wall_1_points.append(wall_1.interpolate(i/SUBDIVISION_AMOUNT, normalized=True))
36+
wall_2_points.append(wall_2.interpolate(i/SUBDIVISION_AMOUNT, normalized=True))
37+
connecting_line = shapely.geometry.LineString([wall_1_points[i], wall_2_points[i]])
38+
centerline_points.append(connecting_line.interpolate(0.5, normalized=True))
39+
centerline = shapely.geometry.LineString(centerline_points)
40+
print(centerline)
41+
42+
features = [0]
43+
geometry = [centerline]
44+
df = {'features': features, 'geometry': geometry}
45+
gdf = gpd.GeoDataFrame(df)
46+
gdf = gdf.set_crs("EPSG:32615")
47+
gdf.to_file('C:\LSDTopoTools\Github\Valley-Centerline\WW\Input\Centerline_interpolation.shp')
48+
print("Exported")

‎valley-centerline.py‎

Lines changed: 101 additions & 81 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,30 @@
1+
import argparse
12
import shapely
23
import geopandas as gpd
34
from geovoronoi import voronoi_regions_from_coords
5+
import networkx as nx
6+
import momepy
7+
import numpy as np
48

59
BUFFER_DISTANCE = 0.0001
610
SUBDIVISION_AMOUNT = 200
11+
'''
12+
parser = argparse.ArgumentParser(description='Find the centerline between two lines.')
13+
parser.add_argument("input_wall_1", help="the first valley wall shapefile")
14+
parser.add_argument("input_wall_2", help="the second valley wall shapefile")
15+
parser.add_argument("output", help="the filepath for the output shapefile")
16+
17+
18+
args = parser.parse_args()
19+
input_1 = gpd.read_file(args.input_wall_1)
20+
input_2 = gpd.read_file(args.input_wall_2)
21+
output = args.output
22+
'''
723

824
# Import line shapefiles
9-
input_1 = gpd.read_file("C:/Users/pmitc/Documents/QGIS/Zumbro/Zumbro_SamplePoints/TestGround/Testwall_1.shp")
10-
input_2 = gpd.read_file("C:/Users/pmitc/Documents/QGIS/Zumbro/Zumbro_SamplePoints/TestGround/Testwall_2.shp")
25+
26+
input_1 = gpd.read_file("C:\LSDTopoTools\Github\Valley-Centerline\WW\Input\WW_ValleyWall_N_10mpoints.shp")
27+
input_2 = gpd.read_file("C:\LSDTopoTools\Github\Valley-Centerline\WW\Input\WW_ValleyWall_S_10mpoints.shp")
1128

1229
wall_1 = shapely.geometry.shape(input_1['geometry'][0])
1330
wall_2 = shapely.geometry.shape(input_2['geometry'][0])
@@ -20,7 +37,7 @@ def subdivide_wall(wall, subdivision_distance):
2037
dist += subdivision_distance
2138
return subdivided_coords
2239

23-
extra_coords = subdivide_wall(wall_1, SUBDIVISION_AMOUNT) + subdivide_wall(wall_2, SUBDIVISION_AMOUNT)
40+
#extra_coords = subdivide_wall(wall_1, SUBDIVISION_AMOUNT) + subdivide_wall(wall_2, SUBDIVISION_AMOUNT)
2441

2542
# Convert valley walls into a collection of points for voronoi algorithm
2643

@@ -36,8 +53,11 @@ def subdivide_wall(wall, subdivision_distance):
3653
start_boundary = shapely.geometry.LineString([wall_1.coords[0], wall_2.coords[-1]])
3754
end_boundary = shapely.geometry.LineString([wall_1.coords[-1], wall_2.coords[0]])
3855

39-
coords = coords + extra_coords
40-
points = shapely.geometry.MultiPoint(coords)
56+
start_pole = start_boundary.interpolate(0.5, normalized=True)
57+
end_pole = end_boundary.interpolate(0.5, normalized=True)
58+
59+
#coords = coords + extra_coords
60+
points = [shapely.Point(i[0], i[1]) for i in coords]
4161

4262
# Generate polygon and buffer from the input walls
4363
buffer = valleypoly.buffer(BUFFER_DISTANCE)
@@ -46,101 +66,101 @@ def subdivide_wall(wall, subdivision_distance):
4666
geometry = [valleypoly]
4767
df = {'features': features, 'geometry': geometry}
4868
gdf = gpd.GeoDataFrame(df)
49-
gdf.to_file('C:/Users/pmitc/Documents/QGIS/Zumbro/Zumbro_SamplePoints/TestGround/ValleyPoly.shp')
69+
gdf = gdf.set_crs("EPSG:32615")
70+
gdf.to_file('C:\LSDTopoTools\Github\Valley-Centerline\WW\Input\ValleyPoly.shp')
5071

5172
features = [0]
5273
geometry = [buffer]
5374
df = {'features': features, 'geometry': geometry}
5475
gdf = gpd.GeoDataFrame(df)
55-
gdf.to_file('C:/Users/pmitc/Documents/QGIS/Zumbro/Zumbro_SamplePoints/TestGround/ValleyBuffer.shp')
56-
76+
gdf = gdf.set_crs("EPSG:32615")
77+
gdf.to_file('C:\LSDTopoTools\Github\Valley-Centerline\WW\Input\ValleyBuffer.shp')
78+
# print(coords)
5779
# Create Voronoi Polygons
5880
region_polys, region_pts = voronoi_regions_from_coords(points, buffer)
59-
81+
# print(type(region_polys))
82+
problem_polys = []
83+
for key, value in region_polys.items():
84+
if value.geom_type == 'MultiPolygon':
85+
problem_polys.append(key)
6086
#Export the voronoi polygons
87+
for key in problem_polys:
88+
region_polys.pop(key)
6189

6290
features = [i for i in range(len(region_polys))]
6391
geometry = [geom for geom in region_polys.values()]
6492
df = {'features': features, 'geometry': geometry}
6593
gdf = gpd.GeoDataFrame(df)
66-
gdf.to_file('C:/Users/pmitc/Documents/QGIS/Zumbro/Zumbro_SamplePoints/TestGround/Voronoi_Demo.shp')
94+
gdf = gdf.set_crs("EPSG:32615")
95+
96+
gdf.to_file('C:\LSDTopoTools\Github\Valley-Centerline\WW\Input\Voronoi.shp')
97+
6798

6899
print("Voronoi Polygons created.")
69-
voronoi_edges = []
70100

101+
102+
########################################
103+
# Find potential start/end edges #
104+
########################################
105+
edges = np.array([])
71106
for poly in region_polys.values():
72-
if poly.geom_type == 'MultiPolygon':
73-
# If the buffer distance is too small (<0.0006 in this case), some of the voronoi outputs are MultiPolygons.
74-
# I'm just discarding them here, which seems to work alright, but I'm wondering if there are cases when this would give a bad result.
75-
pass
76-
else:
77-
for i in range(len(poly.exterior.coords) - 1):
78-
edge = shapely.geometry.LineString((poly.exterior.coords[i], poly.exterior.coords[i+1]))
79-
repeat = False
80-
81-
# Edges will have duplicates if they're the border between polygons. If this is the case, we don't want to add the duplicates.
82-
for test in voronoi_edges:
83-
if edge.coords[:] == test.coords[:] or edge.coords[:] == test.coords[::-1]:
84-
# In all the cases so far, it seems that the coordinates have been flipped, so I'm not sure if this will ever be the case.
85-
repeat = True
86-
# Get rid of edges that are duplicates or outside of the valley
87-
if edge.within(valleypoly) and not repeat:
88-
voronoi_edges.append(edge)
89-
# Get the first and last segments, which should intersect the lines connecting the two valley walls
90-
if not(edge.intersects(wall_1) or edge.intersects(wall_2)):
91-
if edge.intersects(start_boundary):
92-
start_edge = edge
93-
print('Found start edge!')
94-
if edge.intersects(end_boundary):
95-
end_edge = edge
96-
print('Found end edge!')
97-
98-
# Export voronoi edges (for testing purposes)
99-
features = [i for i in range(len(voronoi_edges))]
100-
geometry = [geom for geom in voronoi_edges]
107+
for i in range(len(poly.exterior.coords) - 1):
108+
edge = shapely.geometry.LineString((poly.exterior.coords[i], poly.exterior.coords[i+1]))
109+
if not(edge.intersects(wall_1) or edge.intersects(wall_2)):
110+
edges = np.append(edges, edge)
111+
112+
features = [i for i in range(len(edges))]
113+
geometry = edges
101114
df = {'features': features, 'geometry': geometry}
102115
gdf = gpd.GeoDataFrame(df)
103-
gdf.to_file('C:/Users/pmitc/Documents/QGIS/Zumbro/Zumbro_SamplePoints/TestGround/Voronoi_Edges_Demo.shp')
104-
105-
tried_edges = []
106-
107-
# Recursive solution for getting only the centerline from beginning to end
108-
def get_centerline_path(edge, from_edge):
109-
viable_edges = []
110-
global centerlineedges
111-
global tried_edges
112-
# Add self to the list of tried edges
113-
tried_edges.append(edge)
114-
if edge.touches(end_edge):
115-
# If we've found the end, add ourselves to the list and let the edge that found us know
116-
centerlineedges.append(edge)
117-
centerlineedges.append(end_edge)
118-
print("We've reached the end!")
119-
return True
120-
# See which neighbors haven't been checked yet
121-
for other_edge in voronoi_edges:
122-
if edge.touches(other_edge) and other_edge not in tried_edges:
123-
# Don't add the other edge if it's one connected to the parent edge
124-
if from_edge is None or not(from_edge.touches(other_edge)):
125-
viable_edges.append(other_edge)
126-
for other_edge in viable_edges:
127-
if get_centerline_path(other_edge, edge):
128-
# If one of our viable edges is on the path toward the end, add ourselves to the list and let the edge that found us know
129-
centerlineedges.append(edge)
130-
return True
131-
return False
132-
133-
centerlineedges = []
134-
get_centerline_path(start_edge, None)
135-
136-
# Export centerline
137-
centerline = shapely.geometry.MultiLineString(centerlineedges)
138-
centerline_merged = shapely.ops.linemerge(centerline)
139-
print(centerline)
140-
features = [0]
141-
geometry = [centerline_merged]
116+
edges_gdf = gdf.set_crs("EPSG:32615")
117+
edges_gdf.to_file('C:\LSDTopoTools\Github\Valley-Centerline\WW\Input\edges.shp')
118+
print("Exported edges")
119+
120+
features = [0, 1]
121+
geometry = [start_pole, end_pole]
142122
df = {'features': features, 'geometry': geometry}
143123
gdf = gpd.GeoDataFrame(df)
144-
gdf = gdf.set_crs("EPSG:32615")
145-
gdf.to_file('C:/Users/pmitc/Documents/QGIS/Zumbro/Zumbro_SamplePoints/TestGround/Centerline.shp')
124+
poles_gdf = gdf.set_crs("EPSG:32615")
125+
poles_gdf.to_file('C:\LSDTopoTools\Github\Valley-Centerline\WW\Input\poles.shp')
146126
print("Exported")
127+
128+
129+
130+
########################################
131+
# Find nearest edges to poles #
132+
########################################
133+
nearest_edges = gpd.sjoin_nearest(poles_gdf, edges_gdf).merge(edges_gdf, left_on="index_right", right_index=True)
134+
cols = list(nearest_edges)
135+
print(cols)
136+
print("Found nearest edges!")
137+
features = [i for i in range(len(nearest_edges['geometry_y']))]
138+
geometry = nearest_edges['geometry_y']
139+
df = {'features': features, 'geometry': geometry}
140+
nearest_edges = gpd.GeoDataFrame(df)
141+
nearest_edges.to_file('C:\LSDTopoTools\Github\Valley-Centerline\WW\Input\start_end.shp')
142+
143+
########################################
144+
# Find shortest path from start to end #
145+
########################################
146+
147+
graph = momepy.gdf_to_nx(edges_gdf)
148+
149+
start = nearest_edges['geometry'][0].coords[0]
150+
end = nearest_edges['geometry'][1].coords[0]
151+
path = nx.shortest_path(graph, source=start, target=end)
152+
print('Path found')
153+
print(path[0:5])
154+
path = shapely.LineString(path)
155+
156+
########################################
157+
# Export centerline #
158+
########################################
159+
160+
features = [0]
161+
geometry = path
162+
df = {'features': features, 'geometry': geometry}
163+
centerline_gdf = gpd.GeoDataFrame(df)
164+
centerline_gdf = centerline_gdf.set_crs("EPSG:32615")
165+
centerline_gdf.to_file('C:\LSDTopoTools\Github\Valley-Centerline\WW\Input\centerline.shp')
166+
print("Exported")

0 commit comments

Comments
 (0)