downloader::download("https://geoportal.statistics.gov.uk/Docs/Boundaries/County_and_unitary_authorities_(E+W)_2012_Boundaries_(Generalised_Clipped).zip", "data/cuas.zip")
unzip("data/cuas.zip", exdir = "data")
cuas <- raster::shapefile("data/CTYUA_DEC_2012_EW_BGC.shp")
plot(cuas)

Then follow this link to find out how to get data for road traffic ‘accidents’: https://github.com/Robinlovelace/bikeR/blob/master/video-routes/load-stats19.R

After running the code in this file, the challenge is to create a dataset of traffic casualties like this:

library(spatstat)
## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-6. For overview type 'help("mgcv-package")'.
## Loading required package: deldir
## deldir 0.0-22
## 
## spatstat 1.32-0     (nickname: 'Logistical Nightmare') 
## For an introduction to spatstat, type 'beginner'
library(maptools)
## Loading required package: sp
## Checking rgeos availability: TRUE
library(rgeos)
## rgeos version: 0.3-11, (SVN revision 479)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Linking to sp version: 1.1-0 
##  Polygon checking: TRUE
library(raster)
## 
## Attaching package: 'raster'
## 
## The following objects are masked from 'package:spatstat':
## 
##     rotate, shift
## 
## The following object is masked from 'package:nlme':
## 
##     getData
# downloader::download("https://github.com/Robinlovelace/bikeR/raw/master/geodata/bikeWY.geojson", "bikeWY.geojson")
ac <- rgdal::readOGR("../bikeWY.geojson", "OGRGeoJSON")
## OGR data source with driver: GeoJSON 
## Source: "../bikeWY.geojson", layer: "OGRGeoJSON"
## with 966 features
## It has 32 fields
ac <- spTransform(ac, CRS("+init=epsg:27700"))

# estimate 2d density
acp <- as.ppp(ac)
adens <- density.ppp(x = acp, sigma = 50, eps = 50)
plot(adens)

arast <- raster(adens)
plot(arast)

writeRaster(x = arast, filename = "arast.tif", overwrite = T)
dsg <- as(arast, "SpatialGridDataFrame")
dsg <- as.image.SpatialGridDataFrame(dsg)
dcl <- contourLines(dsg, nlevels = 10)
sldf <- ContourLines2SLDF(dcl)
plot(sldf[8,]) # the most intense accident hotspot

h1 <- gPolygonize(sldf[8,])
spChFIDs(h1) <- 1
h2 <- gPolygonize(sldf[7,])
spChFIDs(h1) <- 2
plot(h2)

h3 <- gPolygonize(sldf[6,])
spChFIDs(h3) <- seq(101, 100 + length(h3))
hspots <- spRbind(h3, h2)
h4 <- gPolygonize(sldf[5,])

h5 <- gPolygonize(sldf[3,]) # the right contour to save
length(h5)
## [1] 6
proj4string(h5) <- proj4string(ac)
nacs <- aggregate(ac, h5, length)
nacs <- spTransform(nacs, CRS("+init=epsg:27700"))
nacs$area <- gArea(nacs, byid = T)
plot(nacs)
plot(adens, add = T)
plot(nacs, add = T)