## ----ch7_1--------------------------------------------------------------------
library(terra)
pol <- matrix(c(1.7, 2.6, 5.6, 8.1, 7.2, 3.3, 1.7, 4.9, 7, 7.6, 6.1, 2.7, 2.7, 4.9), ncol=2)
sppol <- vect(pol, "polygons")


## -----------------------------------------------------------------------------
negpol <- rbind(pol[c(1,6:4),], cbind(pol[4,1], 0), cbind(pol[1,1], 0))
spneg <- vect(negpol, "polygons")


## ----polygons-----------------------------------------------------------------

cols <- c('light gray', 'light blue')
plot(sppol, xlim=c(1,9), ylim=c(1,10), col=cols[1], axes=FALSE, xlab='', ylab='', 
      lwd=2, yaxs="i", xaxs="i")
plot(spneg, col=cols[2], add=T)
plot(spneg, add=T, density=8, angle=-45, lwd=1)
segments(pol[,1], pol[,2], pol[,1], 0)
text(pol, LETTERS[1:6], pos=3, col='red', font=4)
arrows(1, 1, 9, 1, 0.1, xpd=T)
arrows(1, 1, 1, 9, 0.1, xpd=T)
text(1, 9.5, 'y axis', xpd=T)
text(10, 1, 'x axis', xpd=T)
legend(6, 9.5, c('"positive" area', '"negative" area'), fill=cols, bty = "n")


## -----------------------------------------------------------------------------
p <- rbind(pol, pol[1,])
x <- p[-1,1] - p[-nrow(p),1]
y <- (p[-1,2] + p[-nrow(p),2]) / 2
sum(x * y)


## -----------------------------------------------------------------------------
crs(sppol) <- '+proj=utm +zone=1'
expanse(sppol)


## -----------------------------------------------------------------------------
if (!require(geodata)) remotes::install_github('rspatial/geodata')
usa <- geodata::gadm(country='USA', level=1, path=".")
usa <- usa[! usa$NAME_1 %in% c('Alaska', 'Hawaii'), ]


## -----------------------------------------------------------------------------
# patience, this takes a while: 
wus <- relate(usa, relation="touches")
rownames(wus) <- colnames(wus) <- usa$NAME_1
wus[1:5,1:5]


## -----------------------------------------------------------------------------
i <- rowSums(wus)
round(100 * table(i) / length(i), 1)


## ----auck1--------------------------------------------------------------------
if (!require("rspat")) remotes::install_github('rspatial/rspat')
library(rspat)
pols <- spat_data("auctb")


## ----auck2, fig.width=8-------------------------------------------------------
par(mai=c(0,0,0,0))
classes <- seq(0,450,50)
cuts <- cut(pols$TB, classes)
n <- length(classes)
cols <- rev(gray(0:n / n))
plot(pols, col=cols[as.integer(cuts)])
legend('bottomleft', levels(cuts), fill=cols)


## -----------------------------------------------------------------------------
wr <- adjacent(pols, type="rook", symmetrical=TRUE)
head(wr)


## ----links, fig.width=6-------------------------------------------------------
v <- centroids(pols)
p1 <- v[wr[,1], ]
p2 <- v[wr[,2], ]

par(mai=c(0,0,0,0))
plot(pols, col='gray', border='blue')
lines(p1, p2, col='red', lwd=2)
points(v)


## -----------------------------------------------------------------------------
wq <- adjacent(pols, "queen", symmetrical=TRUE)


## -----------------------------------------------------------------------------
wd1 <- nearby(pols, distance=1000)
wd25 <- nearby(pols, distance=2500)


## -----------------------------------------------------------------------------
k3 <- nearby(pols, k=3)
k6 <- nearby(pols, k=6)


## -----------------------------------------------------------------------------
d <- delaunay(centroids(pols))


## -----------------------------------------------------------------------------
wrs <- adjacent(pols, "rook", symmetrical=FALSE)
uf <- sort(unique(wrs[,1]))
wr2 <- list()
for (i in 1:length(pols)) {
	lag1 <- wrs[wrs[,1]==i, 2]
  lag2 <- wrs[wrs[,1] %in% lag1, ] 	
  lag2[,1] <- i
	wr2[[i]] <- unique(lag2)
}

wr2 <- do.call(rbind, wr2)



## ----weights, fig.height=12, fig.width=9--------------------------------------
plotit <- function(nb, lab='') {
  plot(pols, col='gray', border='white')
  v <- centroids(pols)
  p1 <- v[nb[,1], ]
  p2 <- v[nb[,2], ]
  lines(p1, p2, col='red', lwd=2)
  points(v)
  text(2659066, 6482808, paste0('(', lab, ')'), cex=1.25)
}

par(mfrow=c(4, 2), mai=c(0,0,0,0))
plotit(wr, 'i')
plotit(wq, 'ii')
plotit(wd1, 'iii')
plotit(wd25, 'iv')
plotit(k3, 'v')
plotit(k6, 'vi')
plot(pols, col='gray', border='white')
lines(d, col="red")
text(2659066, 6482808, '(vii)', cex=1.25)
plotit(wr2, 'viii')


## -----------------------------------------------------------------------------
n <- length(pols)


## -----------------------------------------------------------------------------
y <- pols$TB
ybar <- mean(y)


## -----------------------------------------------------------------------------
dy <- y - ybar
g <- expand.grid(dy, dy)
yiyj <- g[,1] * g[,2]


## -----------------------------------------------------------------------------
yi <- rep(dy, each=n)
yj <- rep(dy)
yiyj <- yi * yj


## -----------------------------------------------------------------------------
pm <- matrix(yiyj, ncol=n)
round(pm[1:6, 1:9])


## -----------------------------------------------------------------------------
wm <- adjacent(pols, "rook", pairs=FALSE)
wm[1:9, 1:11]
pmw <- pm * wm
round(pmw[1:9, 1:11])


## -----------------------------------------------------------------------------
spmw <- sum(pmw) 
spmw


## -----------------------------------------------------------------------------
smw <- sum(wm)
sw  <- spmw / smw


## -----------------------------------------------------------------------------
vr <- n / sum(dy^2)


## -----------------------------------------------------------------------------
MI <- vr * sw
MI


## -----------------------------------------------------------------------------
autocor(y, wm)


## -----------------------------------------------------------------------------
  EI <- -1/(n-1)
EI


## -----------------------------------------------------------------------------
I <- autocor(pols$TB, wm)
nsim <- 99
mc <- sapply(1:nsim, function(i) autocor(sample(pols$TB), wm))

P <- 1 - sum((I > mc) / (nsim+1))
P


## -----------------------------------------------------------------------------
n <- length(pols)
ms <- cbind(id=rep(1:n, each=n), y=rep(y, each=n), value=as.vector(wm * y))


## -----------------------------------------------------------------------------
ms <- ms[ms[,3] > 0, ]


## -----------------------------------------------------------------------------
ams <- aggregate(ms[,2:3], list(ms[,1]), FUN=mean)
ams <- ams[,-1]
colnames(ams) <- c('y', 'spatially lagged y')
head(ams)


## ----ngb----------------------------------------------------------------------
plot(ams)
reg <- lm(ams[,2] ~ ams[,1])
abline(reg, lwd=2)
abline(h=mean(ams[,2]), lt=2)
abline(v=ybar, lt=2)


## ----ngb2---------------------------------------------------------------------
coefficients(reg)[2]