Density estimation along networks

Section 8.2 Density estimation along networks

One way of estimating the intensity of crime along spatial networks is to use kernel density estimation methods. You could, as we did in Chapter 7, use standard kernel density estimation methods, but this ignores the underlying structure in the data. As Jeremy Gelb indicates in the vignette for spNetwork, there are a number of problems with this cruder approach:

🔗

Calculating density values for locations outside the network is meaningless and the Euclidean distance underestimates the real distance between two objects on the network. Moreover, networks are not isotropic spaces. In other words, it is not possible to move in every direction, but only along the edges of the network ([Gelb, 2021]).
🔗

🔗

Smoothing through kernel density estimation in linear networks is an area of active development and a number of approaches have been proposed over recent years, some of which provide a pragmatic (heuristic) solution, but "their properties are not well understood, statistical insight is lacking," and "computational cost is high" ([McSwiggan et al., 2016], p. 325). The existing software implementation of the "traditional" algorithms that have been developed can deal with small datasets, but may become prohibitive when dealing with larger datasets. More recent algorithms aim to provide faster solutions for larger datasets without seriously compromising performance (see [Rakshit et al., 2019]). You can currently do this in R with at least two separate packages spatstat and the newer spNetwork. For reasons of space, given that it provides a more general framework for point pattern analysis and a wider selection of (more principled) approaches, we will only cover spatstat here.

🔗

In spatstat a point pattern on a linear network is represented by an object of class lpp (stands for point pattern in a linear network). The functions lpp() and as.lpp() convert raw data into an object of class lpp. The two needed arguments are a matrix or data frame with the coordinates of the points (or alternatively a ppp object that we discussed in Chapter 7) and a linear network of class linnet.

🔗

The linnet object will represent the network of a set of line segments, containing information about each of these segments (edges) and its vertices (nodes). You can generate a linnet object from a SpatialLines (sp) object using as.linnet.SpatialLines(). If you need to process this network before it is ready, you can do so with some functions from spatstat or SpnetPrep as preferable (and only then create your linnet object).

🔗

We will show here how to transform a SpatialLines object into linnet and how, then, we generate a lpp object.

🔗

# transform to projected coordinates (British National Grid)
heswall_lines <- st_transform(heswall_lines, crs = 27700)
# make into sp object
heswall_roads_sp <- as(heswall_lines, 'Spatial')
# turn into line network (linnet)
heswall_roads_linnet <- as.linnet.SpatialLines(heswall_roads_sp)

Now let’s get some crime data to link to this line network. We can find this from the open data portal data.police.uk and choose to get data from Merseyside police, who cover this area. For the sake of practice here, one year’s worth of data is included with the datasets available with this textbook (see Preamble). The file is merseyside_asb.csv.

🔗

# read in asb data
shoplifting <- read_csv("data/merseyside_shoplifting.csv")

# transform to spatial object
shoplifting_sf <- st_as_sf(shoplifting,
                           coords = c("longitude", "latitude"),
                           crs = 4326)
# transform to match projection of lines
shoplifting_sf <- st_transform(shoplifting_sf, crs = st_crs(heswall_lines))
# plot them together to have a look
ggplot() +
  geom_sf(data = heswall_lines) +
  geom_sf(data = shoplifting_sf, colour = "red") +
  theme_void()

A simple line street map of Heswall is now accompanied by red dots spreading far above and to the right as well. A concentration of red dots with a vertical gap towards the middle of the figure probably corresponds to Liverpool. — Figure 8.2.1. All shoplifting incidents in Merseyside from May 2019 to April 2020
🔗

We can see that Merseyside Police cover a lot more than just Heswall, so we can use the bounding box around the roads to subset our crime data:

🔗

# create a bounding box around the heswall lines geometry
heswall_bb <- st_bbox(heswall_lines) %>% st_as_sfc()
# subset using the bounding box
heswall_shoplifting <- st_intersection(shoplifting_sf, heswall_bb)

Now that we have the line network and the shoplifting incidents in Heswall, we can create our lpp object. First we create a ppp (point pattern) object, following the steps we learned in Chapter 7 on spatial point patterns:

🔗

# create a window
window <- as.owin(heswall_roads_linnet)
# extract coordinates
shoplifting_coords <- matrix(unlist(heswall_shoplifting$geometry), ncol = 2, byrow = T)
# make into ppp
shoplifting_ppp <- ppp(x = shoplifting_coords[,1], y = shoplifting_coords[,2],
               window = window, check = T)

Now that we have the ppp object, we can use it to create our lpp object (notice that, as noted above, there are other ways to input your points when creating a lpp object) and then plot it. There are 264 locations, some of which represent up to 6 incidents. This is likely the result of geomasking due to privacy concerns ([Tompson et al., 2015]).

🔗

# now use this and previously created linnet to make lpp object
shoplifting_lpp <- lpp(shoplifting_ppp, heswall_roads_linnet, window = window)

Now we can calculate the density and visualise. Like we did before, we can estimate the density of data points along the networks using Kernel estimation with density.lpp(). This is an area of spatstat that has seen significant development in the latest release. Depending on the arguments passed to the function, a different algorithm will be employed to estimate the density.

🔗

# calculate density
d100 <- density.lpp(shoplifting_lpp, 100, finespacing=FALSE)
# plot
plot(d100)

A pixelated street map of Heswall is titled ’d100’. A legend to the right matches a gradient from dark blue to red to yellow, with values from zero to nearly point zero eight. The lines corresponding to the streets on the map are themselves shaded, mostly dark blue. A few red areas can be seen, along with one yellow pixel. — Figure 8.2.2. Network density map with defaults
🔗

plot(d100, style="width")

