CODEWITHSUNDEEP
pythonnetworkxopenstreetmapgeopandasosmnx
ratsrule23
ratsrule23

Reputation: 41

Finding neighbouring edges in gdf from OSM road network

I have extracted all the edges from an osm road network the following way:

#load gdf of amsterdam from osm
gdf = ox.geocode_to_gdf("Drechterland, Netherlands")

#get a graph of the union of their boundaries
G = ox.graph_from_polygon(gdf.unary_union, network_type = "drive")

#project graph
graph = ox.project_graph(G)

#get edges
nodes, edges = ox.graph_to_gdfs(graph)

The edges gdf contains features of each edge, e.g. there is a column for the highway type and max speed. I would like to have additional columns that contain those features for edges that are within a perimeter of k meters around the edge (either around the central point of the edge or, e.g. k meters around each point of the edge).

I know that it is possible to get the nearest points or edges from a certain point, but is there a way to get edges that fall within a perimeter of another edge?

Upvotes: 1

Views: 632

Answers (2)

ratsrule23
ratsrule23

Reputation: 41

Matthews answer is great to get the roads that are within the perimeter. Since my ultimate goal was to get features of the surrounding edges, I used a slightly different approach. I borrowed pieces from Jeets answer here. I generated columns that give the number of edges within the buffered area around the edge, the total length of all edges and the most common highway type.

buffered_edges = edges.buffer(1000) 

total_length = []
nr_edges = []
highway_types = []

#loop over buffered edges
for i in range(0, len(buffered_edges)):
    row = buffered_edges[i:(i+1)].unary_union #take the first buffer zone

    #find all edges that intersect this buffer zone
    intersectants = edges["geometry"].intersection(row) 
    intersecting_edges = edges[~intersectants.is_empty]
    
    #most common highway type in buffer zone
    highway = intersecting_edges["highway"].value_counts().idxmax()
    highway_types.append(highway)
    
    #total edge length in buffer zone
    Length = intersecting_edges["length"].sum()
    total_length.append(Length)
   
    #number of edges in buffer zone
    Nr = len(intersecting_edges)
    nr_edges.append(Nr)


#create columns
edges["common_highway_bufferzone"] = highway_types
edges["total_length_bufferzone"] = total_length
edges["nr_edges_bufferzone"] = nr_edges

Upvotes: 0

Matthew Borish
Matthew Borish

Reputation: 3096

This will add a new column called sjoin_name with arrays of edges within the buffer distance for each edge. In short, we're just buffering each edge and then doing a spatial join for each edge against all the buffers. This is computationally expensive, so it may take a while to run on a large data set. Change the buffer distance and crs codes as needed.

import numpy as np
import geopandas as gpd
import pandas as pd


def find_roads(row, gdf_buff, my_crs=6559):    
        
        gs = gpd.GeoSeries(row['geometry'])

        row=gpd.GeoDataFrame(geometry=gs)
        row.crs=my_crs
        
        gdf_buff.columns = ['sjoin_name', 'geometry']

        gdf_sjoin = gpd.sjoin(row, gdf_buff)

        if gdf_sjoin.shape[0] > 0:
            return gdf_sjoin['sjoin_name'].values
        else:
            return None
        
gdf = gpd.read_file('edges.shp')
gdf = gdf.to_crs(6559)

gdf_buff = gdf.copy(deep=True)

gdf_buff['geometry'] = gdf_buff['geometry'].buffer(500)

gdf['sjoin_name'] = gdf.apply(lambda row : find_roads(row, gdf_buff[['name', 'geometry']]), axis=1)

Upvotes: 1

Related Questions