Proactive tasking map development

[WillSERP]

Ideas, opinions and advice sought. I've produced a leaflet map showing various measure to help guide proactive deployment of road safety resources. Looking to add further functionality though shiny but a previous project ground to a halt when I was unable to get the shiny app to publish on shiny.io. Looking at using RInno to create an executable file from the R-shiny app that can be shared and used offline. Other development ideas include adding peak time analysis graphs showing hourly rates for each day of the week for each collision sub type.
Current output is too big to upload so here is a link https://essexhighways-my.sharepoint.com/:u:/g/personal/william_cubbin_essexhighways_org/EZqA1RQLcdlIrEmzqpZg8RgBuKBYNY2-nAxmqLXkVsKcQw?e=VKp5LI

library(tidyverse)
#> Warning: package 'tidyverse' was built under R version 3.5.3
#> Warning: package 'ggplot2' was built under R version 3.5.3
#> Warning: package 'tibble' was built under R version 3.5.3
#> Warning: package 'tidyr' was built under R version 3.5.3
#> Warning: package 'readr' was built under R version 3.5.3
#> Warning: package 'purrr' was built under R version 3.5.3
#> Warning: package 'dplyr' was built under R version 3.5.3
#> Warning: package 'stringr' was built under R version 3.5.3
#> Warning: package 'forcats' was built under R version 3.5.3
library(sf)
#> Warning: package 'sf' was built under R version 3.5.3
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(stats19)
#> Warning: package 'stats19' was built under R version 3.5.3
#> Data provided under OGL v3.0. Cite the source and link to:
#> www.nationalarchives.gov.uk/doc/open-government-licence/version/3/
library(leaflet)
#> Warning: package 'leaflet' was built under R version 3.5.3
library(lubridate)
#> Warning: package 'lubridate' was built under R version 3.5.3
#> 
#> Attaching package: 'lubridate'
#> The following object is masked from 'package:base':
#> 
#>     date
library(rgdal)
#> Warning: package 'rgdal' was built under R version 3.5.3
#> Loading required package: sp
#> Warning: package 'sp' was built under R version 3.5.3
#> rgdal: version: 1.4-6, (SVN revision 841)
#>  Geospatial Data Abstraction Library extensions to R successfully loaded
#>  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
#>  Path to GDAL shared files: C:/Users/william.cubbin/Documents/R/R-3.5.2/library/rgdal/gdal
#>  GDAL binary built with GEOS: TRUE 
#>  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
#>  Path to PROJ.4 shared files: C:/Users/william.cubbin/Documents/R/R-3.5.2/library/rgdal/proj
#>  Linking to sp version: 1.3-1
library(dplyr)
library(sp)
library(data.table)
#> Warning: package 'data.table' was built under R version 3.5.3
#> 
#> Attaching package: 'data.table'
#> The following objects are masked from 'package:lubridate':
#> 
#>     hour, isoweek, mday, minute, month, quarter, second, wday,
#>     week, yday, year
#> The following objects are masked from 'package:dplyr':
#> 
#>     between, first, last
#> The following object is masked from 'package:purrr':
#> 
#>     transpose
library(KernSmooth)
#> Warning: package 'KernSmooth' was built under R version 3.5.3
#> KernSmooth 2.23 loaded
#> Copyright M. P. Wand 1997-2009
library(geofabric)
library(raster)
#> Warning: package 'raster' was built under R version 3.5.3
#> 
#> Attaching package: 'raster'
#> The following object is masked from 'package:data.table':
#> 
#>     shift
#> The following object is masked from 'package:dplyr':
#> 
#>     select
#> The following object is masked from 'package:tidyr':
#> 
#>     extract

##### DOWNLOAD ROADS DATA ##########
#get osm data.... only necessary if you need to update roads layer
#once downloaded save to working directory using #saveRDS(roads, file = "roads.rds")#
##remotes::install_github("itsleeds/geofabrik")
##region_name = "essex"
##osm_data = get_geofabrik(name = region_name)
##head(osm_data)
##head(osm_data$geometry)
##table(osm_data$highway)
##roads = osm_data %>%
  ##filter(highway %in% c("motorway", "primary", "secondary", "tertiary", "trunk", 
                        ##"motorway_link", "primary_link", "secondary_link", "tertiary_link"))

