Chapter 6 Point Pattern Analysis Using R

library(tidyverse)
library(GISTools)
library(sp)
library(rgeos)
library(tmap)
library(tmaptools)

6.1 Introduction

6.2 What is special about spatial?

6.2.1 Point Patterns

6.3 Techniques for Point Patterns Using R

6.3.1 Kernel Density Estimates

6.3.2 Kernel Density Estimation Using R

# Load GISTools (for the data) and tmap (for the mapping)
require(GISTools)
require(tmap)

# Get the data
data(newhaven)
# look at it
# select 'view' mode
tmap_mode('view')
# Create the map of blocks and incidents
tm_shape(blocks) + tm_borders() + tm_shape(breach) +
  tm_dots(col='navyblue')

# Function to choose bandwidth according Bowman and Azzalini / Scott's rule
# for use with <smooth_map> in <tmaptools>

choose_bw <- function(spdf) {
  X <- coordinates(spdf)
  sigma <- c(sd(X[,1]),sd(X[,2]))  * (2 / (3 * nrow(X))) ^ (1/6)
  return(sigma/1000)
}

library(tmaptools)
tmap_mode('view')
breach_dens <- smooth_map(breach,cover=blocks, bandwidth = choose_bw(breach))

  |                                                        
  |                                                  |   0%
  |                                                        
  |=====                                             |  10%
  |                                                        
  |===============                                   |  30%
  |                                                        
  |=========================                         |  50%
  |                                                        
  |===================================               |  70%
  |                                                        
  |=============================================     |  90%
  |                                                        
  |==================================================| 100%

tm_shape(breach_dens$raster) + tm_raster()

tmap_mode('view')
tm_shape(blocks)+ tm_borders(alpha=0.5) +
  tm_shape(breach_dens$iso) + tm_lines(col='darkred',lwd=2) 

6.3.3 Further Uses of Kernel Density Estimation

# R Kernel Density comparison - first make sure the New Haven data is available
require(GISTools)
data(newhaven)

tmap_mode('plot')
# Create the KDEs for the two data sets:
contours <- seq(0,1.4,by=0.2)

brn_dens <- smooth_map(burgres.n,cover=blocks, breaks=contours, 
                       style='fixed',
                       bandwidth = choose_bw(burgres.n))

  |                                                        
  |                                                  |   0%
  |                                                        
  |=====                                             |  10%
  |                                                        
  |===============                                   |  30%
  |                                                        
  |=========================                         |  50%
  |                                                        
  |===================================               |  70%
  |                                                        
  |=============================================     |  90%
  |                                                        
  |==================================================| 100%

brf_dens <- smooth_map(burgres.f,cover=blocks, breaks=contours, 
                       style='fixed',
                       bandwidth = choose_bw(burgres.f))

  |                                                        
  |                                                  |   0%
  |                                                        
  |=====                                             |  10%
  |                                                        
  |===============                                   |  30%
  |                                                        
  |=========================                         |  50%
  |                                                        
  |===================================               |  70%
  |                                                        
  |=============================================     |  90%
  |                                                        
  |==================================================| 100%

# Create the maps and store them in variables
dn <- tm_shape(blocks) + tm_borders() + 
  tm_shape(brn_dens$polygons) + tm_fill(col='level',alpha=0.8) +
  tm_layout(title="Non-Forced Burglaries")
df <- tm_shape(blocks) + tm_borders() + 
  tm_shape(brf_dens$polygons) + tm_fill(col='level',alpha=0.8) + 
  tm_layout(title="Forced Burglaries")

tmap_arrange(dn,df)

Figure 6.1: KDE Maps to Compare Forced and Non-Forced Burglary Patterns

6.3.4 Hexagonal Binning Using R

install.packages("fMultivar",depend=TRUE)

hexbin_map <- function(spdf, ...) {
  hbins <- fMultivar::hexBinning(coordinates(spdf),...)

  # Hex binning code block
  # Set up the hexagons to plot,  as polygons
  u <- c(1, 0, -1, -1, 0, 1)
  u <- u * min(diff(unique(sort(hbins$x))))
  v <- c(1,2,1,-1,-2,-1)
  v <- v * min(diff(unique(sort(hbins$y))))/3

  # Construct each polygon in the sp model 
  hexes_list <- vector(length(hbins$x),mode='list')
  for (i in 1:length(hbins$x)) {
   pol <- Polygon(cbind(u + hbins$x[i], v + hbins$y[i]),hole=FALSE)
   hexes_list[[i]] <- Polygons(list(pol),i) }
  
  # Build the spatial polygons data frame
  hex_cover_sp <- SpatialPolygons(hexes_list,proj4string=CRS(proj4string(spdf)))
  hex_cover <- SpatialPolygonsDataFrame(hex_cover_sp,
                                      data.frame(z=hbins$z),match.ID=FALSE)
  # Return the result
  return(hex_cover)
}

