# Chapter 6 Point Pattern Analysis Using R

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

## 6.3 Techniques for Point Patterns Using R

### 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)

### 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'))

## 6.4 Using The $$K$$-function in R

# K-function code block
require(spatstat)
# Obtain the bramble cane data
data(bramblecanes)
plot(bramblecanes)
kf <- Kest(bramblecanes,correction='border')
# Plot it
plot(kf)
# Code block to produce k-function with envelope
# Envelope function
kf.env <- envelope(bramblecanes,Kest,correction="border")
# Plot it
plot(kf.env)
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)
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)

## 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)
clenv.bramble <- envelope(bramblecanes,Lcross,i=0,j=1,correction='border')
plot(clenv.bramble)
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) 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) # 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) ## 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) 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) predmat3 <- matrix(krig.est$var1.pred,length(ux),length(uy))
persp(predmat3,box=FALSE)

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

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_scale_bar(position=c("right","top"))
evgm <- variogram(fulmar~1,fulmar.spdf,
plot(evgm,model=fvgm)