############# ENDS #############


######## IMPORT DATA ##########

#Set working directory
setwd("C:/Users/william.cubbin/Desktop/R-projects/Routes_project")

#import AccsMap data from current working directory.
#Currently must be accidents table in standard AccsMap format saved as csv, 
#but with young driver count added to column 15.
#Look at importing all 3 tables as a .xls and getting R to filter accidents on young driver involvement based on vehs tab
Rawdata <- read.csv("AccsMap.csv", header = TRUE)

#new table with fields of interest only
Crashes <- Rawdata[,c(1,2,4,5,6,7,8,9,10,11,12,13,14,15,16,17,21,22,28,52,56,60,64,68,72)]

#import crime data
Crime <- read.csv("Crime.csv", header = TRUE)

#import STORM data
Storm <- read.csv("STORM.csv", header = TRUE)

#prepare data for heatmapping
#convert collision data from eastings and northing to lat and long without creating multipoint vector
latlongcoords = (cbind(Crashes$easting , Crashes$northing))
latlong = SpatialPointsDataFrame(latlongcoords,data = Crashes ,proj4string = CRS("+init=epsg:27700"))
latlong_SP = spTransform(latlong,CRS("+init=epsg:4326"))
heatmap_all = as.data.table(latlong_SP)


############ HEAT MAPPING ROUTINE #####################
#make contours (coords.x1 = longitude, coords.x2 = latitude)
kde <- bkde2D(heatmap_all[ , list(coords.x1, coords.x2)],
              bandwidth=c(.0035, .0035), gridsize = c(3000,3000))  
KDE_all <- raster(list(x=kde$x1 , y=kde$x2 , z=kde$fhat))

#set low density cells as NA so we can make them transparent with the colorNumeric function
KDE_all@data@values[which(KDE_all@data@values < 1)] <- NA

#create pal function for coloring the raster
palRaster <- colorNumeric("Reds", domain = KDE_all@data@values,
                          na.color = "transparent")


#repeat heatmapping for data filtered by:
#Pedestrians, cyclists, P2Ws, children, OAPS, Speed, drink/drugs

heatmap_p2w <- as.data.table(heatmap_all %>%
  filter(p2w > 0))

heatmap_cycles <- as.data.table(heatmap_all %>%
  filter(cycles > 0))

heatmap_pedestrian <- as.data.table(heatmap_all %>%
  filter(pedestrian > 0))

heatmap_YCD <- as.data.table(heatmap_all %>%
  filter(YCD > 0))

heatmap_DDI <- as.data.table(heatmap_all %>%
  filter(CF1 == "501. Impaired by alcohol" |
           CF2 == "501. Impaired by alcohol" | 
           CF3 == "501. Impaired by alcohol" |
           CF4 == "501. Impaired by alcohol" |
           CF5 == "501. Impaired by alcohol" |
           CF6 == "501. Impaired by alcohol" |
           CF1 == "502. Impaired by drugs (illicit or medicinal)" |
           CF2 == "502. Impaired by drugs (illicit or medicinal)" | 
           CF3 == "502. Impaired by drugs (illicit or medicinal)" |
           CF4 == "502. Impaired by drugs (illicit or medicinal)" |
           CF5 == "502. Impaired by drugs (illicit or medicinal)" |
           CF6 == "502. Impaired by drugs (illicit or medicinal)" ))


#P2W heatmap
kdep2w <- bkde2D(heatmap_p2w[ , list(coords.x1, coords.x2)],
              bandwidth=c(.0035, .0035), gridsize = c(3000,3000))  
KDE_p2w <- raster(list(x=kdep2w$x1 , y=kdep2w$x2 , z=kdep2w$fhat))
KDE_p2w@data@values[which(KDE_p2w@data@values < 1)] <- NA
palRaster_p2w <- colorNumeric("Reds", domain = KDE_p2w@data@values,
                          na.color = "transparent")


#Cyclist heatmap
kdecycles <- bkde2D(heatmap_cycles[ , list(coords.x1, coords.x2)],
                 bandwidth=c(.0035, .0035), gridsize = c(3000,3000))  
KDE_cycles <- raster(list(x=kdecycles$x1 , y=kdecycles$x2 , z=kdecycles$fhat))
KDE_cycles@data@values[which(KDE_cycles@data@values < 1)] <- NA
palRaster_cycles <- colorNumeric("Reds", domain = KDE_cycles@data@values,
                              na.color = "transparent")


#Pedestrian heatmap
kdepeds <- bkde2D(heatmap_pedestrian[ , list(coords.x1, coords.x2)],
                 bandwidth=c(.0035, .0035), gridsize = c(3000,3000))  
KDE_peds <- raster(list(x=kdepeds$x1 , y=kdepeds$x2 , z=kdepeds$fhat))
KDE_peds@data@values[which(KDE_peds@data@values < 1)] <- NA
palRaster_peds <- colorNumeric("Reds", domain = KDE_peds@data@values,
                              na.color = "transparent")


#Young car driver heatmap
kdeYCD <- bkde2D(heatmap_YCD[ , list(coords.x1, coords.x2)],
                 bandwidth=c(.0035, .0035), gridsize = c(3000,3000))  
KDE_YCD <- raster(list(x=kdeYCD$x1 , y=kdeYCD$x2 , z=kdeYCD$fhat))
KDE_YCD@data@values[which(KDE_YCD@data@values < 1)] <- NA
palRaster_YCD <- colorNumeric("Reds", domain = KDE_YCD@data@values,
                              na.color = "transparent")


#Drink/drug impaired heatmap
kdeDDI <- bkde2D(heatmap_DDI[ , list(coords.x1, coords.x2)],
                 bandwidth=c(.0035, .0035), gridsize = c(3000,3000))  
KDE_DDI <- raster(list(x=kdeDDI$x1 , y=kdeDDI$x2 , z=kdeDDI$fhat))
KDE_DDI@data@values[which(KDE_DDI@data@values < 1)] <- NA
palRaster_DDI <- colorNumeric("Reds", domain = KDE_DDI@data@values,
                              na.color = "transparent")


#Crime layer heat mapping
latlongcds_crime = (cbind(Crime$Easting , Crime$Northing))
latlong_crime = SpatialPointsDataFrame(latlongcds_crime, data = Crime ,proj4string = CRS("+init=epsg:27700"))
latlong_SP_crime = spTransform(latlong_crime,CRS("+init=epsg:4326"))
heatmap_crime = as.data.table(latlong_SP_crime)

kdecrime <- bkde2D(heatmap_crime[ , list(coords.x1, coords.x2)],
                   bandwidth=c(.0035, .0035), gridsize = c(3000,3000))  
KDE_crime <- raster(list(x=kdecrime$x1 , y=kdecrime$x2 , z=kdecrime$fhat))

KDE_crime@data@values[which(KDE_crime@data@values < 1)] <- NA

palRaster_crime <- colorNumeric("Reds", domain = KDE_crime@data@values,
                                na.color = "transparent")


############# ENDS #############


############# COLLISION DOT MAPPING ROUTINE############
#Define geometry
Crashes = st_as_sf (Crashes, coords = c("easting", "northing"), crs = 27700)

#add text field for severity description
severity = c(1, 2, 3)
sev_text = c("Fatal", "Serious", "Slight")
severity_lookup = data.frame(severity, sev_text)
Crashes = Crashes %>%
    left_join (severity_lookup, by = "severity") 

#add formatted date field
Crashes$DateVal = dmy(Crashes$Date)

#Colour coding for collision markers
severity <- Crashes
getColor <- function(Crashes) {
  sapply(Crashes$sev_text, function(sev_text) {
    if(sev_text == "Fatal") { "red" } 
    else if(sev_text == "Serious") { "blue" }
    else if(sev_text == "Slight") { "green" }
    else { "white" }
  })
}

#create labels
Crashes$Label <- paste("Date:", Crashes$Date,
                       "| Time:", Crashes$Time,
                       "| No. Vehicles:", Crashes$vehicles,
                       "| No. Casualties:", Crashes$casualties,
                       "| Pedestrians:", Crashes$pedestrian,
                       "| Children:", Crashes$children,
                       "| OAPs:", Crashes$oaps,
                       "| Cyclists:", Crashes$cycles,
                       "| P2Ws:", Crashes$p2w,
                       "| CF1:", Crashes$CF1,
                       "| Weather:", Crashes$Weather,
                       "| Speed limit:", Crashes$speed_lim,"mph",
                       "| Junction:", Crashes$Junct_Det)


#Copy Crashes and filter to show just speed related
Speed <- Crashes %>%
  filter(CF1 == "306. Exceeding speed limit" |
           CF2 == "306. Exceeding speed limit" | 
           CF3 == "306. Exceeding speed limit" |
           CF4 == "306. Exceeding speed limit" |
           CF5 == "306. Exceeding speed limit" |
           CF6 == "306. Exceeding speed limit" |
           CF1 == "808. Careless/Reckless/In a hurry" |
           CF2 == "808. Careless/Reckless/In a hurry" | 
           CF3 == "808. Careless/Reckless/In a hurry" |
           CF4 == "808. Careless/Reckless/In a hurry" |
           CF5 == "808. Careless/Reckless/In a hurry" |
           CF6 == "808. Careless/Reckless/In a hurry" |
           CF1 == "307. Travelling too fast for conditions" |
           CF2 == "307. Travelling too fast for conditions" | 
           CF3 == "307. Travelling too fast for conditions" |
           CF4 == "307. Travelling too fast for conditions" |
           CF5 == "307. Travelling too fast for conditions" |
           CF6 == "307. Travelling too fast for conditions" |
           CF1 == "410. Loss of control" |
           CF2 == "410. Loss of control" | 
           CF3 == "410. Loss of control" |
           CF4 == "410. Loss of control" |
           CF5 == "410. Loss of control" |
           CF6 == "410. Loss of control" |
           CF1 == "601. Aggressive driving" |
           CF2 == "601. Aggressive driving" | 
           CF3 == "601. Aggressive driving" |
           CF4 == "601. Aggressive driving" |
           CF5 == "601. Aggressive driving" |
           CF6 == "601. Aggressive driving")


#Copy Crashes and filter to show just P2W
P2W <- Crashes %>%
  filter(p2w > 0)

#Copy Crashes and filter to show just cyclist
cycles <- Crashes %>%
  filter(cycles > 0)

#Copy Crashes and filter to show just pedestrian
pedestrian <- Crashes %>%
  filter(pedestrian > 0)

#Copy Crashes and filter to show just young driver
YCD <- Crashes %>%
  filter(YCD > 0)

#Copy Crashes and filter to show Drink/Drug Impaired
DDI <- Crashes %>%
  filter(CF1 == "501. Impaired by alcohol" |
           CF2 == "501. Impaired by alcohol" | 
           CF3 == "501. Impaired by alcohol" |
           CF4 == "501. Impaired by alcohol" |
           CF5 == "501. Impaired by alcohol" |
           CF6 == "501. Impaired by alcohol" |
           CF1 == "502. Impaired by drugs (illicit or medicinal)" |
           CF2 == "502. Impaired by drugs (illicit or medicinal)" | 
           CF3 == "502. Impaired by drugs (illicit or medicinal)" |
           CF4 == "502. Impaired by drugs (illicit or medicinal)" |
           CF5 == "502. Impaired by drugs (illicit or medicinal)" |
           CF6 == "502. Impaired by drugs (illicit or medicinal)" )
           


######## CRASH PER KM BY ROAD SECTION ##############
#import roads data from rds file
roads = readRDS(file = "roads.rds")

#convert roads to osgb format
roads_osgb = st_transform(roads, 27700)

#calculate road lengths column 
roads_osgb$length = st_length(roads_osgb$geometry)

#create 20m buffer around roads_osgb
polyroads = st_buffer(roads_osgb, dist = 20, endCapStyle = "FLAT", joinStyle = "BEVEL")

#count number of crashes in each polygon and add this as a new column in the polyroads table
polyroads$crash_count <- lengths(st_intersects(polyroads, Crashes))

#add a column for crashes per km
polyroads$crash_per_km = 1000 * polyroads$crash_count / polyroads$length  

###### ENDS #######


#repeat steps above to create crashes per km for speed related only
polyroadspeed = st_buffer(roads_osgb, dist = 20, endCapStyle = "FLAT", joinStyle = "BEVEL")
polyroadspeed$crash_count <- lengths(st_intersects(polyroadspeed, Speed))
polyroadspeed$crash_per_km = 1000 * polyroadspeed$crash_count / polyroadspeed$length  
polyroadspeed <- st_transform(polyroadspeed, 4326)


#repeat steps above to create crashes per km for storm data
Storm = st_as_sf (Storm, coords = c("Easting", "Northing"), crs = 27700)
polyroadstorm = st_buffer(roads_osgb, dist = 20, endCapStyle = "FLAT", joinStyle = "BEVEL")
polyroadstorm$count <- lengths(st_intersects(polyroadstorm, Storm))
polyroadstorm$inc_per_km = 1000 * polyroadstorm$count / polyroadstorm$length  
polyroadstorm <- st_transform(polyroadstorm, 4326)



##### PREPARE DATA FOR LEAFLET #########

#set CRS for all map elements to Lat/Long so they work with leaflet
polyroads_plot <- st_transform(polyroads, 4326)
polyroadspeed <- st_transform(polyroadspeed, 4326)
Crashes <- st_transform(Crashes, 4326)
Speed <- st_transform(Speed, 4326)
P2W <- st_transform(P2W, 4326)
cycles <- st_transform(cycles, 4326)
pedestrian <- st_transform(pedestrian, 4326)
YCD <- st_transform(YCD, 4326)
DDI <- st_transform(DDI, 4326)



#define map key for road lengths, select bins to give best ability to discriminate between different road sections
pal = colorBin("RdYlBu", domain = polyroads_plot$crash_per_km, bins = c(0,1,2,5,8,12,20,50,100,500,2000), 
               pretty = TRUE, reverse = TRUE)

#repeat for STORM data
palstorm = colorBin("RdYlBu", domain = polyroadstorm$inc_per_km, bins = c(0,1,2,5,8,12,20,50,100,500,2000), 
               pretty = TRUE, reverse = TRUE)


### DEV IDEAS #####
#Add date slider if possible? 
#OR add interactive features using shiny and compile to a sharable app using RInno
#Output peak time analysis graphs for various selected subsets
#Create function to determine distance from each collision to nearest crime within +-1 hour of time.
###########


###### LEAFLET MAP ROUTINE #######
map <- leaflet(data = polyroads_plot) %>%
    addTiles() %>%
    setView(lng = 0.58, lat = 51.78 , zoom = 10) %>%
  
    addPolygons(col = ~pal(crash_per_km), opacity = 0.8, group = "Crashes per km (all)") %>%
    addPolygons(data = polyroadspeed, col = ~pal(crash_per_km), opacity = 0.8, group = "Crashes per km (speed)") %>%
    addPolygons(data = polyroadstorm, col = ~palstorm(inc_per_km), opacity = 0.8, group = "STORM road incs per km") %>%
  
    addCircleMarkers(data = Crashes, radius = 7, 
                     color = getColor(severity), popup = ~htmlEscape(Label),
                     group = "All Collisions") %>%
    addCircleMarkers(data = Speed, radius = 7, 
                    color = getColor(severity), popup = ~htmlEscape(Label),
                    group = "Speed related collisions") %>%
  
    addCircleMarkers(data = P2W, radius = 7, 
                   color = getColor(severity), popup = ~htmlEscape(Label),
                   group = "P2W collisions") %>%
    addCircleMarkers(data = cycles, radius = 7, 
                   color = getColor(severity), popup = ~htmlEscape(Label),
                   group = "Cyclist collisions") %>%
    addCircleMarkers(data = pedestrian, radius = 7, 
                   color = getColor(severity), popup = ~htmlEscape(Label),
                   group = "Pedestrian collisions") %>%
    addCircleMarkers(data = YCD, radius = 7, 
                   color = getColor(severity), popup = ~htmlEscape(Label),
                   group = "Young car drivers") %>%
    addCircleMarkers(data = DDI, radius = 7, 
                   color = getColor(severity), popup = ~htmlEscape(Label),
                   group = "Drink/drug driver collisions") %>%
    
  
    #add heat maps
    addRasterImage(KDE_all, colors = palRaster, opacity = .7, group = "Heatmap (all)") %>%
    addRasterImage(KDE_p2w, colors = palRaster_p2w, opacity = .7, group = "Heatmap (P2W)") %>%
    addRasterImage(KDE_cycles, colors = palRaster_cycles, opacity = .7, group = "Heatmap (Cyclists)") %>%
    addRasterImage(KDE_peds, colors = palRaster_peds, opacity = .7, group = "Heatmap (Pedestrian)") %>%
    addRasterImage(KDE_YCD, colors = palRaster_YCD, opacity = .7, group = "Heatmap (Young car drivers)") %>%
    addRasterImage(KDE_DDI, colors = palRaster_DDI, opacity = .7, group = "Heatmap (Drink/drug impaired drivers)") %>%

    addRasterImage(KDE_crime, colors = palRaster_crime, opacity = .7, group = "Heatmap (Crime)") %>%
   
    #add legends  
    addLegend("topright", colors = c("red", "blue", "green"), 
            labels = c("Fatal", "Serious", "Slight"), title = "Crash severity")%>%  
    addLegend("topright", pal = pal, values = ~crash_per_km, title = "Road incident density") %>%
    addLegend("topright", pal = palRaster, values = KDE_all@data@values, title = "Heat map density")%>%


#Layer control dialogue box
    addLayersControl(
      overlayGroups = c("Crashes per km (all)", "Crashes per km (speed)", "STORM road incs per km",
                        "Heatmap (all)", "Heatmap (P2W)", "Heatmap (Cyclists)","Heatmap (Pedestrian)",
                        "Heatmap (Young car drivers)","Heatmap (Drink/drug impaired drivers)", "Heatmap (Crime)", 
                        "All Collisions", "Speed related collisions","P2W collisions", "Cyclist collisions",
                        "Pedestrian collisions", "Young car drivers", "Drink/drug driver collisions"),
                      position = "topleft",
                      options = layersControlOptions(collapsed = FALSE,)
                      ) %>%
      hideGroup("Speed related collisions") %>%
      hideGroup("Crashes per km (speed)") %>%
      hideGroup("STORM road incs per km") %>%
      hideGroup("Heatmap (all)") %>%
      hideGroup("Heatmap (P2W)") %>%
      hideGroup ("Heatmap (Cyclists)") %>%
      hideGroup ("Heatmap (Pedestrian)") %>%
      hideGroup ("Heatmap (Young car drivers)") %>%
      hideGroup ("Heatmap (Drink/drug impaired drivers)") %>%
      hideGroup ("Heatmap (Crime)") %>%
      hideGroup("P2W collisions") %>%
      hideGroup("Cyclist collisions") %>%
      hideGroup("Pedestrian collisions") %>%
      hideGroup("Young car drivers") %>%
      hideGroup("Drink/drug driver collisions") %>%
      hideGroup("All Collisions")
#> Error in htmlEscape(Label): could not find function "htmlEscape"

#plot map
map
#> function (.x, .f, ...) 
#> {
#>     .f <- as_mapper(.f, ...)
#>     .Call(map_impl, environment(), ".x", ".f", "list")
#> }
#> <bytecode: 0x0000000018a330a8>
#> <environment: namespace:purrr>

^{Created on 2020-01-30 by the reprex package (v0.3.0)}

[Robinlovelace]

Very impressive Will, many thanks for sharing. No detailed comments now, other than that the outputs make me think that implementing a method for assigning crash points to non-overlapping road sections and junctions, as @agila5 has demonstrated in ropensci/stplanr#362, should be bumped up the priority list. In fact I had a go at implementing this method and plan to push an implementation to a branch of stplanr for testing soon. Will keep you posted!

[agila5]

Hi Will! I just checked your example and I agree with Robin's suggestions. Moreover I think you could also apply some of the methods described here to combine and split street segments. I know that the code there is quite bad but I will update it asap.

I also tried to reproduce your example but I don't have several .csv or .RDS files. I noticed the error at the end (i.e. Error in htmlEscape(Label): could not find function "htmlEscape") and I think you should just update your R version and run install.packages("htmltools").

[layik]

Hi @WillSERP, relevant and great work but I know it is not an issue for the package.

Could you possibly share the data you generated in another format? I am not good with shiny so can't (reasonable time) extract it from the file you shared. And from Andea's comment I did not have a go at running the code either.

Kindly close the issue if you think it is not an issue for the package, but we can carry on discussing your needs.