Spatial Point Patterns: Methodology and Applications with R
by Adrian Baddeley, Ege Rubak and Rolf Turner


R code for chapter 5

Below is the R code used to generate results and figures in chapter 5. The code is in a rather raw format extracted from the book manuscript files – please read the instructions for use if you haven’t done so yet. You can download the script here.

### R code from vignette source '05pointproc.Rnw'
## Copyright (C) Adrian Baddeley, Ege Rubak and Rolf Turner

###################################################
### code chunk number 1: 05pointproc.Rnw:9-10
###################################################
source("R/startup.R")


###################################################
### code chunk number 2: 05pointproc.Rnw:76-79
###################################################
set.seed(19171026)
Xlist <- rMatClust(20, 0.2, 5, nsim=10)
nX <- length(Xlist)


###################################################
### code chunk number 3: Unit5x2.Rnw:4-5
###################################################
zeromargins()


###################################################
### code chunk number 4: 05pointproc.Rnw:99-102
###################################################
plot(Xlist, main="", main.panel="", equal.scales=TRUE, 
     pch=16, cex=0.75, 
     nrows=2, mar.panel=0, hsep=1, vsep=1)


###################################################
### code chunk number 5: potatoes
###################################################
load("data/potatoes.rda") # provides potato-shaped windows A and B 


###################################################
### code chunk number 6: matclex
###################################################
set.seed(13042024)
X <- rpoispp(50)
W <- as.owin(X)


###################################################
### code chunk number 7: Unit.Rnw:3-5
###################################################
newplot(6, 0.7)
setmargins(0)


###################################################
### code chunk number 8: 05pointproc.Rnw:205-209
###################################################
PCEX <- 1.0
plot(as.owin(X), main="", type="n")
plot(A, add=TRUE, lwd=2, col="grey")
plot(X, add=TRUE, pch=16, cex=PCEX)


###################################################
### code chunk number 9: 05pointproc.Rnw:273-275
###################################################
set.seed(1234560)
X1 <- runifpoint(1, A) ## cheating: ensures X falls in A.


###################################################
### code chunk number 10: Unit2.Rnw:3-5
###################################################
newplot(12.5, 0.9)
setmargins(0)


###################################################
### code chunk number 11: 05pointproc.Rnw:279-280
###################################################
setmargins(0)


###################################################
### code chunk number 12: 05pointproc.Rnw:284-294
###################################################
par(mfrow=c(1,2))
PCEX <- 1.0
plot(square(1), main="")
plot(X1, add=TRUE, pch=16, cex=PCEX)
segments(X1$x,    0, X1$x, 1,    lty=2, lwd=3)
segments(   0, X1$y,    1, X1$y, lty=2, lwd=3)
plot(square(1), main="")
plot(A, add=TRUE, lwd=2, col="grey")
plot(X1, add=TRUE, pch=16, cex=PCEX)
par(mfrow=c(1,1))


###################################################
### code chunk number 13: 05pointproc.Rnw:391-392
###################################################
Bin8list <- runifpoint(8, nsim=10)


###################################################
### code chunk number 14: Unit5x2.Rnw:4-5
###################################################
zeromargins()


###################################################
### code chunk number 15: 05pointproc.Rnw:397-400
###################################################
plot(Bin8list, main="", main.panel="", equal.scales=TRUE, 
     pch=16, cex=0.75, 
     nrows=2, mar.panel=0, hsep=1, vsep=1)


###################################################
### code chunk number 16: 05pointproc.Rnw:444-448
###################################################
set.seed(10982178)
PoisList <- rpoispp(50, nsim=10)
Pois8List <- rpoispp(8, nsim=10)
Pois8List[[7]] <- runifpoint(8) # cheat for purposes of illustration.


###################################################
### code chunk number 17: Unit5x2.Rnw:4-5
###################################################
zeromargins()


###################################################
### code chunk number 18: 05pointproc.Rnw:454-457
###################################################
plot(PoisList, main="", main.panel="", equal.scales=TRUE, 
     pch=16, cex=0.75, 
     nrows=2, mar.panel=0, hsep=1, vsep=1)


###################################################
### code chunk number 19: Unit5x2.Rnw:4-5
###################################################
zeromargins()


###################################################
### code chunk number 20: 05pointproc.Rnw:510-513
###################################################
plot(Pois8List, main="", main.panel="", equal.scales=TRUE, 
     pch=16, cex=0.75, 
     nrows=2, mar.panel=0, hsep=1, vsep=1)


