R terra
SCALE AND DISTANCE
Create data
set.seed(0)
xy <- cbind(x=runif(1000, 0, 100), y=runif(1000, 0, 100))
head(xy)## x y
## [1,] 89.66972 26.55367
## [2,] 26.55087 53.08088
## [3,] 37.21239 68.48609
## [4,] 57.28534 38.32834
## [5,] 90.82078 95.49880
## [6,] 20.16819 11.83566income <- (runif(1000) * abs((xy[,1] - 50) * (xy[,2] - 50))) / 500
head(income)## [1] 0.7201148 0.1259653 0.4572771 0.1474312 1.6259342 0.4370467Income Plot
par(mfrow=c(1,3), las=1)
plot(sort(income), col=rev(terrain.colors(1000)), pch=20, cex=.75, ylab='income')
Income Summary Plot
hist(income, main='', col=rev(terrain.colors(10)), xlim=c(0,5), breaks=seq(0,5,0.5))
Plot lat long
plot(xy,
xlim = c(0, 100),
ylim = c(0, 100),
cex = income,
col = rev(terrain.colors(50))[pmin(50, 10 * (income + 1))],
main = "Scatter plot of XY")
Gene Coefficient
Income inequality is often expressed with the Gini coefficient.
n <- length(income)
G <- (2 * sum(sort(income) * 1:n)/sum(income) - (n + 1)) / n
G ## [1] 0.5814548## [1] 0.5814548 Create Districts
- Assume that the household data was grouped by some kind of census districts.
- Create different districts, in this case rectangular raster cells, and compute mean income for each district.
library(terra) ## terra 1.7.62 ## terra 1.7.83Create vector
v <- vect(xy)
v## class : SpatVector
## geometry : points
## dimensions : 1000, 0 (geometries, attributes)
## extent : 0.1314657, 99.99306, 0.06052661, 99.88775 (xmin, xmax, ymin, ymax)
## coord. ref. :v$income <- income
v## class : SpatVector
## geometry : points
## dimensions : 1000, 1 (geometries, attributes)
## extent : 0.1314657, 99.99306, 0.06052661, 99.88775 (xmin, xmax, ymin, ymax)
## coord. ref. :
## names : income
## type : <num>
## values : 0.7201
## 0.126
## 0.4573Create different district
district 1
r1 <- rast(ncol=1, nrow=4, xmin=0, xmax=100, ymin=0, ymax=100)
r1 <- rasterize(v, r1, "income", mean)
r1## class : SpatRaster
## dimensions : 4, 1, 1 (nrow, ncol, nlyr)
## resolution : 100, 25 (x, y)
## extent : 0, 100, 0, 100 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## name : mean
## min value : 0.2787015
## max value : 0.9272025district 2
r2 <- rast(ncol=4, nrow=1, xmin=0, xmax=100, ymin=0, ymax=100)
r2 <- rasterize(v, r2, "income", mean) district 3
r3 <- rast(ncol=2, nrow=2, xmin=0, xmax=100, ymin=0, ymax=100)
r3 <- rasterize(v, r3, "income", mean) district 4
r4 <- rast(ncol=3, nrow=3, xmin=0, xmax=100, ymin=0, ymax=100)
r4 <- rasterize(v, r4, "income", mean) district 5
r5 <- rast(ncol=5, nrow=5, xmin=0, xmax=100, ymin=0, ymax=100)
r5 <- rasterize(v, r5, "income", mean) district 6
r6 <- rast(ncol=10, nrow=10, xmin=0, xmax=100, ymin=0, ymax=100)
r6 <- rasterize(v, r6, "income", mean) Plot district data
par(mfrow=c(2,3), las=1)
plot(r1);
plot(r2);
plot(r3);
plot(r4);
plot(r5);
plot(r6) 
Plot
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)) 
Distance Matrix
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 ## [,1] [,2]
## A 40 43
## B 101 1
## C 111 54
## D 104 65
## E 60 22
## F 20 2Plot Distance
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]) 
distfunction
dis <- dist(pts)
dis ## A B C D E
## B 74.06079
## C 71.84706 53.93515
## D 67.67570 64.07027 13.03840
## E 29.00000 46.06517 60.20797 61.52235
## F 45.61798 81.00617 104.80935 105.00000 44.72136Check 1st point using using Pythagoras’ theorem.
sqrt((40-101)^2 + (43-1)^2) ## [1] 74.06079## [1] 74.06079 transform a distance matrix into a normal matrix
D <- as.matrix(dis)
D## A B C D E F
## A 0.00000 74.06079 71.84706 67.67570 29.00000 45.61798
## B 74.06079 0.00000 53.93515 64.07027 46.06517 81.00617
## C 71.84706 53.93515 0.00000 13.03840 60.20797 104.80935
## D 67.67570 64.07027 13.03840 0.00000 61.52235 105.00000
## E 29.00000 46.06517 60.20797 61.52235 0.00000 44.72136
## F 45.61798 81.00617 104.80935 105.00000 44.72136 0.00000round(D) ## A B C D E F
## A 0 74 72 68 29 46
## B 74 0 54 64 46 81
## C 72 54 0 13 60 105
## D 68 64 13 0 62 105
## E 29 46 60 62 0 45
## F 46 81 105 105 45 0Distance for longitude/latitude coordinates
gdis <- distance(pts, lonlat=TRUE)
gdis ## 1 2 3 4 5
## 2 7614198
## 3 5155577 5946748
## 4 4581656 7104895 1286094
## 5 2976166 5011592 5536367 5737063
## 6 4957298 9013726 9894640 9521864 4859627Adjacency
a <- D < 50
a ## A B C D E F
## A TRUE FALSE FALSE FALSE TRUE TRUE
## B FALSE TRUE FALSE FALSE TRUE FALSE
## C FALSE FALSE TRUE TRUE FALSE FALSE
## D FALSE FALSE TRUE TRUE FALSE FALSE
## E TRUE TRUE FALSE FALSE TRUE TRUE
## F TRUE FALSE FALSE FALSE TRUE TRUEdiag(a) <- NA
Adj50 <- a * 1
Adj50 ## A B C D E F
## A NA 0 0 0 1 1
## B 0 NA 0 0 1 0
## C 0 0 NA 1 0 0
## D 0 0 1 NA 0 0
## E 1 1 0 0 NA 1
## F 1 0 0 0 1 NATwo nearest neighbours
cols <- apply(D, 1, order)
# we need to transpose the result
cols <- t(cols)
cols <- cols[, 2:3]
cols ## [,1] [,2]
## A 5 6
## B 5 3
## C 4 2
## D 3 5
## E 1 6
## F 5 1row-column pairs that we want (rowcols)
rowcols <- cbind(rep(1:6, each=2), as.vector(t(cols)))
rowcols## [,1] [,2]
## [1,] 1 5
## [2,] 1 6
## [3,] 2 5
## [4,] 2 3
## [5,] 3 4
## [6,] 3 2
## [7,] 4 3
## [8,] 4 5
## [9,] 5 1
## [10,] 5 6
## [11,] 6 5
## [12,] 6 1Use these pairs as indices to change the values in matrix Ak3
Ak3 <- Adj50 * 0
Ak3[rowcols] <- 1
Ak3 ## A B C D E F
## A NA 0 0 0 1 1
## B 0 NA 1 0 1 0
## C 0 1 NA 1 0 0
## D 0 0 1 NA 1 0
## E 1 0 0 0 NA 1
## F 1 0 0 0 1 NAWeights matrix
W <- 1 / D
round(W, 4) ## A B C D E F
## A Inf 0.0135 0.0139 0.0148 0.0345 0.0219
## B 0.0135 Inf 0.0185 0.0156 0.0217 0.0123
## C 0.0139 0.0185 Inf 0.0767 0.0166 0.0095
## D 0.0148 0.0156 0.0767 Inf 0.0163 0.0095
## E 0.0345 0.0217 0.0166 0.0163 Inf 0.0224
## F 0.0219 0.0123 0.0095 0.0095 0.0224 InfW[!is.finite(W)] <- NA
rtot <- rowSums(W, na.rm=TRUE)
rtot ## A B C D E F
## 0.09860117 0.08170418 0.13530597 0.13285878 0.11141516 0.07569154W <- W / rtot
rowSums(W, na.rm=TRUE) ## A B C D E F
## 1 1 1 1 1 1colSums(W, na.rm=TRUE) ## A B C D E F
## 0.9784548 0.7493803 1.2204900 1.1794393 1.1559273 0.7163082Spatial influence for polygons
p <- vect(system.file("ex/lux.shp", package="terra")) wr <- adjacent(p, "rook", pairs=FALSE)
dim(wr) ## [1] 12 12i <- rowSums(wr)
i ## 1 2 3 4 5 6 7 8 9 10 11 12
## 3 6 4 2 3 3 3 4 4 3 5 6round(100 * table(i) / length(i), 1) ## i
## 2 3 4 5 6
## 8.3 41.7 25.0 8.3 16.7par(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], ]
plot(p, col="gray", border="blue")
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) 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") 