GeoPandas road intersections (crossroads)

Question

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.

Answer 1

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

Answer 2

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

#read multiline with geopandas
shapefile = gpd.read_file('C:/path/to/test.shp')

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