###################################################
### code chunk number 21: 05pointproc.Rnw:541-558
###################################################
pb <- function(i, y, ..., main) {
  plot(as.owin(y), main=main, add=TRUE, show.all=TRUE, col=grey(0.4))
  R <- as.rectangle(y)
  vec <- unlist(centroid.owin(R))-0.5
  AA <- shift(A, vec)
  BB <- shift(B, vec)
  plot(AA, add=TRUE, col=grey(0.7), border=NA)
  plot(BB, add=TRUE, col=grey(0.7), border=NA)
}
pe <- function(i, y, ...) {
  R <- as.rectangle(y)
  vec <- unlist(centroid.owin(R))-0.5
  AA <- shift(A, vec)
  BB <- shift(B, vec)
  text(centroid.owin(AA), labels=npoints(y[AA]), cex=1.25)
  text(centroid.owin(BB), labels=npoints(y[BB]), cex=1.25)
}


###################################################
### code chunk number 22: 05pointproc.Rnw:561-563
###################################################
set.seed(26011788)
idsim <- solapply(1:3, function(z) { rpoispp(100) })


###################################################
### code chunk number 23: Unit3.Rnw:3-5
###################################################
newplot(19,0.9)
zeromargins()


###################################################
### code chunk number 24: 05pointproc.Rnw:569-573
###################################################
plot(idsim, panel.begin=pb, panel.end=pe, equal.scales=TRUE, 
     main="", main.panel="", 
     mar.panel=c(0, 0.2, 0.2, 0), hsep=1,
     pch=21, bg="white", cols="white")


###################################################
### code chunk number 25: Unit.Rnw:3-5
###################################################
newplot(6, 0.7)
setmargins(0)


###################################################
### code chunk number 26: 05pointproc.Rnw:604-607
###################################################
plot(W, col="grey", main="")
plot(X, add=TRUE, pch=16, cols="white", cex=1.5)
plot(quadratcount(X, 5), add=TRUE, cex=2.5)


###################################################
### code chunk number 27: 05pointproc.Rnw:609-612
###################################################
plot(W, col="grey", main="")
plot(X, add=TRUE, pch=16, cols="white", cex=1.5)
plot(quadratcount(X, 15), add=TRUE, cex=1.2)


###################################################
### code chunk number 28: 05pointproc.Rnw:656-667
###################################################
N <- 20
QN <- quadrats(X, N)
TN <- tiles(QN)
ZN <- quadratcount(X, N)
AN <- pixellate(rebound.owin(A, W), dimyx=N)
AN <- t(as.matrix(AN)[N:1,])
hit <- as.vector(AN > 0.2/N^2)
cen <- lapply(TN, centroid.owin)
x0hit <- unlist(lapply(cen, getElement, "x"))[hit]
y0hit <- unlist(lapply(cen, getElement, "y"))[hit]
ZNhit <- t(as.matrix(ZN))[t(matrix(hit, N, N, byrow=TRUE))]


###################################################
### code chunk number 29: Unit.Rnw:3-5
###################################################
newplot(6, 0.7)
setmargins(0)


###################################################
### code chunk number 30: 05pointproc.Rnw:674-678
###################################################
plot(W, col="grey", main="")
plot(X, add=TRUE, pch=16, cols="white", cex=1.1)
text(x0hit, y0hit, pmin(1, ZNhit))
plot(A, add=TRUE)


###################################################
### code chunk number 31: 05pointproc.Rnw:721-722
###################################################
setmargins(4,4,1,1)


###################################################
### code chunk number 32: 05pointproc.Rnw:726-734
###################################################
par(mfrow=c(1,2))
kk <- 0:6
d1 <- dpois(kk, 0.6)
d2 <- dpois(kk, 1.7)
yr <- range(d1, d2)
barplot(d1, names.arg=kk, ylab="Probability", ylim=yr, xlab="Count")
barplot(d2, names.arg=kk, ylab="Probability", ylim=yr, xlab="Count")
par(mfrow=c(1,1))


###################################################
### code chunk number 33: 05pointproc.Rnw:801-809
###################################################
set.seed(1848)
X <- rpoispp(100)
retain <- (runif(npoints(X)) < 0.4)
Y <- X[retain]
ThinDemo <- solist(layered(X, plotargs=list(pch=16)),
                   layered(X[retain], X[!retain],
                           plotargs=list(list(pch=16), list(pch=1))),
                   layered(Y, plotargs=list(pch=16)))


###################################################
### code chunk number 34: Unit3.Rnw:3-5
###################################################
newplot(19,0.9)
zeromargins()


###################################################
### code chunk number 35: 05pointproc.Rnw:814-816
###################################################
plot(ThinDemo, main="", main.panel="", nrows=1, 
     equal.scales=TRUE, mar.panel=0, hsep=1)


###################################################
### code chunk number 36: 05pointproc.Rnw:838-842
###################################################
set.seed(1967)
X <- rpoispp(40)
Y <- rpoispp(60)
Z <- superimpose(X,Y, W=square(1))


###################################################
### code chunk number 37: Unit3.Rnw:3-5
###################################################
newplot(19,0.9)
zeromargins()