tmap_mode('view')
breach_hex <- hexbin_map(breach,bins=20)
tm_shape(breach_hex) + 
  tm_fill(col='z',title='Count',alpha=0.7)

hexprop_map <- function(spdf, ...) {
  hbins <- fMultivar::hexBinning(coordinates(spdf),...)

  # Hex binning code block
  # Set up the hexagons to plot,  as polygons
  u <- c(1, 0, -1, -1, 0, 1)
  u <- u * min(diff(unique(sort(hbins$x))))
  v <- c(1,2,1,-1,-2,-1)
  v <- v * min(diff(unique(sort(hbins$y))))/3

  scaler <- sqrt(hbins$z/max(hbins$z))
  # Construct each polygon in the sp model 
  hexes_list <- vector(length(hbins$x),mode='list')
  for (i in 1:length(hbins$x)) {
   pol <- Polygon(cbind(u*scaler[i] + hbins$x[i], v*scaler[i] + hbins$y[i]),hole=FALSE)
   hexes_list[[i]] <- Polygons(list(pol),i) }
  
  # Build the spatial polygons data frame
  hex_cover_sp <- SpatialPolygons(hexes_list,proj4string=CRS(proj4string(spdf)))
  hex_cover <- SpatialPolygonsDataFrame(hex_cover_sp,
                                      data.frame(z=hbins$z),match.ID=FALSE)
  # Return the result
  return(hex_cover)
}

tmap_mode('plot')
breach_prop <- hexprop_map(breach,bins=20)
tm_shape(blocks) + tm_borders(col='grey') +
  tm_shape(breach_prop) + 
  tm_fill(col='indianred',alpha=0.7) + 
  tm_layout("Breach of Peace Incidents",title.position=c('right','top'))

Figure 6.2: Hexagonal Binning of residential burglaries

6.3.5 Second Order Analysis of Point Patterns

6.4 Using The \(K\)-function in R

# K-function code block
# Load the spatstat package
require(spatstat)
# Obtain the bramble cane data
data(bramblecanes)
plot(bramblecanes)

Figure 6.3: Bramble cane locations

kf <- Kest(bramblecanes,correction='border')
# Plot it
plot(kf)

Figure 6.4: Ripley’s \(k\) function plot

# Code block to produce k-function with envelope
# Envelope function
kf.env <- envelope(bramblecanes,Kest,correction="border")
# Plot it
plot(kf.env)

Figure 6.5: \(k\) function with envelope

mad.test(bramblecanes,Kest,verbose=FALSE)

    Maximum absolute deviation test of CSR
    Monte Carlo test based on 99 simulations
    Summary function: K(r)
    Reference function: theoretical
    Alternative: two.sided
    Interval of distance values: [0, 0.25] units (one
    unit = 9 metres)
    Test statistic: Maximum absolute deviation
    Deviation = observed minus theoretical

data:  bramblecanes
mad = 0.016159, rank = 1, p-value = 0.01

dclf.test(bramblecanes,Kest,verbose=FALSE)

    Diggle-Cressie-Loosmore-Ford test of CSR
    Monte Carlo test based on 99 simulations
    Summary function: K(r)
    Reference function: theoretical
    Alternative: two.sided
    Interval of distance values: [0, 0.25] units (one
    unit = 9 metres)
    Test statistic: Integral of squared absolute
    deviation
    Deviation = observed minus theoretical

data:  bramblecanes
u = 3.3372e-05, rank = 1, p-value = 0.01

6.4.1 The \(L\)-function

# Code block to produce k-function with envelope
# Envelope function
lf.env <- envelope(bramblecanes,Lest,correction="border")
# Plot it
plot(lf.env)

Figure 6.6: \(L\) function with envelope

mad.test(bramblecanes,Lest,verbose=FALSE)

    Maximum absolute deviation test of CSR
    Monte Carlo test based on 99 simulations
    Summary function: L(r)
    Reference function: theoretical
    Alternative: two.sided
    Interval of distance values: [0, 0.25] units (one
    unit = 9 metres)
    Test statistic: Maximum absolute deviation
    Deviation = observed minus theoretical

data:  bramblecanes
mad = 0.017759, rank = 1, p-value = 0.01

6.5 The \(G\)-function

# Code block to produce G-function with envelope
# Envelope function
gf.env <- envelope(bramblecanes,Gest,correction="border")
# Plot it
plot(gf.env)

Figure 6.7: \(G\) function with envelope

