## ---- echo=FALSE, include=FALSE-----------------------------------------------
library(knitr)
library(rspatial)
opts_chunk$set(fig.width = 5, fig.height = 5, fig.cap = '', collapse  = TRUE)


## ----getData, echo=TRUE-------------------------------------------------------
if (!require("rspatial")) remotes::install_github('rspatial/rspatial')


## ---- message=FALSE-----------------------------------------------------------
library(rspatial)
city <- sp_data("city")
crime <- sp_data("crime.rds")


## ---- pp1, fig.width=8, fig.height=4------------------------------------------
par(mai=c(0,0,0,0))
plot(city, col='light blue')
points(crime, col='red', cex=.5, pch='+')


## -----------------------------------------------------------------------------
tb <- sort(table(crime$CATEGORY))[-1]
tb


## -----------------------------------------------------------------------------
xy <- coordinates(crime)
dim(xy)
xy <- unique(xy)
dim(xy)
head(xy)


## -----------------------------------------------------------------------------
# mean center
mc <- apply(xy, 2, mean)
# standard distance
sd <- sqrt(sum((xy[,1] - mc[1])^2 + (xy[,2] - mc[2])^2) / nrow(xy))


## ---- pp2, fig.width=8, fig.height=4.5----------------------------------------
plot(city, col='light blue')
points(crime, cex=.5)
points(cbind(mc[1], mc[2]), pch='*', col='red', cex=5)

# make a circle
bearing <- 1:360 * pi/180
cx <- mc[1] + sd * cos(bearing)
cy <- mc[2] + sd * sin(bearing)
circle <- cbind(cx, cy)
lines(circle, col='red', lwd=2)


## ---- message=FALSE-----------------------------------------------------------
CityArea <- area(city)
dens <- nrow(xy) / CityArea


## -----------------------------------------------------------------------------
r <- raster(city) 
res(r) <- 1000
r


## ---- pp3, fig.width=8, fig.height=4------------------------------------------
r <- rasterize(city, r)
plot(r)
quads <- as(r, 'SpatialPolygons')
plot(quads, add=TRUE)
points(crime, col='red', cex=.5)


## ---- pp4, fig.width=8, fig.height=4------------------------------------------
nc <- rasterize(coordinates(crime), r, fun='count', background=0)
plot(nc)
plot(city, add=TRUE)


## ---- pp5, fig.width=8, fig.height=4------------------------------------------
ncrimes <- mask(nc, r)
plot(ncrimes)
plot(city, add=TRUE)


## ---- pp6---------------------------------------------------------------------
f <- freq(ncrimes, useNA='no')
head(f)
plot(f, pch=20)


## -----------------------------------------------------------------------------
# number of quadrats
quadrats <- sum(f[,2])
# number of cases
cases <- sum(f[,1] * f[,2])
mu <- cases / quadrats
mu


## -----------------------------------------------------------------------------
ff <- data.frame(f)
colnames(ff) <- c('K', 'X')
ff$Kmu <- ff$K - mu
ff$Kmu2 <- ff$Kmu^2
ff$XKmu2 <- ff$Kmu2 * ff$X
head(ff)


## -----------------------------------------------------------------------------
s2 <- sum(ff$XKmu2) / (sum(ff$X)-1)
s2


## -----------------------------------------------------------------------------
VMR <- s2 / mu 
VMR


## -----------------------------------------------------------------------------
d <- dist(xy)
class(d)


## -----------------------------------------------------------------------------
dm <- as.matrix(d)
dm[1:5, 1:5]
diag(dm) <- NA
dm[1:5, 1:5]


## -----------------------------------------------------------------------------
dmin <- apply(dm, 1, min, na.rm=TRUE)
head(dmin)


## -----------------------------------------------------------------------------
mdmin <- mean(dmin)


## -----------------------------------------------------------------------------
wdmin <- apply(dm, 1, which.min)


## ---- pp10, fig.width=8, fig.height=4.5---------------------------------------
plot(city)
points(crime, cex=.1)
ord <- rev(order(dmin))

far25 <- ord[1:25]
neighbors <- wdmin[far25]

points(xy[far25, ], col='blue', pch=20)
points(xy[neighbors, ], col='red')

# drawing the lines, easiest via a loop
for (i in far25) {
	lines(rbind(xy[i, ], xy[wdmin[i], ]), col='red')
}


## ---- pp11--------------------------------------------------------------------
max(dmin)
# get the unique distances (for the x-axis)
distance <- sort(unique(round(dmin)))
# compute how many cases there with distances smaller that each x
Gd <- sapply(distance, function(x) sum(dmin < x)) 
# normalize to get values between 0 and 1
Gd <- Gd / length(dmin)
plot(distance, Gd)
# using xlim to exclude the extremes
plot(distance, Gd, xlim=c(0,500))


## -----------------------------------------------------------------------------
stepplot <- function(x, y, type='l', add=FALSE, ...) {
	x <- as.vector(t(cbind(x, c(x[-1], x[length(x)]))))
	y <- as.vector(t(cbind(y, y)))
  if (add) {
     lines(x,y, ...)    
  } else {
	   plot(x,y, type=type, ...)
  }
}


## ---- pp12--------------------------------------------------------------------
stepplot(distance, Gd, type='l', lwd=2, xlim=c(0,500))


## ---- pp13--------------------------------------------------------------------
# get the centers of the 'quadrats' (raster cells)
p <- rasterToPoints(r)
# compute distance from all crime sites to these cell centers
d2 <- pointDistance(p[,1:2], xy, longlat=FALSE)