###################################################
### code chunk number 38: 05pointproc.Rnw:847-850
###################################################
plot(solist(X, Y, Z), 
     main="", main.panel="", nrows=1, 
     equal.scales=TRUE, mar.panel=0, hsep=1)


###################################################
### code chunk number 39: 05pointproc.Rnw:879-880 (eval = FALSE)
###################################################
## plot(rpoispp(50))


###################################################
### code chunk number 40: 05pointproc.Rnw:889-890 (eval = FALSE)
###################################################
## plot(rpoispp(100, nsim=9))


###################################################
### code chunk number 41: 05pointproc.Rnw:910-911 (eval = FALSE)
###################################################
## plot(rpoispp(100, win=letterR))


###################################################
### code chunk number 42: Rletter.Rnw:3-5
###################################################
newplot(4, 0.3)
setmargins(0)


###################################################
### code chunk number 43: 05pointproc.Rnw:917-918
###################################################
plot(rpoispp(100, win=letterR), main="", pch=16)


###################################################
### code chunk number 44: 05pointproc.Rnw:1056-1058
###################################################
lambda <- function(x,y) { 100 * (x^2+y) }
X <- rpoispp(lambda, win=square(1))


###################################################
### code chunk number 45: 05pointproc.Rnw:1060-1064
###################################################
lamim <- as.im(lambda, square(1))
L <- layered(lamim, lamim,
             plotargs=list(list(),
                           list(.plot="contour", col="white")))


###################################################
### code chunk number 46: Unit2r.Rnw:3-5
###################################################
newplot(12.5, 0.9)
setmargins(0,0,0,1)


###################################################
### code chunk number 47: 05pointproc.Rnw:1071-1074
###################################################
plot(solist(L, X),
     main="", main.panel="", equal.scales=TRUE, 
     pch=16, hsep=1, mar.panel=0.2)


###################################################
### code chunk number 48: 05pointproc.Rnw:1119-1121
###################################################
set.seed(3332)
Rad <- 0.15


###################################################
### code chunk number 49: 05pointproc.Rnw:1123-1124
###################################################
Xmc <- rMatClust(kappa=8, r=0.15, mu=5)


###################################################
### code chunk number 50: 05pointproc.Rnw:1126-1129
###################################################
Ymc <- attr(Xmc, "parents")
WB <- as.owin(Ymc)
W <- as.owin(Xmc)


###################################################
### code chunk number 51: Unit3.Rnw:3-5
###################################################
newplot(19,0.9)
zeromargins()


###################################################
### code chunk number 52: 05pointproc.Rnw:1134-1135
###################################################
setmargins(0.5)


###################################################
### code chunk number 53: 05pointproc.Rnw:1139-1154
###################################################
CEX <- 1.2
opa <- par(mfrow=c(1,3), mar=rep(0.5,4))
plot(WB, type="n", main="")
plot(W, add=TRUE)
plot(Ymc, add=TRUE, pch=16, cex=CEX)

plot(WB, type="n", main="")
plot(W, add=TRUE)
plot(Ymc, pch=16,add=TRUE, cex=CEX)
for(i in 1:npoints(Ymc)) plot(disc(Rad, Ymc[i]), add=TRUE, lty=3)
plot(Xmc, add=TRUE, cex=CEX)

plot(WB, type="n", main="")
plot(W, add=TRUE)
plot(Xmc, add=TRUE, cex=CEX)


###################################################
### code chunk number 54: 05pointproc.Rnw:1189-1190
###################################################
XmI <- rMaternI(30, 0.1)


###################################################
### code chunk number 55: 05pointproc.Rnw:1192-1195
###################################################
R <- 0.1
P <- rpoispp(30)
kill <- (nndist(P) < R)


###################################################
### code chunk number 56: Unit3.Rnw:3-5
###################################################
newplot(19,0.9)
zeromargins()


###################################################
### code chunk number 57: 05pointproc.Rnw:1201-1210
###################################################
par(mfrow=c(1,3))
# left
plot(P, main="")
# middle
plot(P[kill], pch=4, main="")
points(P[!kill], pch=16)
# right
plot(P[!kill], main="", pch=16)
par(mfrow=c(1,1))


###################################################
### code chunk number 58: 05pointproc.Rnw:1231-1232
###################################################
XmII <- rMaternII(30, 0.1)


###################################################
### code chunk number 59: 05pointproc.Rnw:1234-1239
###################################################
age <- runif(npoints(P))
D <- pairdist(P)
close <- (D < R)
earlier <- outer(age, age, "<")
killyoung <- apply(close & earlier, 1, any)


###################################################
### code chunk number 60: Unit3.Rnw:3-5
###################################################
newplot(19,0.9)
zeromargins()


