# How do I find where a line intersects itself?

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

``````   geometry
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.

update 2021:

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

``````import geopandas as gpd

# before
gdf.plot()
``````

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
geom.append(a)
geom.append(b)

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()
gdf.plot(ax=ax)
endpoints.plot(ax=ax)

``````

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])
union.plot()
``````
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)
lm.plot()
``````
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()
gdf.plot(ax=ax)
endpoints.plot(ax=ax)
``````

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

fig, ax = plt.subplots()
gdf.plot(ax=ax)
sjoin[sjoin['index_right'].isna()].plot(ax=ax)

``````

There you go! now we have the points.

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
points=[]
for i in gdf.id:
print(i)
# check overlap
feature = gdf[gdf['id']==i]['geometry'][i]
overlap_feature = gdf[gdf['id']!=i]['geometry'].unary_union
intersects = feature.intersection(overlap_feature)
points.append(intersects)
points
``````
• now, make a `GeoDataFrame` out of the `points`

``````intersections = gpd.GeoDataFrame(
{"id": [n for n,i in enumerate(points)]},
crs={'init':'epsg:4326'},
geometry=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)
gdf.plot(ax=ax,zorder=1)
``````

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.reset_index(inplace=True)
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()
points=[]
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))
n+=1
intersections2
``````

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.

• 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. Commented Aug 30, 2020 at 9:21
• sounds like the nearest neighbor analysis. identifying the nearest point, and joining them. look at here 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

def findcrossings(line1, line2):
crossings = []
if line1.crosses(line2):
crossings.append(line1.intersection(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
segments.append(segment)
for l1, l2 in combinations(segments,2): #For all combinations of segments
res = findcrossings(l1, l2) #Find crossings
if len(res)>0:
results[i].extend(res)

#results
#defaultdict(list,
#            {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])
#gdf.to_file('/home/bera/Desktop/crossing_lines_self_intersections.shp')
``````