## ----getDatapp---------------------------------------------------------------- if (!require("rspat")) remotes::install_github("rspatial/rspat") library(rspat) city <- spat_data("city") crime <- spat_data("crime") ## ----pp1, fig.width=8, fig.height=4------------------------------------------- plot(city, col="light blue") points(crime, col="red", cex=.5, pch="+") ## ----------------------------------------------------------------------------- tb <- sort(table(crime$CATEGORY))[-1] tb ## ----------------------------------------------------------------------------- xy <- crds(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 <- expanse(city) dens <- nrow(xy) / CityArea ## ----------------------------------------------------------------------------- r <- rast(city, res=1000) ## ----pp3, fig.width=8, fig.height=4------------------------------------------- r <- rasterize(city, r) plot(r) quads <- as.polygons(r) plot(quads, add=TRUE) points(crime, col='red', cex=.5) ## ----pp4, fig.width=8, fig.height=4------------------------------------------- nc <- rasterize(crime, r, fun=function(i){length(i)}, 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) 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--------------------------------------------------------------------- c# get the centers of the 'quadrats' (raster cells) p <- as.points(r) # compute distance from all crime sites to these cell centers d2 <- distance(p, crime) d2 <- as.matrix(d2) # 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", ylim=c(0,1.1)) 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------------------------------------------------------------ cityOwin <- as.owin(sf::st_as_sf(city)) class(cityOwin) cityOwin ## ----------------------------------------------------------------------------- pts <- terra::crds(crime) head(pts) ## ----pp20--------------------------------------------------------------------- p <- ppp(pts[,1], pts[,2], window=cityOwin) class(p) p plot(p) ## ----pp21--------------------------------------------------------------------- ds <- density(p) class(ds) plot(ds, main='crime density') ## ----------------------------------------------------------------------------- nrow(pts) r <- rast(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 <- spat_data("census2000.rds") ## ----------------------------------------------------------------------------- census$area <- expanse(census) census$area <- census$area/27878400 census$dens <- census$POP2000 / census$area ## ----------------------------------------------------------------------------- p <- terra::crds(centroids(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(sf::st_as_sf(win)) pp <- ppp(p[,1], p[,2], window=owin, marks=census$dens) pp ## ----------------------------------------------------------------------------- sp <- vect(p, crs=crs(win)) i <- relate(sp, win, "intersects") i <- which(!i) 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(sf::st_as_sf(city)) xy <- terra::crds(crime) mpp <- ppp(xy[,1], xy[,2], window = w, marks=as.factor(crime$fcat)) ## ----pp27, fig.width=9-------------------------------------------------------- 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, ds) KS.arson KS.drunk <- cdf.test(spp$'Drunk in Public', ds) 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)