###################################################
### code chunk number 61: 05pointproc.Rnw:1245-1254
###################################################
par(mfrow=c(1,3))
# left
plot(P, main="")
# middle
plot(P[killyoung], pch=4, main="")
points(P[!killyoung], pch=16)
# right
plot(P[!killyoung], main="", pch=16)
par(mfrow=c(1,1))


###################################################
### code chunk number 62: 05pointproc.Rnw:1288-1289
###################################################
Xs <- rSSI(0.05, n = 200)


###################################################
### code chunk number 63: Unit.Rnw:3-5
###################################################
newplot(6, 0.7)
setmargins(0)


###################################################
### code chunk number 64: 05pointproc.Rnw:1296-1297
###################################################
plot(Xs, pch=16, main="")


###################################################
### code chunk number 65: 05pointproc.Rnw:1481-1484
###################################################
W <- square(1)
R2 <- grow.rectangle(W, 2, 0.5)
X <- rpoispp(50,win=R2)


###################################################
### code chunk number 66: 05pointproc.Rnw:1490-1494
###################################################
plot(R2, type="n", main="")
plot(X[W], pch=16,add=TRUE)
plot(X[complement.owin(W, R2)], add=TRUE)
plot(W, add=TRUE)


###################################################
### code chunk number 67: Unit2.Rnw:3-5
###################################################
newplot(12.5, 0.9)
setmargins(0)


###################################################
### code chunk number 68: 05pointproc.Rnw:1530-1541
###################################################
set.seed(7891095)
ncell <- 16
X <- rSSI(sqrt(2)/ncell,12)
XX <- layered(X, plotargs=list(pch=16))
Y <- quadratcount(X, ncell)
is0 <- (as.vector(t(as.table(Y))) == 0)
cols <- ifelse(is0, "grey", "black")
YY <- layered(as.owin(Y), Y, 
              plotargs=list(list(), list(show.tiles=FALSE, col=cols)))
plot(solist(XX, YY), main="", main.panel="", equal.scales=TRUE,
     mar.panel=rep(0.2,4), hsep=1)


###################################################
### code chunk number 69: 05pointproc.Rnw:1587-1588
###################################################
GR <- rotate(gordon, 0.95)


###################################################
### code chunk number 70: GordonRot.Rnw:3-5
###################################################
newplot(6.6, 0.85)
setmargins(0)


###################################################
### code chunk number 71: 05pointproc.Rnw:1594-1596
###################################################
plot(Window(GR), col="grey", main="")
plot(GR, add=TRUE, pch=16)


###################################################
### code chunk number 72: 05pointproc.Rnw:1825-1826
###################################################
setmargins(0.1)


###################################################
### code chunk number 73: 05pointproc.Rnw:1830-1836
###################################################
Y <- bei.extra$elev
par(mfrow=c(1,2))
plot(Y, main="", axes=FALSE)
contour(Y, add=TRUE, nlevels=6)
persp(Y, main="", theta=20, phi=30, border=NA, shade=0.6, zlab="Elevation")
par(mfrow=c(1,1))


###################################################
### code chunk number 74: 05pointproc.Rnw:1881-1887
###################################################
Y <- rotate(copper[["SouthLines"]], pi/2)
Z <- distmap(Y, eps=0.25)
ZY <- layered(Z, Y, Window(Y), 
              plotargs=list(list(ribbon=FALSE, col=lighttodark),
                            list(lwd=2, col="white"),
                            list()))


###################################################
### code chunk number 75: CopperS2.Rnw:3-5
###################################################
newplot(6, 0.95)
setmargins(0)


###################################################
### code chunk number 76: 05pointproc.Rnw:1893-1896
###################################################
plot(solist(Y, ZY), ncols=1,
     main="", main.panel="", equal.scales=TRUE, halign=TRUE,
     mar.panel=0, vsep=0.5)


###################################################
### code chunk number 77: 05pointproc.Rnw:2021-2027
###################################################
swp <- rescale(swedishpines)
lam <- predict(ppm(swp ~ polynom(x,y,2)))
K <- Kest(swp)
Ki <- Kinhom(swp, lam)
coco <- layered(swp, lam,
                plotargs=list(list(cols="grey"), list(.plot="contour")))


###################################################
### code chunk number 78: Unit2.Rnw:3-5
###################################################
newplot(12.5, 0.9)
setmargins(0)


###################################################
### code chunk number 79: 05pointproc.Rnw:2033-2035
###################################################
plot(solist(swp, coco), pch=16, main="", main.panel="", 
     mar.panel=0, hsep=1)


###################################################
### code chunk number 80: fv2.Rnw:3-5
###################################################
newplot(12, 0.95)
setmargins(0.5+c(3,3,0,1))


###################################################
### code chunk number 81: 05pointproc.Rnw:2047-2048
###################################################
plot(anylist(K, Ki), main="", main.panel="", legend=FALSE, xlim=c(0, 2))