I am using Python 3.7 with Shapely and GeoPandas.

I have a big line of 181,000 points, and would like to find all the points where the line intersects itself. It does so a lot. I don't need a new point at the precise intersection, just one of the existing points which is closest.

I have been writing code to loop through the points and find other points close by using.

for i,point in gdf.iterrows():
    gdf[gdf.geometry.intersects(point.buffer(10) == True].index.tolist()

Where gdf is a geopandas GeoDataFrame where each row is a point from the line. (eg it looks like this:)

0  POINT (-47.91000 -15.78000)
1  POINT (-47.92000 -15.78000)

But surely there is a way to do this using existing functions?

My way is very slow and records many duplicates at each intersection, so will require more code to reduce each intersection to one point.

2 Answers 2


update 2021:

a more elegant way using unary_union and linemerge. you can download the notebook here.

  1. read the file
import geopandas as gpd

# before
gdf = gpd.read_file('selfintersects.geojson')

enter image description here

  1. let's check the endpoints
def get_endpoints(gdf):
    from shapely.geometry import Point
    startpoint = gdf.geometry.apply(lambda x: x.coords[0])
    endpoint = gdf.geometry.apply(lambda x: x.coords[-1])

    startpoints = [Point(i) for i in startpoint]
    endpoints = [Point(i) for i in endpoint]

    return startpoints, endpoints

def create_endpoints(startpoints, endpoints):
    geom = []
    for a,b in zip(startpoints, endpoints):
        from shapely.geometry import Point

    endpoints = gpd.GeoDataFrame({'id': range(0, len(geom))}, crs=gdf.crs, geometry=geom)
    return endpoints

startpoints, endpoints = get_endpoints(gdf)
endpoints = create_endpoints(startpoints, endpoints)
import matplotlib.pyplot as plt

fig, ax = plt.subplots()

enter image description here

  1. union it to merge all lines into one geometry. Note: unary_union will take time if your data is large!
union_geom = gdf.unary_union
union = gpd.GeoDataFrame({'id':[0]}, crs=gdf.crs, geometry=[gdf.unary_union])
  1. and then explode it!
from shapely.ops import linemerge

lm = gpd.GeoDataFrame({'id':[0]}, crs=gdf.crs, geometry=[linemerge(union_geom)]).explode().reset_index(drop=True)
  1. let's check the endpoint of the exploded union.
startpoints, endpoints = get_endpoints(lm)
endpoints = create_endpoints(startpoints, endpoints)

# cleansing with snap
from shapely.ops import snap
endpoints['geometry'] = endpoints.geometry.apply(lambda x: snap(x, union_geom, 0.00001))

fig, ax = plt.subplots()

enter image description here

  1. filter out the dangles
sjoin = endpoints.sjoin(gdf, how='left')

fig, ax = plt.subplots()

enter image description here

There you go! now we have the points.

DEPRECATED answer from 2020:

Here's how I did it

  1. slice the first feature
  2. make a unary_union of the rest of the feature
  3. do line intersections using shapely
  4. you'll get one point of intersection.
  5. now repeat for the second, third, fourth, and so on.

here's the example.

  • suppose a geodataframe (gdf) of 6 lines like this GeoJSON

  • then, apply this code to the gdf. This is returning the geometry of the intersections

# the points of intersections will be appended here
for i in gdf.id:
    # check overlap
    feature = gdf[gdf['id']==i]['geometry'][i]
    overlap_feature = gdf[gdf['id']!=i]['geometry'].unary_union
    intersects = feature.intersection(overlap_feature)
  • now, make a GeoDataFrame out of the points

intersections = gpd.GeoDataFrame(
    {"id": [n for n,i in enumerate(points)]},
  • here's the plot of the result

import matplotlib.pyplot as plt
fig,ax = plt.subplots()
intersections.plot(color="r", ax=ax,zorder=2)

enter image description here

the intersections data frame has Point and MultiPoint geometries. But there's a problem here... the points are intersecting. here's how to delete the overlapping points

from shapely.geometry import Point

# convert the multipoints into points 
intersections['ispoint'] = intersections['geometry'].apply(lambda x: isinstance(x, Point)) #backup
is_point = intersections[intersections.ispoint] #check if it's point
was_multipoint = intersections[~intersections.ispoint].explode().reset_index() # converting the multipoint into points 

# now appending both data frames.
now_point = is_point.append(was_multipoint)
now_point = now_point[['id','geometry']]
now_point['id'] = now_point.index
# ok, now_point contains all intersections, but the points are still overlapping each other

# delete overlapping points
intersections2 = now_point.copy()
n= 0
for i in intersections2.id:
    # check overlap
    feature = intersections2[intersections2['id']==i]['geometry'][i]
    overlap_feature = intersections2[intersections2['id']!=i]['geometry'].unary_union

    # IF the point is intersecting with other points, delete the point!
    if feature.intersects(overlap_feature):
        intersections2.drop(i, inplace=True)
    print(n, feature.intersects(overlap_feature))

the result is the same, but the intersection points won't overlap each other. here's the plot, and there are 6 row of dataframe, I checked.

edit: note, using `unary_union` means that if we have a large dataset, this may be RAM consuming.

enter image description here

  • Nice! Thanks, I ended up using a method like this for another problem. For my original problem however, I wanted to return one of the points closest the the intersection from the original dataset, not a new point at the intersection. This is because I have elevation (z) data at those original points, which i want to compare at the cross points.
    – Aaron
    Commented Aug 30, 2020 at 9:21
  • sounds like the nearest neighbor analysis. identifying the nearest point, and joining them. look at here
    – sutan
    Commented Aug 30, 2020 at 10:31

This splits each line at each vertice and then use crosses on all combinations of split lines:

import geopandas as gpd
from shapely.geometry import LineString
from itertools import combinations
from collections import defaultdict

df = gpd.read_file('/home/bera/Desktop/crossing_lines.shp')

def findcrossings(line1, line2):
    crossings = []
    if line1.crosses(line2):
    return crossings

results = defaultdict(list)
indices = df.index.to_list()
for i in indices: #For each row/feature/line
    line = df.geometry.iloc[i] #Fetch the geometry
    verts = [v for v in line.coords] #List all vertices
    segments = []
    for p1, p2 in zip(verts, verts[1:]): #For each pair of vertices (x1,y1) , (x2,y2)
        segment = LineString([p1, p2]) #create a line segment
    for l1, l2 in combinations(segments,2): #For all combinations of segments
        res = findcrossings(l1, l2) #Find crossings
        if len(res)>0:

#            {0: [<shapely.geometry.point.Point at 0x7f5c83f356d0>,
#              <shapely.geometry.point.Point at 0x7f5c83f35990>,
#              <shapely.geometry.point.Point at 0x7f5c83f2be50>,
#              <shapely.geometry.point.Point at 0x7f5c83f2b650>]})
#Export results
#gdf = gpd.GeoDataFrame(geometry=results[0])

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.