在与其他LineString相交处形状拆分LineStrings

这里有一个更通用的方法:计算所有点(直线的起点和终点+要拆分的点)沿直线的距离,按这些点排序,然后按正确的顺序生成线段。一起发挥作用:

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

我喜欢乔里斯的做法。不幸的是,我在尝试使用它时遇到了一个关键的困难:如果字符串中有两个点位于同一坐标,那么它们的投影就不明确了。两者都将得到相同的投影值并进行排序。

如果你有一条从同一点开始和结束的路径,这一点尤为明显。结束点得到一个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 .

返回值是 MultiLineString ,您可以找到组件 LineString s正在使用它 geoms 方法

这是我尝试通过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的方法很好,但如果你试图传递一个点列表,其中一些点实际上并不与直线相交,就会出错,在我的例子中,这是因为我从一个多条直线的列表中预先计算出一个相交点列表。

我能够通过在继续执行该函数之前将输入点列表预过滤为仅实际相交的点来解决这个问题。对于大的点列表来说,它不会有效,但在我的情况下,我的列表总是很小,所以它对我来说已经足够好了。如果没有点与线相交,它也会起作用,在这种情况下,它会短路返回原始线作为列表(为了一致性)

我最初使用 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