6.6 Looking at Marked Point Patterns

6.7 Cross-\(L\)-function analysis in R

marks(bramblecanes)

  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [28] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [55] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [82] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[109] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[136] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[163] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[190] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[217] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[244] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[271] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[298] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[325] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[352] 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[379] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[406] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[433] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[460] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[487] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[514] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[541] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[568] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[595] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[622] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[649] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[676] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[703] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[730] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
[757] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[784] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[811] 2 2 2 2 2 2 2 2 2 2 2 2 2
Levels: 0 1 2

cl.bramble <- Lcross(bramblecanes,i=0,j=1,correction='border')
plot(cl.bramble)

Figure 6.8: Cross-\(L\) function for levels 0 and 1 of the bramble cane data

clenv.bramble <- envelope(bramblecanes,Lcross,i=0,j=1,correction='border')
plot(clenv.bramble)

Figure 6.9: Cross-\(L\) function envelope for levels 0 and 1 of the bramble cane data

dclf.test(bramblecanes,Lcross,i=0,j=1,correction='border',verbose=FALSE)

    Diggle-Cressie-Loosmore-Ford test of CSR
    Monte Carlo test based on 99 simulations
    Summary function: L["0", "1"](r)
    Reference function: theoretical
    Alternative: two.sided
    Interval of distance values: [0, 0.25] units (one
    unit = 9 metres)
    Test statistic: Integral of squared absolute
    deviation
    Deviation = observed minus theoretical

data:  bramblecanes
u = 4.3982e-05, rank = 1, p-value = 0.01

6.8 Interpolation of Point Patterns With Continuous Attributes

6.8.1 Nearest Neighbour Interpolation

#
# Original code from Carson Farmer
# http://www.carsonfarmer.com/2009/09/voronoi-polygons-with-r/
# Subject to minor stylistic modifications
#

require(deldir)
require(sp)

# Modified Carson Farmer code
voronoipolygons = function(layer) {
  crds <- layer@coords
  z <- deldir(crds[,1], crds[,2])
  w <- tile.list(z)
  polys <- vector(mode='list', length=length(w))
  for (i in seq(along=polys)) {
    pcrds <- cbind(w[[i]]$x, w[[i]]$y)
    pcrds <- rbind(pcrds, pcrds[1,])
    polys[[i]] <- Polygons(list(Polygon(pcrds)),
                           ID=as.character(i))
  }
  SP <- SpatialPolygons(polys)
  voronoi <- SpatialPolygonsDataFrame(SP, 
                data=data.frame(x=crds[,1], 
                                y=crds[,2], 
                                layer@data,
                                row.names=sapply(slot(SP, 'polygons'), 
                                                 function(x) slot(x, 'ID'))))
  proj4string(voronoi) <- CRS(proj4string(layer))
  return(voronoi)
}

6.8.2 A Look at the Data

library(gstat)
library(tmap)
data(fulmar)

fulmar.spdf <- SpatialPointsDataFrame(cbind(fulmar$x,fulmar$y),
                                      fulmar)
fulmar.spdf <- fulmar.spdf[fulmar.spdf$year==1999,]
proj4string(fulmar.spdf) <- CRS("+init=epsg:32631")
fulmar.voro <- voronoipolygons(fulmar.spdf)

tmap_mode('plot')
fpt <- tm_shape(fulmar.spdf) + tm_dots(size=0.1)
fvr <- tm_shape(fulmar.voro) + tm_borders()
tmap_arrange(fpt,fvr)

Figure 6.10: Fulmar Sighting Transects (LHS=Points; RHS=Voronoi Diagram)

library(gstat)
library(GISTools)
tmap_mode('view')
sh <- shading(breaks=c(5,15,25,35),
              cols=brewer.pal(5,'Purples'))
tm_shape(fulmar.voro) + tm_fill(col='fulmar',style='fixed',breaks=c(0,5,15,25,35),
                                alpha=0.6,title='Fulmar Density')

6.8.3 Inverse Distance Weighting (IDW)

6.8.4 Computing IDW with the gstat Package

library(maptools) # Required package
library(GISTools) # Required package
library(gstat) # Set up the gstat package
# Define a sample grid then use it as a set of points 
# to estimate fulmar density via IDW, with alpha=1
s.grid <- spsample(fulmar.voro,type='regular',n=6000)
idw.est <- gstat::idw(fulmar~1,fulmar.spdf,
                      newdata=s.grid,idp=1.0)

tmap_mode('view')
tm_shape(idw.est) + tm_dots(col='var1.pred',border.col=NA,alpha=0.7)

