## -----------------------------------------------------------------------------
set.seed(0)
xy <- cbind(x=runif(1000, 0, 100), y=runif(1000, 0, 100))
income <- (runif(1000) * abs((xy[,1] - 50) * (xy[,2] - 50))) / 500


## ----incplot, fig.width=8, fig.height=4---------------------------------------
par(mfrow=c(1,3), las=1)
plot(sort(income), col=rev(terrain.colors(1000)), pch=20, cex=.75, ylab='income')
hist(income, main='', col=rev(terrain.colors(10)),  xlim=c(0,5), breaks=seq(0,5,0.5))
plot(xy, xlim=c(0,100), ylim=c(0,100), cex=income, col=rev(terrain.colors(50))[10*(income+1)])


## -----------------------------------------------------------------------------
n <- length(income)
G <- (2 * sum(sort(income) * 1:n)/sum(income) - (n + 1)) / n
G


## -----------------------------------------------------------------------------
library(terra)
v <- vect(xy)
v$income <- income
r1 <- rast(ncol=1, nrow=4, xmin=0, xmax=100, ymin=0, ymax=100)
r1 <- rasterize(v, r1, "income", mean)

r2 <- rast(ncol=4, nrow=1, xmin=0, xmax=100, ymin=0, ymax=100)
r2 <- rasterize(v, r2, "income", mean)

r3 <- rast(ncol=2, nrow=2, xmin=0, xmax=100, ymin=0, ymax=100)
r3 <- rasterize(v, r3, "income", mean)

r4 <- rast(ncol=3, nrow=3, xmin=0, xmax=100, ymin=0, ymax=100)
r4 <- rasterize(v, r4, "income", mean)

r5 <- rast(ncol=5, nrow=5, xmin=0, xmax=100, ymin=0, ymax=100)
r5 <- rasterize(v, r5, "income", mean)

r6 <- rast(ncol=10, nrow=10, xmin=0, xmax=100, ymin=0, ymax=100)
r6 <- rasterize(v, r6, "income", mean)



## ----incdistplot, fig.height=5, fig.width=9-----------------------------------
par(mfrow=c(2,3), las=1)
plot(r1); plot(r2); plot(r3); plot(r4); plot(r5); plot(r6)


## ----inchist------------------------------------------------------------------
par(mfrow=c(1,3), las=1)
hist(r4, col=rev(terrain.colors(10)), xlim=c(0,5), breaks=seq(0, 5, 0.5))
hist(r5, main="", col=rev(terrain.colors(10)), xlim=c(0,5), breaks=seq(0, 5, 0.5))
hist(r6, main="", col=rev(terrain.colors(10)), xlim=c(0,5), breaks=seq(0, 5, 0.5))


## -----------------------------------------------------------------------------
A <- c(40, 43)
B <- c(101, 1)
C <- c(111, 54)
D <- c(104, 65)
E <- c(60, 22)
F <- c(20, 2)
pts <- rbind(A, B, C, D, E, F)
pts


## ----pointlabs----------------------------------------------------------------
plot(pts, xlim=c(0,120), ylim=c(0,120), pch=20, cex=2, col='red', xlab='X', ylab='Y', las=1)
text(pts+5, LETTERS[1:6])


## -----------------------------------------------------------------------------
dis <- dist(pts)
dis


## -----------------------------------------------------------------------------
sqrt((40-101)^2 + (43-1)^2)


## -----------------------------------------------------------------------------
D <- as.matrix(dis)
round(D)


## -----------------------------------------------------------------------------
gdis <- distance(pts, lonlat=TRUE)
gdis


## -----------------------------------------------------------------------------
a <-  D < 50
a


## -----------------------------------------------------------------------------
diag(a) <- NA
Adj50 <- a * 1
Adj50


## -----------------------------------------------------------------------------
cols <- apply(D, 1, order)
# we need to transpose the result
cols <- t(cols)


## -----------------------------------------------------------------------------
cols <- cols[, 2:3]
cols


## -----------------------------------------------------------------------------
rowcols <- cbind(rep(1:6, each=2), as.vector(t(cols)))
head(rowcols)


## -----------------------------------------------------------------------------
Ak3 <- Adj50 * 0
Ak3[rowcols] <- 1
Ak3


## -----------------------------------------------------------------------------
W <- 1 / D
round(W, 4)


## -----------------------------------------------------------------------------
W[!is.finite(W)] <- NA


## -----------------------------------------------------------------------------
rtot <- rowSums(W, na.rm=TRUE)
# this is equivalent to
# rtot <- apply(W, 1, sum, na.rm=TRUE)
rtot


## -----------------------------------------------------------------------------
W <- W / rtot
rowSums(W, na.rm=TRUE)


## -----------------------------------------------------------------------------
colSums(W, na.rm=TRUE)


## -----------------------------------------------------------------------------
p <- vect(system.file("ex/lux.shp", package="terra"))


## -----------------------------------------------------------------------------
wr <- adjacent(p, "rook", pairs=FALSE)
dim(wr)
wr[1:6,1:11]


## -----------------------------------------------------------------------------
i <- rowSums(wr)
i


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


## ----links, fig.width=6-------------------------------------------------------
par(mai=c(0,0,0,0))
plot(p, col="gray", border="blue")
nb <- adjacent(p, "rook")
v <- centroids(p)
p1 <- v[nb[,1], ]
p2 <- v[nb[,2], ]
lines(p1, p2, col="red", lwd=2)


## -----------------------------------------------------------------------------
wd10 <- nearby(v, distance=10000)
wd25 <- nearby(v, distance=25000)


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


## ----weights, fig.height=5, fig.width=9---------------------------------------
plotit <- function(nb, lab='') {
  plot(p, col='gray', border='white')
  v <- centroids(p)
  p1 <- v[nb[,1], ,drop=FALSE]
  p2 <- v[nb[,2], ,drop=FALSE]
  lines(p1, p2, col="red", lwd=2)
  text(6.3, 50.1, paste0('(', lab, ')'), cex=1.25)
}

par(mfrow=c(1, 3), mai=c(0,0,0,0))
plotit(nb, "adjacency")
plotit(wd25, "25 km")
plotit(k3, "k=3")