## -----------------------------------------------------------------------------
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------------------------------------------------------------------
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(raster)
r1 <- raster(ncol=1, nrow=4, xmn=0, xmx=100, ymn=0, ymx=100, crs=NA)
r1 <- rasterize(xy, r1, income, mean)

r2 <- raster(ncol=4, nrow=1, xmn=0, xmx=100, ymn=0, ymx=100, crs=NA)
r2 <- rasterize(xy, r2, income, mean)

r3 <- raster(ncol=2, nrow=2, xmn=0, xmx=100, ymn=0, ymx=100, crs=NA)
r3 <- rasterize(xy, r3, income, mean)

r4 <- raster(ncol=3, nrow=3, xmn=0, xmx=100, ymn=0, ymx=100, crs=NA)
r4 <- rasterize(xy, r4, income, mean)

r5 <- raster(ncol=5, nrow=5, xmn=0, xmx=100, ymn=0, ymx=100, crs=NA)
r5 <- rasterize(xy, r5, income, mean)

r6 <- raster(ncol=10, nrow=10, xmn=0, xmx=100, ymn=0, ymx=100, crs=NA)
r6 <- rasterize(xy, r6, income, mean)



## ----incdistplot, fig.height = 7, 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, main='', 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)


## -----------------------------------------------------------------------------
library(raster)
gdis <- pointDistance(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)


## -----------------------------------------------------------------------------
library(raster)
p <- shapefile(system.file("external/lux.shp", package="raster"))


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


## -----------------------------------------------------------------------------
wr <- poly2nb(p, row.names=p$ID_2, queen=FALSE)
wr
wm <- nb2mat(wr, style='B', zero.policy = TRUE)
dim(wm)


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


## -----------------------------------------------------------------------------
i <- rowSums(wm)
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')
xy <- coordinates(p)
plot(wr, xy, col='red', lwd=2, add=TRUE)


## -----------------------------------------------------------------------------
wd10 <- dnearneigh(xy, 0, 10)
wd25 <- dnearneigh(xy, 0, 25, longlat=TRUE)


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


## -----------------------------------------------------------------------------
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=11, fig.width=9-------------------------------------
plotit <- function(nb, lab='') {
  plot(p, col='gray', border='white')
  plot(nb, xy, add=TRUE, pch=20)
  text(6.3, 50.1, paste0('(', lab, ')'), cex=1.25)
}

par(mfrow=c(2, 3), mai=c(0,0,0,0))
plotit(wr, 'adjacency')
plotit(wr2, 'lag-2 adj.')
plotit(wd10, '10 km')
plotit(wd25, '25 km')
plotit(k3, 'k=3')
plotit(k6, 'k=6')