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


## -----------------------------------------------------------------------------
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)
library(raster)
sppol <- spPolygons(pol)


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


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


## -----------------------------------------------------------------------------
# make sure that the coordinates are interpreted as planar (not longitude/latitude)
crs(sppol) <- '+proj=utm +zone=1'
area(sppol)


## -----------------------------------------------------------------------------
library(raster)
usa <- raster::getData('GADM', country='USA', level=1)
usa <- usa[! usa$NAME_1 %in% c('Alaska', 'Hawaii'), ]


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


## -----------------------------------------------------------------------------
# patience, this takes a while: 
wus <- poly2nb(usa, row.names=usa$OBJECTID, queen=FALSE)
wus
wmus <- nb2mat(wus, style='B', zero.policy = TRUE)
dim(wmus)


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


## ---- auck1-------------------------------------------------------------------
if (!require("rspatial")) remotes::install_github('rspatial/rspatial')
library(rspatial)
pols <- sp_data("auctb.rds")


## ---- 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 <- poly2nb(pols, row.names=pols$Id, queen=FALSE)
class(wr)
summary(wr)


## ---- links, fig.width=6------------------------------------------------------
par(mai=c(0,0,0,0))
plot(pols, col='gray', border='blue')
xy <- coordinates(pols)
plot(wr, xy, col='red', lwd=2, add=TRUE)


## -----------------------------------------------------------------------------
wm <- nb2mat(wr, style='B')
dim(wm)


## -----------------------------------------------------------------------------
wr[1:6]
wm[1:6,1:11]


## -----------------------------------------------------------------------------
wq <- poly2nb(pols, row.names=pols$Id, queen=TRUE)


## -----------------------------------------------------------------------------
wd1 <- dnearneigh(xy, 0, 1000)
wd25 <- dnearneigh(xy, 0, 2500)


## -----------------------------------------------------------------------------
k3 <- knn2nb(knearneigh(xy, k=3))
k6 <- knn2nb(knearneigh(xy, k=6))


## -----------------------------------------------------------------------------
library(deldir)
d <- deldir(xy[,1], xy[,2], suppressMsge=TRUE)


## -----------------------------------------------------------------------------
wr2 <- wr
for (i in 1:length(wr)) {
	lag1 <- wr[[i]]
	lag2 <- wr[lag1]
	lag2 <- sort(unique(unlist(lag2)))
	lag2 <- lag2[!(lag2 %in% c(wr[[i]], i))] 
	wr2[[i]] <- lag2
}


## ---- weights, fig.height=12, fig.width=9-------------------------------------
  plotit <- function(nb, lab='') {
  plot(pols, col='gray', border='white')
  plot(nb, xy, add=TRUE, pch=20)
  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')
plot(d, wlines='triang', add=TRUE, pch=20)
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, n)
yiyj <- yi * yj


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


## -----------------------------------------------------------------------------
pmw <- pm * wm
wm[1:6, 1:9]
round(pmw[1:6, 1:9])


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


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


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


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


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


## -----------------------------------------------------------------------------
ww <- nb2listw(wr, style='B')
ww


## -----------------------------------------------------------------------------
moran(pols$TB, ww, n=length(ww$neighbours), S0=Szero(ww))


## -----------------------------------------------------------------------------
Szero(ww)


## -----------------------------------------------------------------------------
sum(pmw != 0)


## -----------------------------------------------------------------------------
moran.test(pols$TB, ww, randomisation=FALSE)


## -----------------------------------------------------------------------------
moran.mc(pols$TB, ww, nsim=99)


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


## -----------------------------------------------------------------------------
rwm <- mat2listw(wm, style='W')
# Checking if rows add up to 1
mat <- listw2mat(rwm)
apply(mat, 1, sum)[1:15]


## ---- mplot-------------------------------------------------------------------
moran.plot(y, rwm)