python - 在与其他 LineStrings 的交点处匀称分割 LineStrings

Question

我有一组 LineStrings 与其他 LineStrings 相交，我想在这些交点处将 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 终点配对。如果交叉点以与起点和终点相同的顺序沿线列出，则此程序有效。我想在某些情况下情况并非总是如此？有没有办法推广这种方法，或者有更好的方法吗？

score 6 · Accepted Answer

这是一种更通用的方法：计算所有点（线的起点和终点+要分割的点）沿线的距离，按这些点排序，然后以正确的顺序生成线段。一起在一个函数中：

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 http://stackoverflow.com/questions/6618515/sorting-list-based-on-values-from-another-list
    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

唯一的问题是，这不会保留原始行的角（但是您在问题中的示例也没有这样做，所以我不知道这是否是一个要求。这是可能的，但会使它更复杂一些）

更新保持原始行中的角完整的版本（我的方法是还保留一个 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 http://stackoverflow.com/questions/6618515/sorting-list-based-on-values-from-another-list    
    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

score 2 · Accepted Answer

我喜欢joris的方法。不幸的是，我在尝试使用它时遇到了一个严重的困难：如果线串在同一坐标上有两个点，那么它们的投影是不明确的。两者都将获得相同的投影值并一起排序。

如果您有一条在同一点开始和结束的路径，这一点尤其明显。结束点的投影为 0 并在开始时排序，这会抛出整个算法，因为它期望最后的“cuts”值为“1”。

这是一个适用于 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))

是的，真的就是这么简单。这里的问题是这些点必须完全在线。如果不是，请按此答案将它们对齐。

返回值为 a MultiLineString，您可以LineString使用其geoms方法获取组件 s。

score 1 · Accepted Answer

这是我尝试通过 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
>>>

score 0 · Accepted Answer

@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 http://stackoverflow.com/questions/6618515/sorting-list-based-on-values-from-another-list
    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

python - 在与其他 LineStrings 的交点处匀称分割 LineStrings

4 回答 4

Related

Reference