Geographic Information Systems Asked by Aaron W on April 9, 2021
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.
Here's how I did it
unary_union
of the rest of the featureintersections
using shapely
here's the example.
gdf
) of 6 lines like this GeoJSONgdf
. 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
GeoDataFrame
out of the points
intersections = gpd.GeoDataFrame(
{"id": [n for n,i in enumerate(points)]},
crs={'init':'epsg:4326'},
geometry=points
)
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
Answered by sutan on April 9, 2021
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):
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')
Answered by BERA on April 9, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP