This is a small task which I completed as a first repository on GitHub (https://github.com/routineactivity/hot_routes).
Also known as graduated street segments. This is a basic step-by-step on how to extract a street network from Open Street Map features and open-source crime point data, and use a spatial join to count the number of offences along small street segments. This is output as a basic leaflet map.
For further information on the background and history of this method, see the links below:
Packages required for this project
# osm for retrieving open street map roads polylines
if(!require("osmdata")) install.packages("osmdata")
# spatial features
if(!require("sf")) install.packages("sf")
# spatial data
if(!require("sp")) install.packages("sp")
# GIS Tools
if(!require("GISTools")) install.packages("GISTools")
# nice interactive maps
if(!require("leaflet")) install.packages("leaflet")
# data wrangling
if(!require("tidyverse")) install.packages("tidyverse")
# data wranglingOSM features available for geospatial analysis can be viewed here https://wiki.openstreetmap.org/wiki/Map_features
We can see which features are available, the index url above for osm identifies that we can find road data under the categories of highway.
Results show the various features contained within highway. Of interest for our crime by street segments, or hot routes, we are most likely interested in:
# see which sub-categories fall within the highway features
available_tags("highway")## [1] "bridleway" "bus_guideway" "bus_stop"
## [4] "busway" "construction" "corridor"
## [7] "crossing" "cycleway" "elevator"
## [10] "emergency_access_point" "emergency_bay" "escape"
## [13] "footway" "give_way" "living_street"
## [16] "milestone" "mini_roundabout" "motorway"
## [19] "motorway_junction" "motorway_link" "passing_place"
## [22] "path" "pedestrian" "platform"
## [25] "primary" "primary_link" "proposed"
## [28] "raceway" "residential" "rest_area"
## [31] "road" "secondary" "secondary_link"
## [34] "service" "services" "speed_camera"
## [37] "steps" "stop" "street_lamp"
## [40] "tertiary" "tertiary_link" "toll_gantry"
## [43] "track" "traffic_mirror" "traffic_signals"
## [46] "trailhead" "trunk" "trunk_link"
## [49] "turning_circle" "turning_loop" "unclassified"For this example I am using the city of Baltimore, Maryland, US.
bbox_1 <- getbb("Baltimore") %>%
# get bounding box of Baltimore
opq() %>%
# add feature
add_osm_feature("highway")
# check out the structure
str(bbox_1)## List of 4
## $ bbox : chr "39.1972328,-76.7112977,39.37204,-76.5296757"
## $ prefix : chr "[out:xml][timeout:25];\n(\n"
## $ suffix : chr ");\n(._;>;);\nout body;"
## $ features: chr " [\"highway\"]"
## - attr(*, "class")= chr [1:2] "list" "overpass_query"
## - attr(*, "nodes_only")= logi FALSEWe want a polyline of the street network in Baltimore.
# create a simple features object from bbox_1
b_md_roads <- osmdata_sf(bbox_1)
# see object classes, will show a simple features collection by number of points, lines, polygons etc
class(b_md_roads)## [1] "list" "osmdata" "osmdata_sf"# we want the polylines
b_md_roadsPL <- b_md_roads$osm_linesThe osm datasets are very wide, but we are only interested in a small number of columns for this process:
# subset data
br <- b_md_roadsPL[,c("osm_id", "name", "highway", "geometry")]
# add length of segments, useful if wanting to calculate observations per 400m/ft for example
br$length <- as.numeric(st_length(br$geometry))
# subset the highway features of interest
baltimore_roads <- subset(br, highway=="motorway" | highway=="primary" | highway=="residential" | highway=="secondary" | highway=="tertiary" | highway=="track" | highway=="trunk")We are using BPD Part 1 crime data which is available open source via url.
# import data from url
bpd <- read.csv(url("https://opendata.arcgis.com/datasets/3eeb0a2cbae94b3e8549a8193717a9e1_0.csv?outSR=%7B%22latestWkid%22%3A2248%2C%22wkid%22%3A102685%7D"))
# extract date from the datetime column
bpd$date <- as.Date(substr(bpd$CrimeDateTime,1,10))
# filter group of serious violent crimes since 2018
bpd_violence <- bpd %>%
filter(Description =="AGG. ASSAULT" | Description =="HOMICIDE" | Description =="SHOOTING") %>%
filter(date >"2018-02-28")Saving the filtered sample of data.
save(bpd_violence, file = "data_bpd_sample.RData")
# load("data_bpd_sample.RData") if opendata link changes# create simple features spatial data and include projected coordinates
bpd_violenceSF <- st_transform(st_as_sf(bpd_violence, coords = c("Longitude", "Latitude"), crs=4326), crs = 6487)
# create a buffer area around the street segments, this example is 10m using projected coordinates
# note that lat/lon degrees have been set as coordinate reference (4326), this is same crs as the crime point data
br_crime_buff <- st_transform(st_buffer(bpd_violenceSF, 10), crs = 4326)# intersect street segments intersecting buffered crime points
bpd_vio_seg_sf <- st_join(baltimore_roads, br_crime_buff, join = st_intersects)
#create table grouped by osm_id for crime counts, using Homicide only in this example
counts_by_seg <- bpd_vio_seg_sf %>%
filter(!is.na(Neighborhood)) %>%
filter(Description == "HOMICIDE") %>%
group_by(osm_id, name, highway, length) %>%
count()
#arbitrarily choosing counts of street segments with more than 1 homicide
counts_by_seg2 <- counts_by_seg %>%
filter(n > 1)summary(counts_by_seg2$n)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.000 2.000 2.000 2.831 3.000 9.000# set colour palette
pal <- colorNumeric(
palette = brewer.pal(3, "Reds"),
domain = c(0, 5, 9))
# create leaflet map
leaflet() %>%
addProviderTiles("CartoDB.DarkMatter") %>%
addPolylines(
data = counts_by_seg2,
col = pal(counts_by_seg2$n),
opacity = 1) %>%
addLegend(pal = pal, values = counts_by_seg2$n, title = "Homicide", position = "bottomleft")