在与其他线串的交叉点处整齐地分割线串
Shapely Split LineStrings at Intersections with other LineStrings
我有一组 LineString，它们与其他 LineString 相交，我想在这些交点处将 LineString 拆分为单独的线段。我有一个解决方案，但我认为这不是最好的方法。
假设我们正在处理一个 LineString：
>>> import shapely
>>> from shapely.geometry import *
>>> import geopandas as gpd
>>>
>>> MyLine=LineString([(0,0),(5,0),(10,3)])
>>> MyLine
<shapely.geometry.linestring.LineString object at 0x1277EEB0>
>>>
与此 LineString 相交的 2 条线：
>>> IntersectionLines=gpd.GeoSeries([LineString([(2,1.5),(3,-.5)]), LineString([(5,.5),(7,.5)])])
>>> IntersectionLines
0 LINESTRING (2 1.5, 3 -0.5)
1 LINESTRING (5 0.5, 7 0.5)
dtype: object
>>>
看起来像这样：
我可以得到交点如下：
>>> IntPoints=MyLine.intersection(IntersectionLines.unary_union)
>>> IntPointCoords=[x.coords[:][0] for x in IntPoints]
>>> IntPointCoords
[(2.75, 0.0), (5.833333333333333, 0.5)]
>>>
然后我提取起点和终点，并用这些点和将用于形成线段的交点创建对：
>>> StartPoint=MyLine.coords[0]
>>> EndPoint=MyLine.coords[-1]
>>> SplitCoords=[StartPoint]+IntPointCoords+[EndPoint]
>>>
>>> def pair(list):
... for i in range(1, len(list)):
... yield list[i-1], list[i]
...
>>>
>>> SplitSegments=[LineString(p) for p in pair(SplitCoords)]
>>> SplitSegments
[<shapely.geometry.linestring.LineString object at 0x127F7A90>, <shapely.geometry.linestring.LineString object at 0x127F7570>, <shapely.geometry.linestring.LineString object at 0x127F7930>]
>>>
>>> gpd.GeoSeries(SplitSegments)
0 LINESTRING (0 0, 2.75 0)
1 LINESTRING (2.75 0, 5.833333333333333 0.5)
2 LINESTRING (5.833333333333333 0.5, 10 3)
dtype: object
>>>
但是，我的方法有一个问题是我知道哪个交点应该与 LineString 起点连接，哪个交点应该与 LineString 终点配对。如果交点以与起点和终点相同的顺序沿线列出，则该程序有效。我想在某些情况下情况并非总是如此？有没有一种方法可以推广这种方法，或者是否有更好的方法？
这里有一个更通用的方法：计算所有点（线的起点和终点+要分割的点）沿线的距离，按这些点排序，然后生成线段在正确的顺序。一起参加一个活动：
def cut_line_at_points(line, points):
# First coords of line (start + end)
coords = [line.coords[0], line.coords[-1]]
# Add the coords from the points
coords += [list(p.coords)[0] for p in points]
# Calculate the distance along the line for each point
dists = [line.project(Point(p)) for p in coords]
# sort the coords based on the distances
# see
coords = [p for (d, p) in sorted(zip(dists, coords))]
# generate the Lines
lines = [LineString([coords[i], coords[i+1]]) for i in range(len(coords)-1)]
return lines
在您的示例中应用此函数：
In [13]: SplitSegments = cut_line_at_points(MyLine, IntPoints)
In [14]: gpd.GeoSeries(SplitSegments)
Out[14]:
0 LINESTRING (0 0, 2.75 0)
1 LINESTRING (2.75 0, 5.833333333333333 0.5)
2 LINESTRING (5.833333333333333 0.5, 10 3)
dtype: object
唯一的问题是这不会保留原始线的角（但问题中的示例也没有这样做，所以我不知道这是否是一个要求。这是可能的，但要做到有点复杂）
Update 保持原线角完整的版本（我的做法是还保留一个0/1的列表，表示一个坐标是要分割还是不是):
def cut_line_at_points(line, points):
# First coords of line
coords = list(line.coords)
# Keep list coords where to cut (cuts = 1)
cuts = [0] * len(coords)
cuts[0] = 1
cuts[-1] = 1
# Add the coords from the points
coords += [list(p.coords)[0] for p in points]
cuts += [1] * len(points)
# Calculate the distance along the line for each point
dists = [line.project(Point(p)) for p in coords]
# sort the coords/cuts based on the distances
# see
coords = [p for (d, p) in sorted(zip(dists, coords))]
cuts = [p for (d, p) in sorted(zip(dists, cuts))]
# generate the Lines
#lines = [LineString([coords[i], coords[i+1]]) for i in range(len(coords)-1)]
lines = []
for i in range(len(coords)-1):
if cuts[i] == 1:
# find next element in cuts == 1 starting from index i + 1
j = cuts.index(1, i + 1)
lines.append(LineString(coords[i:j+1]))
return lines
应用于示例：
In [3]: SplitSegments = cut_line_at_points(MyLine, IntPoints)
In [4]: gpd.GeoSeries(SplitSegments)
Out[4]:
0 LINESTRING (0 0, 2.75 0)
1 LINESTRING (2.75 0, 5 0, 5.833333333333333 0.5)
2 LINESTRING (5.833333333333333 0.5, 10 3)
dtype: object
这是我尝试调整 joris 的函数，以便也包括线段的角。这还不能完美地工作，因为除了包括包含角的合并段之外，它还包括原始的未合并段。
def cut_line_at_points(line, points):
#make the coordinate list all of the coords that define the line
coords=line.coords[:]
coords += [list(p.coords)[0] for p in points]
dists = [line.project(Point(p)) for p in coords]
coords = [p for (d, p) in sorted(zip(dists, coords))]
lines = [LineString([coords[i], coords[i+1]]) for i in range(len(coords)-1)]
#Now go through the lines and merge together as one segment if there is no point interrupting it
CleanedLines=[]
for i,line in enumerate(lines):
if i<>len(lines)-1:
LinePair=[line,lines[i+1]]
IntPoint= LinePair[0].intersection(LinePair[1])
if IntPoint not in points:
CleanedLine=shapely.ops.linemerge(LinePair)
else:
CleanedLine=line
else:
CleanedLine=line
CleanedLines.append(CleanedLine)
return CleanedLines
>>> SplitSegments = cut_line_at_points(MyLine, IntPoints)
>>> gpd.GeoSeries(SplitSegments)
0 LINESTRING (0 0, 2.75 0)
1 LINESTRING (2.75 0, 5 0, 5.833333333333333 0.5)
2 LINESTRING (5 0, 5.833333333333333 0.5)
3 LINESTRING (5.833333333333333 0.5, 10 3)
dtype: object
>>>
我喜欢 joris 的方法。不幸的是，我 运行 在尝试使用它时遇到了一个严重的困难：如果线串有两个点在相同的坐标，它们的投影是不明确的。两者将获得相同的投影值并一起排序。
如果您的路径在同一点开始和结束，这一点尤其明显。结束点的投影为 0 并在开始时进行排序，这会抛出整个算法，因为它期望 "cuts" 末尾的值为“1”。
这是一个适用于 shapely 1.6.1 的解决方案：
import shapely.ops
from shapely.geometry import MultiPoint
def cut_linestring_at_points(linestring, points):
return shapely.ops.split(linestring, MultiPoint(points))
是的，就是这么简单。这里要注意的是点 必须 正好在线上。如果不是，将它们对齐到线，如 this answer.
return 值是一个 MultiLineString，您可以使用其 geoms 方法获取组件 LineStrings。
@joris 的方法非常好，但是如果您尝试将点列表传递给它，它就会出错，其中一些点实际上并不与直线相交，在我的情况下发生这种情况是因为我预先计算了一个交点列表来自许多行的列表。
我能够通过将输入点列表预过滤到仅那些在继续执行该函数之前实际相交的点来修复它。对于大量的点列表来说效率不高，但就我而言，我的列表总是很小，所以对我来说已经足够了。如果没有点与线相交，它也可以工作，并且在这种情况下只会短路返回原始线作为列表（为了保持一致性）
我最初使用 line.intersects(point) 但它总是返回 False，可能是由于插值精度。
def cut_line_at_points(line, points):
# Filter out any points that are not on the line
# 0.01 is arbitrary, make it smaller for more precision
points = [point for point in points if line.distance(point) < 0.01]
if not points:
return [line]
# First coords of line
coords = list(line.coords)
# Keep list coords where to cut (cuts = 1)
cuts = [0] * len(coords)
cuts[0] = 1
cuts[-1] = 1
# Add the coords from the points
coords += [list(p.coords)[0] for p in points]
cuts += [1] * len(points)
# Calculate the distance along the line for each point
dists = [line.project(Point(p)) for p in coords]
# sort the coords/cuts based on the distances
# see
coords = [p for (d, p) in sorted(zip(dists, coords))]
cuts = [p for (d, p) in sorted(zip(dists, cuts))]
# generate the Lines
# lines = [LineString([coords[i], coords[i+1]]) for i in range(len(coords)-1)]
lines = []
for i in range(len(coords) - 1):
if cuts[i] == 1:
# find next element in cuts == 1 starting from index i + 1
j = cuts.index(1, i + 1)
lines.append(LineString(coords[i:j + 1]))
return lines
我有一组 LineString，它们与其他 LineString 相交，我想在这些交点处将 LineString 拆分为单独的线段。我有一个解决方案，但我认为这不是最好的方法。
假设我们正在处理一个 LineString：
>>> import shapely
>>> from shapely.geometry import *
>>> import geopandas as gpd
>>>
>>> MyLine=LineString([(0,0),(5,0),(10,3)])
>>> MyLine
<shapely.geometry.linestring.LineString object at 0x1277EEB0>
>>>
与此 LineString 相交的 2 条线：
>>> IntersectionLines=gpd.GeoSeries([LineString([(2,1.5),(3,-.5)]), LineString([(5,.5),(7,.5)])])
>>> IntersectionLines
0 LINESTRING (2 1.5, 3 -0.5)
1 LINESTRING (5 0.5, 7 0.5)
dtype: object
>>>
看起来像这样：
我可以得到交点如下：
>>> IntPoints=MyLine.intersection(IntersectionLines.unary_union)
>>> IntPointCoords=[x.coords[:][0] for x in IntPoints]
>>> IntPointCoords
[(2.75, 0.0), (5.833333333333333, 0.5)]
>>>
然后我提取起点和终点，并用这些点和将用于形成线段的交点创建对：
>>> StartPoint=MyLine.coords[0]
>>> EndPoint=MyLine.coords[-1]
>>> SplitCoords=[StartPoint]+IntPointCoords+[EndPoint]
>>>
>>> def pair(list):
... for i in range(1, len(list)):
... yield list[i-1], list[i]
...
>>>
>>> SplitSegments=[LineString(p) for p in pair(SplitCoords)]
>>> SplitSegments
[<shapely.geometry.linestring.LineString object at 0x127F7A90>, <shapely.geometry.linestring.LineString object at 0x127F7570>, <shapely.geometry.linestring.LineString object at 0x127F7930>]
>>>
>>> gpd.GeoSeries(SplitSegments)
0 LINESTRING (0 0, 2.75 0)
1 LINESTRING (2.75 0, 5.833333333333333 0.5)
2 LINESTRING (5.833333333333333 0.5, 10 3)
dtype: object
>>>
但是，我的方法有一个问题是我知道哪个交点应该与 LineString 起点连接，哪个交点应该与 LineString 终点配对。如果交点以与起点和终点相同的顺序沿线列出，则该程序有效。我想在某些情况下情况并非总是如此？有没有一种方法可以推广这种方法，或者是否有更好的方法？
这里有一个更通用的方法：计算所有点（线的起点和终点+要分割的点）沿线的距离，按这些点排序，然后生成线段在正确的顺序。一起参加一个活动：
def cut_line_at_points(line, points):
# First coords of line (start + end)
coords = [line.coords[0], line.coords[-1]]
# Add the coords from the points
coords += [list(p.coords)[0] for p in points]
# Calculate the distance along the line for each point
dists = [line.project(Point(p)) for p in coords]
# sort the coords based on the distances
# see
coords = [p for (d, p) in sorted(zip(dists, coords))]
# generate the Lines
lines = [LineString([coords[i], coords[i+1]]) for i in range(len(coords)-1)]
return lines
在您的示例中应用此函数：
In [13]: SplitSegments = cut_line_at_points(MyLine, IntPoints)
In [14]: gpd.GeoSeries(SplitSegments)
Out[14]:
0 LINESTRING (0 0, 2.75 0)
1 LINESTRING (2.75 0, 5.833333333333333 0.5)
2 LINESTRING (5.833333333333333 0.5, 10 3)
dtype: object
唯一的问题是这不会保留原始线的角（但问题中的示例也没有这样做，所以我不知道这是否是一个要求。这是可能的，但要做到有点复杂）
Update 保持原线角完整的版本（我的做法是还保留一个0/1的列表，表示一个坐标是要分割还是不是):
def cut_line_at_points(line, points):
# First coords of line
coords = list(line.coords)
# Keep list coords where to cut (cuts = 1)
cuts = [0] * len(coords)
cuts[0] = 1
cuts[-1] = 1
# Add the coords from the points
coords += [list(p.coords)[0] for p in points]
cuts += [1] * len(points)
# Calculate the distance along the line for each point
dists = [line.project(Point(p)) for p in coords]
# sort the coords/cuts based on the distances
# see
coords = [p for (d, p) in sorted(zip(dists, coords))]
cuts = [p for (d, p) in sorted(zip(dists, cuts))]
# generate the Lines
#lines = [LineString([coords[i], coords[i+1]]) for i in range(len(coords)-1)]
lines = []
for i in range(len(coords)-1):
if cuts[i] == 1:
# find next element in cuts == 1 starting from index i + 1
j = cuts.index(1, i + 1)
lines.append(LineString(coords[i:j+1]))
return lines
应用于示例：
In [3]: SplitSegments = cut_line_at_points(MyLine, IntPoints)
In [4]: gpd.GeoSeries(SplitSegments)
Out[4]:
0 LINESTRING (0 0, 2.75 0)
1 LINESTRING (2.75 0, 5 0, 5.833333333333333 0.5)
2 LINESTRING (5.833333333333333 0.5, 10 3)
dtype: object
这是我尝试调整 joris 的函数，以便也包括线段的角。这还不能完美地工作，因为除了包括包含角的合并段之外，它还包括原始的未合并段。
def cut_line_at_points(line, points):
#make the coordinate list all of the coords that define the line
coords=line.coords[:]
coords += [list(p.coords)[0] for p in points]
dists = [line.project(Point(p)) for p in coords]
coords = [p for (d, p) in sorted(zip(dists, coords))]
lines = [LineString([coords[i], coords[i+1]]) for i in range(len(coords)-1)]
#Now go through the lines and merge together as one segment if there is no point interrupting it
CleanedLines=[]
for i,line in enumerate(lines):
if i<>len(lines)-1:
LinePair=[line,lines[i+1]]
IntPoint= LinePair[0].intersection(LinePair[1])
if IntPoint not in points:
CleanedLine=shapely.ops.linemerge(LinePair)
else:
CleanedLine=line
else:
CleanedLine=line
CleanedLines.append(CleanedLine)
return CleanedLines
>>> SplitSegments = cut_line_at_points(MyLine, IntPoints)
>>> gpd.GeoSeries(SplitSegments)
0 LINESTRING (0 0, 2.75 0)
1 LINESTRING (2.75 0, 5 0, 5.833333333333333 0.5)
2 LINESTRING (5 0, 5.833333333333333 0.5)
3 LINESTRING (5.833333333333333 0.5, 10 3)
dtype: object
>>>
我喜欢 joris 的方法。不幸的是，我 运行 在尝试使用它时遇到了一个严重的困难：如果线串有两个点在相同的坐标，它们的投影是不明确的。两者将获得相同的投影值并一起排序。
如果您的路径在同一点开始和结束，这一点尤其明显。结束点的投影为 0 并在开始时进行排序，这会抛出整个算法，因为它期望 "cuts" 末尾的值为“1”。
这是一个适用于 shapely 1.6.1 的解决方案：
import shapely.ops
from shapely.geometry import MultiPoint
def cut_linestring_at_points(linestring, points):
return shapely.ops.split(linestring, MultiPoint(points))
是的，就是这么简单。这里要注意的是点 必须 正好在线上。如果不是，将它们对齐到线，如 this answer.
return 值是一个 MultiLineString，您可以使用其 geoms 方法获取组件 LineStrings。
@joris 的方法非常好，但是如果您尝试将点列表传递给它，它就会出错，其中一些点实际上并不与直线相交，在我的情况下发生这种情况是因为我预先计算了一个交点列表来自许多行的列表。
我能够通过将输入点列表预过滤到仅那些在继续执行该函数之前实际相交的点来修复它。对于大量的点列表来说效率不高，但就我而言，我的列表总是很小，所以对我来说已经足够了。如果没有点与线相交，它也可以工作，并且在这种情况下只会短路返回原始线作为列表（为了保持一致性）
我最初使用 line.intersects(point) 但它总是返回 False，可能是由于插值精度。
def cut_line_at_points(line, points):
# Filter out any points that are not on the line
# 0.01 is arbitrary, make it smaller for more precision
points = [point for point in points if line.distance(point) < 0.01]
if not points:
return [line]
# First coords of line
coords = list(line.coords)
# Keep list coords where to cut (cuts = 1)
cuts = [0] * len(coords)
cuts[0] = 1
cuts[-1] = 1
# Add the coords from the points
coords += [list(p.coords)[0] for p in points]
cuts += [1] * len(points)
# Calculate the distance along the line for each point
dists = [line.project(Point(p)) for p in coords]
# sort the coords/cuts based on the distances
# see
coords = [p for (d, p) in sorted(zip(dists, coords))]
cuts = [p for (d, p) in sorted(zip(dists, cuts))]
# generate the Lines
# lines = [LineString([coords[i], coords[i+1]]) for i in range(len(coords)-1)]
lines = []
for i in range(len(coords) - 1):
if cuts[i] == 1:
# find next element in cuts == 1 starting from index i + 1
j = cuts.index(1, i + 1)
lines.append(LineString(coords[i:j + 1]))
return lines