2

enter image description here

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.

2

2 Answers 2

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

enter image description here

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")
1
  • Thanks @gene! I'm always learning from folks like you around here.
    – Kartograaf
    Oct 15, 2021 at 20:54
1

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:

enter image description here

Result:

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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