# the remainder is similar to the G function
Fdistance <- sort(unique(round(d2)))
mind <- apply(d2, 1, min)
Fd <- sapply(Fdistance, function(x) sum(mind < x)) 
Fd <- Fd / length(mind)
plot(Fdistance, Fd, type='l', lwd=2, xlim=c(0,3000))


## -----------------------------------------------------------------------------
ef <- function(d, lambda) {
  E <- 1 - exp(-1 * lambda * pi * d^2)
}
expected <- ef(0:2000, dens)


## ---- pp14--------------------------------------------------------------------
plot(distance, Gd, type='l', lwd=2, col='red', las=1, 
    ylab='F(d) or G(d)', xlab='Distance', yaxs="i", xaxs="i")
lines(Fdistance, Fd, lwd=2, col='blue')
lines(0:2000, expected, lwd=2)

legend(1200, .3, 
   c(expression(italic("G")["d"]), expression(italic("F")["d"]), 'expected'),
   lty=1, col=c('red', 'blue', 'black'), lwd=2, bty="n")


## ---- pp15--------------------------------------------------------------------
distance <- seq(1, 30000, 100)
Kd <- sapply(distance, function(x) sum(d < x)) # takes a while
Kd <- Kd / (length(Kd) * dens)
plot(distance, Kd, type='l', lwd=2)


## ---- message=FALSE-----------------------------------------------------------
library(spatstat)


## ---- message=FALSE-----------------------------------------------------------
library(maptools)
cityOwin <- as.owin(city)
class(cityOwin)
cityOwin


## -----------------------------------------------------------------------------
pts <- coordinates(crime)
head(pts)


## ---- pp20, fig.width=8, fig.height=4-----------------------------------------
p <- ppp(pts[,1], pts[,2], window=cityOwin)
class(p)
p
par(mai=c(0,0,0,0))
plot(p)


## ---- pp21, fig.width=8, fig.height=4-----------------------------------------
ds <- density(p)
class(ds)
par(mai=c(0,0,0.5,0.5))
plot(ds, main='crime density')


## -----------------------------------------------------------------------------
nrow(pts)
r <- raster(ds)
s <- sum(values(r), na.rm=TRUE) 
s * prod(res(r))


## -----------------------------------------------------------------------------
str(ds)
sum(ds$v, na.rm=TRUE) * ds$xstep * ds$ystep
p$n


## -----------------------------------------------------------------------------
census <- sp_data("census2000.rds")


## -----------------------------------------------------------------------------
census$area <- area(census)
census$area <- census$area/27878400
census$dens <- census$POP2000 / census$area


## -----------------------------------------------------------------------------
p <- coordinates(census)
head(p)


## -----------------------------------------------------------------------------
win <- aggregate(census)


## ---- pp22--------------------------------------------------------------------
plot(census)
points(p, col='red', pch=20, cex=.25)
plot(win, add=TRUE, border='blue', lwd=3)


## -----------------------------------------------------------------------------
owin <- as.owin(win)
pp <- ppp(p[,1], p[,2], window=owin, marks=census$dens)
pp


## -----------------------------------------------------------------------------
library(rgeos)
sp <- SpatialPoints(p, proj4string=CRS(proj4string(win)))
i <- gIntersects(sp, win, byid=TRUE)
which(!i)


## ---- pp23--------------------------------------------------------------------
plot(census)
points(sp)
points(sp[!i,], col='red', cex=3, pch=20)


## ---- pp24--------------------------------------------------------------------
pp <- ppp(p[i,1], p[i,2], window=owin, marks=census$dens[i])
plot(pp)
plot(city, add=TRUE)


## ---- pp25--------------------------------------------------------------------
s <- Smooth.ppp(pp)
plot(s)
plot(city, add=TRUE)


## ---- pp26, fig.width=9-------------------------------------------------------
par(mfrow=c(2,2), mai=c(0.25, 0.25, 0.25, 0.25))
for (offense in c("Auto Theft", "Drunk in Public", "DUI", "Arson")) {
  plot(city, col='grey')
	acrime <- crime[crime$CATEGORY == offense, ]
	points(acrime, col = "red")
	title(offense)
}


## -----------------------------------------------------------------------------
crime$fcat <- as.factor(crime$CATEGORY)
w <- as.owin(city)
xy <- coordinates(crime)
mpp <- ppp(xy[,1], xy[,2], window = w, marks=crime$fcat)


## ---- pp27, fig.width=9-------------------------------------------------------
par(mai=c(0,0,0,0))
spp <- split(mpp)
plot(spp[1:4], main='')


## ---- pp28, fig.width=9-------------------------------------------------------
plot(density(spp[1:4]), main='')


## -----------------------------------------------------------------------------
spatstat.options(checksegments = FALSE)
ktheft <- Kest(spp$"Auto Theft")
ketheft <- envelope(spp$"Auto Theft", Kest)
ktheft <- Kest(spp$"Arson")
ketheft <- envelope(spp$"Arson", Kest)


## ---- pp29, fig.width=10------------------------------------------------------
par(mfrow=c(1,2))
plot(ktheft)
plot(ketheft)


## -----------------------------------------------------------------------------
KS.arson <- cdf.test(spp$Arson, covariate=ds, test='ks')
KS.arson
KS.drunk <- cdf.test(spp$'Drunk in Public', covariate=ds, test='ks')
KS.drunk


## ---- pp30--------------------------------------------------------------------
kc <- Kcross(mpp, i = "Drunk in Public", j = "Arson")
ekc <- envelope(mpp, Kcross, nsim = 50, i = "Drunk in Public", j = "Arson")
plot(ekc)