dangerous Streets?

[emantzo]

Hello and thank you for the insightful tools. I am not quite familiar with mapping tools in R, and I hope my question is short and relevant. Is there an easy way to group accidents by street, or otherwise extract the info on street name by accident? I would like to investigate which road sections are particularly dangerous.

Many thanks!

[Robinlovelace]

Yes, you can do this. You need to do a geographic join between the streets in the OSM data and the crashes. See ?sf::st_join for details and let us know how you get on!

[emantzo]

Thank you for your reply!

I have finally managed to proceed with the code a bit (using code from this and from related packages), and I am copying here, in case you have the time to have a look or in case anyone else might be trying a similar approach.

In the following part, I try to join crashes and osm roads, in Leeds. It seems that st_join will not work, so I am using st_nearest_feature

library(stats19)
library(ggmap)
library(dplyr)
library(tibble)
library(tidygeocoder)
library(sf)

pkgs = c(
  "pct",
  "osmextract",
  "tmap",
  "stplanr",
  "od"
)

lapply(pkgs, library, character.only = TRUE)[length(pkgs)]

tmap_mode("view")

## roads from osm

# leeds+buffer

leeds_point = tmaptools::geocode_OSM( "leeds")
c_m_coordiantes = rbind( leeds_point$coords,leeds_point$coords)
c_m_od = od::points_to_od(p = c_m_coordiantes, interzone_only = TRUE)
c_m_desire_line = od::odc_to_sf(c_m_od[-(1:2)])[1, ]
leeds_buffer = stplanr::geo_buffer(c_m_desire_line, dist = 5000)
qtm(leeds_buffer)

### osm roads

library(osmextract)
library(osmdata)

osm_lines = oe_get(leeds_buffer, stringsAsFactors = FALSE, quiet = TRUE)

## crash data
crashes = get_stats19(year = 2019, output_format = "sf", lonlat = TRUE)
casualties = get_stats19(year = 2019, type = "casualties")
crashes_combined = inner_join(crashes, casualties)
# table(crashes_combined$casualty_type)
crashes_active = crashes_combined 

crashes_in_area = crashes_active[leeds_buffer, ]

tm_shape(crashes_in_area) +
  tm_dots("accident_severity", popup.vars = c("casualty_type", "accident_severity", "datetime"), palette = "viridis")

cras_samp=crashes_in_area

points.sf <- st_as_sf(cras_samp, coords = c("longitude","latitude"),crs=4326)

## match crashes to osm roads
cras_samp$road=st_nearest_feature(points.sf$geometry,osm_lines)

dist = st_distance(points.sf$geometry,osm_lines[cras_samp$road,], by_element=TRUE)

## using st_join... not working??
dist2=st_join(st_sf(points.sf),osm_lines)

In addition, I tried to join crashes with od routes, in a similar manner (to find out whether more crashes are related to more road usage), and I am not sure if it is ok.

### od data

# get nationwide OD data
od_all <- pct::get_od("west-yorkshire")

centroids_all <- pct::get_centroids_ew() %>% sf::st_transform(4326)

leeds <- pct::pct_regions %>% filter(region_name == "west-yorkshire")
centroids_leeds <- centroids_all[leeds, ]
od_leeds <- od_all %>%
  filter(geo_code1 %in% centroids_leeds$msoa11cd) %>%
  filter(geo_code2 %in% centroids_leeds$msoa11cd)

###

desire_lines <- od2line(od_leeds, centroids_leeds)

desire_lines$distance = as.numeric(st_length(desire_lines))
desire_carshort = dplyr::filter(desire_lines, car_driver > 300 & distance < 5000)

# make routes from od
route_carshort = route(l = desire_carshort, route_fun = route_osrm)

## join routes- osm
route_osm=st_join(route_carshort, osm_lines)

## correct??
## join routes- crash
cras_samp$road2=st_nearest_feature(points.sf$geometry, route_osm)

Many thanks for your time!!

[Robinlovelace]

Hi @emantzo I'm sorry but just got back from holidays and lots of other things going on. Looks like you have made great progress, will try to pick up on this at some point. A cursory glance suggests that you're in the right ballpark, interested to see what the results look like.