The line street map of Heswall, titled ’d100’. Some of the streets seem thicker or as if they had blotches on them. To the right, a legend matches four increasing examples of line widths to the numbers zero, point zero two, point zero four, and point zero six. — Figure 8.2.3. Network density map with defaults with line width
🔗

A lpp object contains the linear network information, the spatial coordinates of the data points, and any number of columns of marks (in this case, the mark is telling us the type of crime we are dealing with). It also contains the local coordinates seg and tp for the data points. The local coordinate seg is an integer identifying the particular street segment the data point is located in. A segment is each of the sections of a street between two vertices (marking the intersection with another segment). The local coordinate tp is a real number between 0 and 1 indicating the position of the point within the segment: tp=0 corresponds to the first endpoint and tp=1 corresponds to the second endpoint.

🔗

The visual inspection of the map suggests that the intensity of crime along the network is not spatially uniform. Alternatively one could use the new specific functions for each of these different algorithms: densityQuick() (see [Rakshit et al., 2019]), densityEqualSplit() (see [Okabe and Sugihara, 2012]), or densityHeat() (see [McSwiggan et al., 2016]).

🔗

The Okabe-Sugihara method was the first to be developed for this kind of applications and can be invoked with densityEqualSplit(). We can specify the kernel (the default is epanechnikov) and the smoothing bandwidth, sigma (the standard deviation of the kernel). Above, we used 100 metres for illustration, but see Chapter 7: Spatial Point Patterns of Crime Events in this book for details on these parameters. There are two versions of this algorithm, the continuous that is excruciatingly slow (particularly for larger datasets) and the discontinuous that is a bit quicker but has less desirable statistical properties. We are using the number of events per point as the weight for computation, and wrapping everything on system.time() to provide a notion of computation time. We will only show the discontinuous version.

🔗

system.time(d_okabe <- densityEqualSplit(shoplifting_lpp,
                                         weights = marks(shoplifting_lpp),
                                         sigma = 100,
                                         continuous = FALSE,
                                         verbose= FALSE))

##    user  system elapsed 
##   9.381   0.154   9.550

And to visualise:

🔗

plot(d_okabe, main="Okabe-Sugihara (disc.)")

A pixelated street map of Heswall is titled ’Okabe-Sugihara (disc.)’. A legend to the right matches a gradient from dark blue to red to yellow, with values from zero to about point one. The lines corresponding to the streets on the map are themselves shaded, mostly dark blue. A few lighter blue areas are visible, reaching purple in a few spots. In addition, the ends of two streets to the top right brighten to yellow. — Figure 8.2.4. Network density map (Okabe-Sugihara (disc.) method)
🔗

[McSwiggan et al. (2016)] subsequently developed an algorithm that is equivalent to the Okabe-Sugihara method using the continuous rule applied to the Gaussian kernel, and therefore has more desirable properties than the equal split discontinuous, but that runs faster (than the Okabe-Sugihara equal split continuous version).

🔗

system.time(d_mcswiggan <- densityHeat.lpp(shoplifting_lpp,
                                    weights = marks(shoplifting_lpp),
                                    sigma = 100,
                                    finespacing=FALSE))

##    user  system elapsed 
##   1.085   0.007   1.093

plot(d_mcswiggan, main="McSwiggan")

Like the previous figure, we have a shaded street map of Heswall titled ’McSwiggan’, with a similar legend. Amidst the dark blue, a few areas of lighter blue appear, as well as several red to purple spots. One area towards the top right reaches yellow. — Figure 8.2.5. Network density map (McSwiggan Gaussian method)
🔗

Finally, as noted above, [Rakshit et al. (2019)] developed a fast algorithm, using a convolution method, that is particularly helpful with large datasets. The analysis of the authors suggest this estimator "is consistent, and its statistical efficiency is only slightly suboptimal" ([Rakshit et al., 2019]).

🔗

system.time(d_rakshit <- densityQuick.lpp(shoplifting_lpp,
                                    weights = marks(shoplifting_lpp),
                                    sigma = 100,
                                    edge2D = TRUE))

##    user  system elapsed 
##   0.250   0.002   0.254

We can now plot the data and compare the results:

🔗

plot(d_rakshit, main="Rakshit")

Once again we have a shaded street map of Heswall titled ’Rakshit’, with a similar legend, which ranges from zero to point two five this time. Most of Heswall is dark blue, with only the ends of two roads near the top right brightening to yellow. — Figure 8.2.6. Network density map (Rakshit method)
🔗

So these are some ways in which we can create a density map of crime along networks. You can compare the three maps to see how choice of method affects conclusions. This is something we must always reflect on, when applying these techniques.

🔗

Subsection 8.2.1 Adding context

Another thing you might like to do is add some context to explore where these hot spots are. Similar to how we did this with our KDE in continuous space in Chapter 7, we can use the raster and leaflet packages here too. One difference, we must also use the as.im() function to make sure our heatmap is of class "im".

🔗

# convert to raster wrapped inside as.im() function
r1 <- raster(as.im(d_mcswiggan))
#remove very low density values
r1[r1 < 0.0001 ] <- NA
# set CRS
crs(r1) <- "EPSG:27700"

#create a colour palette
pal <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"),
                    values(r1), na.color = "transparent")


#and then make map!
leaflet() %>%
  addTiles() %>%
  addRasterImage(r1, colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(r1),
    title = "Shoplifting map")

A leaflet map of Heswall and area, with a legend in the top right labeled 'Shoplifting map' matching a gradient from dark blue at zero to lighter blue at point zero seven. Several blotches of darker blue are spread throughout the city, bleeding into lighter blue a few times.

Figure 8.2.7. Add interactive basemap to explore hot spots

🔗

Prev Top Next