geopandas: match accident point to nearest road

Question

I have the following GeoDataFrames: accidents_c = collection of all accidents in austria as points, streets = roads in austria as a linestrings (open street map data) Both have the crs type epsg:3310

Now i want to match every accident to the nearest road. My first attempt was this:

def nearest_street(accident_point, streets):
  row_canditates=streets.copy()
  nearest_road = None
  min_distance = None
  row_canditates["distance_road"] = row_canditates.apply(lambda row:  accident_point["geometry"].distance(row.geometry),axis=1)
  min_distance = row_canditates["distance_road"].min()
  min_road = row_canditates.loc[row_canditates["distance_road"] == min_distance]
  nearest_road = min_road["osm_id"].values[0]
  return nearest_road, min_distance

accidents_c["nearest_road"], accidents_c["distance_road"] = zip(*accidents_c.apply(nearest_street, streets=roads, axis=1))

This works but takes forever. So i was thinking about a way to make it faster by only including roads which are not more than 3000 meters away from the accident point. For this i used the buffer method. And did this:

def nearest_street(accident_point, streets):
  row_canditates=streets.copy()
  nearest_road = None
  min_distance = None
  buffered_accident = accident_point["geometry"].buffer(2000)
  bounds = buffered_accident.bounds
  x_min, x_max, y_min, y_max  = buffered_accident.bounds
  row_canditates=row_canditates.cx[x_min:x_max, y_min:y_max]

  row_canditates["distance_road"] = row_canditates.apply(lambda row:  accident_point["geometry"].distance(row.geometry),axis=1)
  min_distance = row_canditates["distance_road"].min()
  min_road = row_canditates.loc[row_canditates["distance_road"] == min_distance]
  nearest_road = min_road["osm_id"].values[0]
  return nearest_road, min_distance

accidents_c["nearest_road"], accidents_c["distance_road"] = zip(*accidents_c.apply(nearest_street, streets=roads, axis=1))

This works much faster but the results are worse. Differences are sometimes over 1000 meters between code a and code b. Where do you think is the problem in the code? Do you know any better method to limit the search for the nearest environment?

user_9441404 · Answer

I would expect https://osmnx.readthedocs.io/en/stable/osmnx.html#osmnx.distance.nearest_edges would perform better. Below shows full example of using UK accident data and a UK city. Same approach would work for Austria. Just need reference to accident data and choose a polygon that is area you want to use.

import osmnx as ox
import pandas as pd
import geopandas as gpd

place = "Gloucester"

# get bounding polygon of investigated location
gdf_poly = ox.geocode_to_gdf({"city": place}).loc[
    :, ["geometry", "display_name"]
]

# all uk accidents
df = pd.read_csv(
    "https://data.dft.gov.uk/road-accidents-safety-data/dft-road-casualty-statistics-accident-provisional-mid-year-unvalidated-2021.csv"
)
# accidents within investigated location
gdf = gpd.GeoDataFrame(
    data=df,
    geometry=gpd.points_from_xy(
        df["location_easting_osgr"], df["location_northing_osgr"]
    ),
    crs="EPSG:27700",
).to_crs("epsg:4326").sjoin(gdf_poly).reset_index(drop=True)

# OSMNX graph for investigated location
G = ox.graph_from_polygon(gdf_poly.iloc[0,0], network_type="drive")

# for speed project everything to UTM CRS
G_proj = ox.project_graph(G)
gdf = gdf.to_crs(G_proj.graph["crs"])

# get nodes and edges associated with investigated location
gdf_nodes, gdf_edges = ox.utils_graph.graph_to_gdfs(G_proj)

# find nearest edges (road) to accident points
ne, d = ox.nearest_edges(
    G_proj, X=gdf.geometry.x.values, Y=gdf.geometry.y.values, return_dist=True
)

# reindex accidents by OSM nearest edge
gdf = (
    gdf.set_index(pd.MultiIndex.from_tuples(ne, names=["u", "v", "key"]))
    .assign(distance=d)
    .sort_index()
)

# join accidents to nearest edge, now we have road name etc
gdf.join(gdf_edges.loc[:,["ref","name","highway","maxspeed"]]).loc[:,[
 'accident_year',
 'accident_reference',
 'date',
 'first_road_class',
 'first_road_number',
 'road_type',
 'speed_limit',
 'trunk_road_flag',
 'geometry',
 'ref',
 'name',
 'highway',
 'maxspeed']]

user_6874682 · Answer

i think i found the mistake.

x_min, x_max, y_min, y_max = buffered_accident.bounds

needs to be this one

x_min, y_min, x_max, y_max = buffered_accident.bounds

This was the reason why results were so bad. But thanks for your help Rob Raymond.

geopandas: match accident point to nearest road

2 Answers