tmap_mode('view')
idw.grid <- SpatialPixelsDataFrame(idw.est,data.frame(idw.est)) 
tm_shape(idw.grid) + tm_raster(col='var1.pred',title='Fulmar')

require(raster)
levs <- c(0,2,4,6,8,Inf)
idw.raster <- raster(idw.grid,layer='var1.pred') 
idw.contour <- rasterToContour(idw.raster,levels=levs) 

tmap_mode('view')
tm_shape(idw.contour) + tm_iso()

tmap_mode('view')
idw.levels <- tmaptools:::lines2polygons(fulmar.voro,
                idw.contour,rst=idw.raster,lvls=levs) 
tm_shape(idw.levels) + tm_fill(col='level') 

idw.est2 <- gstat::idw(fulmar~1,fulmar.spdf,
                      newdata=s.grid,idp=2.0)
idw.grid2 <- SpatialPixelsDataFrame(idw.est2,data.frame(idw.est2)) 
tmap_mode('view')
tm_shape(idw.grid2) +  tm_raster(col='var1.pred',title='Fulmar',breaks=levs)

tmap_mode('plot')
idw1 <- tm_shape(idw.grid) +  tm_raster(col='var1.pred',title='Alpha = 1',breaks=levs)
idw2 <- tm_shape(idw.grid2) +  tm_raster(col='var1.pred',title='Alpha = 2',breaks=levs)
tmap_arrange(idw1,idw2)

Figure 6.11: Inverse Distance Weighting Estimate of Fulmar Density (Comparison)

# Extract the distinct x and y coordinates of the grid
# Extract the predicted values and form into a matrix
# of gridded values
ux <- unique(coordinates(idw.est)[,1])
uy <- unique(coordinates(idw.est)[,2])
predmat <- matrix(idw.est$var1.pred,length(ux),length(uy))

Next, the same is done for idw.est2

predmat2 <- matrix(idw.est2$var1.pred,length(ux),length(uy))

par(mfrow=c(1,2),mar=c(0,0,2,0))
persp(predmat,box=FALSE)
persp(predmat2,box=FALSE)

Figure 6.12: 3D Plots of IDW (LHS:\(\alpha=1\); RHS:\(\alpha=2\))

6.9 The Kriging approach

6.9.1 A Brief Introduction to Kriging

6.9.2 Random Functions

6.9.3 Estimating the Semivariogram

require(gstat)
evgm <- variogram(fulmar~1,fulmar.spdf,
                  boundaries=seq(0,250000,l=51))
fvgm <- fit.variogram(evgm,vgm(3,"Mat",100000,1))
plot(evgm,model=fvgm)

Figure 6.13: Kriging semivariogram

krig.est <- krige(fulmar~1,fulmar.spdf,newdata=s.grid,model=fvgm)
krig.grid <- SpatialPixelsDataFrame(krig.est,krig.est@data)
krig.map.est <- tm_shape(krig.grid) + 
  tm_raster(col='var1.pred',breaks=levs,title='Fulmar Density',palette='Reds') + 
  tm_layout(legend.bg.color='white',legend.frame = TRUE)
var.levs <- c(0,3,6,9,12,Inf) 
krig.map.var <- tm_shape(krig.grid) + 
  tm_raster(col='var1.var',breaks=var.levs,title='Estimate Variance',palette='Reds') + 
  tm_layout(legend.bg.color='white',legend.frame = TRUE)
tmap_arrange(krig.map.est,krig.map.var)

Figure 6.14: Kriging estimates of fulmar density (RHS), and associated variance (LHS)

predmat3 <- matrix(krig.est$var1.pred,length(ux),length(uy))
persp(predmat3,box=FALSE)

Figure 6.15: 3D Plot Kriging-based interpolation.

Self Test Question 2

Try fitting an exponential variogram to the Fulmar data, and creating the surface plot, and maps. You may want to look at the help for fit.variogram to find out how to specify alternative variogram models.

6.10 Concluding Remarks

Answer to Self-Test Question 1

tmap_mode('plot')
tm_shape(breach_dens$polygons) + 
  tm_fill(title='Intensity \n(Incidents per square Km.)',col='level') + 
  tm_shape(blocks) + tm_borders() + 
  tm_shape(roads) + tm_lines(lwd=0.5,col='brown') +
  tm_scale_bar(position=c("right","top"))

Figure 6.16: KDE Map for Breaches of The Peace, for Self-Test Question 1

Answer to Self-Test Question 2

evgm <- variogram(fulmar~1,fulmar.spdf,
                  boundaries=seq(0,250000,l=51))
fvgm <- fit.variogram(evgm,vgm(3,"Exp",100000,1))
plot(evgm,model=fvgm)

Figure 6.17: Kriging semivariogram (Exponential Model)