How can I in GeoPandas (!) get the crossroads for a GeomDataFrame containing multilinestrings ? I have tried the solution as presented in

Intersecting lines to get points in geopandas

However, that is not working for me.

I can do it, using ST in PostgreSQL queries and a cursor, but that is a very inefficient solution.

• Can you share some sample data? Oct 14, 2021 at 19:40
• you should be able to download a limited sample dataset from github.com/EricBHK/test Oct 14, 2021 at 20:52

The problem with the `itertools.combination` solution of Kartograaf is that there are duplicate points in the result. The figure above obtained from the file provided shows 5 points but:

``````points = []
for line1, line2 in combinations([line for line in multiline],2):
if line1.intersects(line2):
points.append(line1.intersection(line2))
len(points)
11 # the result is 11 points
``````

Duplicates points (script from python count duplicate in list):

`````` for pt in [elem for ind, elem in enumerate(points,1) if elem not in points[ind:]]:
print("{} {}".format(points.count(pt),pt))

3 POINT (-13.0674742 8.560511)
1 POINT (-13.0815984 8.5574817)
3 POINT (-13.0846766 8.561736700000001)
1 POINT (-13.0813782 8.557494800000001)
3 POINT (-13.0714241 8.5588376)
``````

To eliminate the duplicates, create a new GeoDataFrame with the points list as the geometry columns

``````res = gpd.gpd.GeoDataFrame({'geometry':points})
res

geometry
0   POINT (-13.06747 8.56051)
1   POINT (-13.07142 8.55884)
2   POINT (-13.07142 8.55884)
3   POINT (-13.06747 8.56051)
4   POINT (-13.06747 8.56051)
5   POINT (-13.08468 8.56174)
6   POINT (-13.08468 8.56174)
7   POINT (-13.08160 8.55748)
8   POINT (-13.08468 8.56174)
9   POINT (-13.08138 8.55749)
10  POINT (-13.07142 8.55884)
``````

And eliminate the duplicates

``````res.drop_duplicates('geometry', inplace=True)
res
geometry
0   POINT (-13.06747 8.56051)
1   POINT (-13.07142 8.55884)
5   POINT (-13.08468 8.56174)
7   POINT (-13.08160 8.55748)
9   POINT (-13.08138 8.55749)
``````

The result is 5 points

Now you can save directly the result as a shapefile:

``````res.to_file("intersection_points.shp")
``````

But it is also possible to do it directly without `itertools` using `overlay`

``````shapefile = gpd.read_file("test.shp")
gdf = gpd.overlay(shapefile, shapefile, how='intersection', keep_geom_type=False)
gdf = gdf[gdf.geom_type=='Point'] # select only the points
gdf.drop_duplicates('geometry', inplace=True) # delete duplicates
gdf.geometry
1     POINT (-13.06747 8.56051)
2     POINT (-13.07142 8.55884)
9     POINT (-13.08138 8.55749)
19    POINT (-13.08468 8.56174)
23    POINT (-13.08160 8.55748)
Name: geometry, dtype: geometry
gdf.to_file("intersection_points.shp")
``````
• Thanks @gene! I'm always learning from folks like you around here. Oct 15, 2021 at 20:54

Here is a solution in python using a combination of geopandas, shapely, and fiona:

``````import geopandas as gpd
from shapely.geometry import mapping
from itertools import combinations
import fiona

#extract shapely geometries
multiline = shapefile['geometry']

#create empty list of points
points = []

#find intersections and append to points list
for line1, line2 in combinations([line for line in multiline],2):
if line1.intersects(line2):
points.append(line1.intersection(line2))

# Define a point feature geometry with one attribute
schema = {
'geometry': 'Point',
'properties': {'id': 'int'},
}

# Write a new points shapefile with fiona
with fiona.open('intersection_points.shp', 'w', 'ESRI Shapefile', schema) as c:
n = 1
for i in points:
c.write({
'geometry': mapping(i),
'properties': {'id': n},
})
n = n+1
``````

Original:

Result:

• Why use Fiona here when you can do it all with GeoPandas ?
– gene
Oct 15, 2021 at 8:20