Shapely Split LineStrings aux intersections avec d'autres LineStrings

Question

J'ai un ensemble de LineStrings qui sont croisées par d'autres LineStrings et je veux diviser LineString en segments séparés à ces points d'intersection. J'ai une solution mais je ne pense pas que ce soit la meilleure approche.

Disons que nous avons affaire à un 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> 
>>> 

Et 2 lignes qui se croisent ce LineString:

>>> 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 
>>> 

Il ressemble à ceci:

je peux obtenir le points d'intersection comme suit:

>>> IntPoints=MyLine.intersection(IntersectionLines.unary_union) 
>>> IntPointCoords=[x.coords[:][0] for x in IntPoints] 
>>> IntPointCoords 
[(2.75, 0.0), (5.833333333333333, 0.5)] 
>>> 

Et puis-je obtenir en extraire les points de début et de fin et créer des paires avec ces et les points d'intersection qui seront utilisés pour former des segments de ligne:

>>> 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 
>>> 

Cependant, une question que j'ai avec mon approche est que Je sais quel point d'intersection doit être joint au point de départ LineString et quel point d'intersection doit être associé au point d'extrémité LineString. Ce programme fonctionne si les points d'intersection sont listés le long de la ligne dans le même ordre que le point de départ et le point final. J'imagine qu'il y aurait des situations où ce ne serait pas toujours le cas? Y a-t-il un moyen de généraliser cette approche ou y a-t-il une meilleure approche?

Answer 1

Voici une façon plus générale: calculer la distance le long de la ligne pour tous les points (début et fin de la ligne + points où vous voulez diviser), trier par ces points, puis générer les segments de droite commande. Ensemble, dans un funtion:

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 

L'application de cette fonction sur votre exemple:

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 

La seule chose est que cela ne préserve pas le coin de votre ligne d'origine (mais votre exemple dans la question le fait ni, donc je ne sais pas si cela est une exigence. il serait possible, mais le rendre un peu plus complexe)

---

Mise à jour Une version qui maintient les coins de la ligne intacte d'origine (mon approche est de garder une liste de 0/1 qui indique si un coord doit être scindée ou non):

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 

Appliqué sur l'exemple:

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 

Answer 2

Et voici ma tentative d'adapter la fonction par joris afin que le coin du segment de ligne soit également inclus. Cela ne fonctionne pas encore parfaitement car, en plus d'inclure le segment fusionné qui inclut le coin, il inclut également le segment non fusionné d'origine.

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 
>>> 

Answer 3

J'aime l'approche de joris. Malheureusement, j'ai rencontré une difficulté critique en essayant de l'utiliser: si les lignes ont deux points aux mêmes coordonnées, leurs projections sont ambiguës.Les deux auront la même valeur de projection et seront triés ensemble.

Ceci est particulièrement évident si vous avez un chemin qui commence et se termine au même point. Le point de fin obtient une projection de 0 et est trié au début, ce qui lève l'algorithme entier puisqu'il attend une valeur de "1" à la fin.

est ici une solution qui fonctionne bien faite 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)) 

Oui, il est vraiment aussi simple que cela. Le hic ici est que les points doivent être exactement sur la ligne. Si ce n'est pas le cas, accrochez-les à la ligne comme dans this answer.

La valeur de retour est MultiLineString et vous pouvez obtenir le composant LineString s en utilisant sa méthode geoms.