# -- Load packages -- #
library("Spbsampling")
library("sp")
library("sf")
library("sampling")
library("plot3D")
data("lucas_abruzzo", package = "Spbsampling")
data("meuse", package = "sp")
data("meuse.grid", package = "sp")
load("maps.RData")

# ------------- #
# -- Section 5 - Population plots -- #
par(mfrow = c(1, 2))
abruzzo <- lucas_abruzzo
abruzzo$prov <- as.vector(abruzzo$prov)
abruzzo$lu <- as.vector(abruzzo$lu)
abruzzo$lc <- as.vector(abruzzo$lc)
pr <- provmp[provmp@data$COD_REG == 13, ]
re <- regmp[regmp@data$COD_REG == 13, ]
df <- abruzzo
df$lc[df$lc %in% c("A11", "A12", "A13", "A21", "A22")]                                             <- "A"
df$lc[df$lc %in% c("B11", "B12")]                                                                  <- "B11"
df$lc[df$lc %in% c("B13", "B14", "B15", "B16", "B17", "B18", "B19")]                               <- "B19"
df$lc[df$lc %in% c("B21", "B22", "B23", "B31", "B32", "B33", "B37", "B41", "B42", "B43", "B45")]   <- "B9"
df$lc[df$lc %in% c("B51", "B52", "B53", "B54", "B55")]                                             <- "B5"
df$lc[df$lc %in% c("B71", "B72", "B73", "B74", "B75", "B77", "B81", "B82", "B83")]                 <- "B7"
df$lc[df$lc %in% c("C10", "C21", "C22", "C23", "C31", "C32", "C33")]                               <- "C"
df$lc[df$lc %in% c("D10", "D20", "F10", "F20", "F40", "G10", "G20", "G30")]                        <- "G"
df$lc[df$lc %in% c("E10", "E20", "E30")]                                                           <- "E"

df$lc <- factor(df$lc,
                labels = c("Artif. Land", "Wheat", "Other Cereals",
                           "Fodder Crops", "Permanent Crops",
                           "Other Cropland", "Woodland", "Grassland",
                           "Other Land Cover"))

df$lu[df$lu %in% c("U111", "U112", "U113")]                                    <- "U11"
df$lu[df$lu %in% c("U120")]                                                    <- "U12"
df$lu[df$lu %in% c("U130", "U140", "U210", "U221", "U223", "U224", "U225", "U226", "U227")] <- "U5"
df$lu[df$lu %in% c("U311", "U312", "U313", "U315", "U317", "U318", "U321", "U322", "U330", "U340", "U350", "U361", "U362", "U370")] <- "U3"
df$lu[df$lu %in% c("U410", "U420")]                                           <- "U4"

df$lu   <- factor(df$lu, labels = c("Agriculture", "Forestry", "Urban", "Unused", "Other Land Use"))
df$prov <- factor(df$prov, labels = c("L'Aquila", "Teramo", "Chieti", "Pescara"))

abruzzo <- df
C1 <- (abruzzo$x - 1000000) / 2000 - 416
R1 <- (abruzzo$y) / 2000 - 2313
lancov <- matrix(NA_real_, 67, 73)
cods <- as.numeric(abruzzo$lc)
for (i in 1:nrow(abruzzo)) {
  lancov[R1[i], C1[i]] <- cods[i] 
}
cods <- as.numeric(abruzzo$lu)
lanuse <- matrix(NA_real_, 67, 73)
for (i in 1:nrow(abruzzo)) {
  lanuse[R1[i], C1[i]] <- cods[i] 
}
par(mar = c(0, 0, 0, 0))
plot(provmp, xlim = c(834000, 974000), ylim = c(4628000, 4760000))
image(seq(834000, 978000, 2000), seq(4628000, 4760000, 2000), t(lancov), col = gray.colors(9), add = TRUE)
plot(provmp, lwd = 2, add = TRUE)
plot(regmp, lwd = 4, add = TRUE)
legend(930000, 4810000, title = "Abruzzo - Land Cover",
       c("Artif. Land", "Wheat", "Other Cereals", "Fodder Crops",
         "Permanent Crops", "Other Cropland", "Woodland", "Grassland",
         "Other Land Cover"), fill = gray.colors(11), cex = 0.4, ncol = 1)
box()
data("meuse.grid", package = "sp")
data("meuse.riv", package = "sp")
coordinates(meuse.grid) <- ~x+y
gridded(meuse.grid) <- TRUE
plot(geometry(meuse.grid), asp = 1)
lines(meuse.riv, type = "l", asp = 1)
polygon(meuse.riv,col = "gray")
box()
par(mfrow = c(1, 1))

# ---------------------------------- #
# -- Section 5.1. - Overview of the package use -- #
library("Spbsampling")
data("lucas_abruzzo", package = "Spbsampling")
dis_la <- as.matrix(dist(cbind(lucas_abruzzo$x, lucas_abruzzo$y)))
lucas_abruzzo_sf <- st_as_sf(lucas_abruzzo, coords = c("x", "y"))
dis_la_sf <- st_distance(lucas_abruzzo_sf)
# --------- #
# -- PWD -- #
con <- rep(0, nrow(dis_la))
stand_dist_la_pwd <- stprod(mat = dis_la, con = con)$mat # it takes some minutes
set.seed(42)
s_pwd_la <- pwd(dis = stand_dist_la_pwd, n = 10)$s
s_pwd_la
lucas_abruzzo[s_pwd_la[1, ], c("x", "y")]
pi <- rep(10 / nrow(lucas_abruzzo), nrow(lucas_abruzzo))
sbi(dis = dis_la, pi = pi, s = s_pwd_la[1, ])
set.seed(42)
s2_pwd_la <- pwd(dis = stand_dist_la_pwd, n = 10, nrepl = 2)$s
s2_pwd_la
# --------- #
# -- SWD -- #
con <- rep(1, nrow(dis_la))
stand_dist_la_swd <- stsum(mat = dis_la, con = con)$mat
set.seed(42)
s_swd_la <- swd(dis = stand_dist_la_swd, n = 10)$s
s_swd_la 
pi <- rep(10 / nrow(lucas_abruzzo), nrow(lucas_abruzzo))
sbi(dis = dis_la, pi = pi, s = s_swd_la[1, ])
# --------- #
# -- HPWD -- #
set.seed(42)
s_hpwd_la <- hpwd(dis = stand_dist_la_pwd, n = 10)
s_hpwd_la
pi <- rep(10 / nrow(lucas_abruzzo), nrow(lucas_abruzzo))
sbi(dis = dis_la, pi = pi, s = s_hpwd_la[1, ])
# ---------- #
# -- plot -- #
par(mfrow = c(1, 3))
plot(provmp[provmp@data$COD_REG == 13, ], asp = 1, main = "PWD")
plot(regmp[regmp@data$COD_REG == 13, ], asp = 1, lwd = 2, add = TRUE)
points(cbind(lucas_abruzzo[s_pwd_la, 4] - 1000000, lucas_abruzzo[s_pwd_la, 5]), asp = 1, pch = 19)
plot(provmp[provmp@data$COD_REG == 13, ], asp = 1, main = "SWD")
plot(regmp[regmp@data$COD_REG == 13, ], asp = 1, lwd = 2, add = TRUE)
points(cbind(lucas_abruzzo[s_swd_la, 4] - 1000000, lucas_abruzzo[s_swd_la, 5]), asp = 1, pch = 19)
plot(provmp[provmp@data$COD_REG == 13, ], asp = 1, main = "HPWD")
plot(regmp[regmp@data$COD_REG == 13, ], asp = 1, lwd = 2, add = TRUE)
points(cbind(lucas_abruzzo[s_hpwd_la, 4] - 1000000, lucas_abruzzo[s_hpwd_la, 5]), asp = 1, pch = 19)
par(mfrow = c(1, 1))
# ---------- #
# -- Table 1 -- #
# Warning: it takes time (this is a table about computing time)
N <- seq(1000, 10000, by = 1000)
stprod_mat <- list()
stsum_mat <- list()
stprod_time <- rep(0, length(N))
stsum_time <- rep(0, length(N))
set.seed(42)
dis <- as.matrix(dist(cbind(runif(1000), runif(1000))))
stsum_time[1] <- system.time(stsum(mat = dis, con = rep(1, nrow(dis))))[3]
stprod_time[1] <- system.time(stprod(mat = dis, con = rep(0, nrow(dis))))[3]
set.seed(42)
dis <- as.matrix(dist(cbind(runif(2000), runif(2000))))
stsum_time[2] <- system.time(stsum(mat = dis, con = rep(1, nrow(dis))))[3]
stprod_time[2] <- system.time(stprod(mat = dis, con = rep(0, nrow(dis))))[3]
set.seed(42)
dis <- as.matrix(dist(cbind(runif(3000), runif(3000))))
stsum_time[3] <- system.time(stsum(mat = dis, con = rep(1, nrow(dis))))[3]
stprod_time[3] <- system.time(stprod(mat = dis, con = rep(0, nrow(dis))))[3]
set.seed(42)
dis <- as.matrix(dist(cbind(runif(4000), runif(4000))))
stsum_time[4] <- system.time(stsum(mat = dis, con = rep(1, nrow(dis))))[3]
stprod_time[4] <- system.time(stprod(mat = dis, con = rep(0, nrow(dis))))[3]
set.seed(42)
dis <- as.matrix(dist(cbind(runif(5000), runif(5000))))
stsum_time[5] <- system.time(stsum(mat = dis, con = rep(1, nrow(dis))))[3]
stprod_time[5] <- system.time(stprod(mat = dis, con = rep(0, nrow(dis))))[3]
set.seed(42)
dis <- as.matrix(dist(cbind(runif(6000), runif(6000))))
stsum_time[6] <- system.time(stsum(mat = dis, con = rep(1, nrow(dis))))[3]
stprod_time[6] <- system.time(stprod(mat = dis, con = rep(0, nrow(dis))))[3]
set.seed(42)
dis <- as.matrix(dist(cbind(runif(7000), runif(7000))))
stsum_time[7] <- system.time(stsum(mat = dis, con = rep(1, nrow(dis))))[3]
stprod_time[7] <- system.time(stprod(mat = dis, con = rep(0, nrow(dis))))[3]
set.seed(42)
dis <- as.matrix(dist(cbind(runif(8000), runif(8000))))
stsum_time[8] <- system.time(stsum(mat = dis, con = rep(1, nrow(dis))))[3]
stprod_time[8] <- system.time(stprod(mat = dis, con = rep(0, nrow(dis))))[3]
set.seed(42)
dis <- as.matrix(dist(cbind(runif(9000), runif(9000))))
stsum_time[9] <- system.time(stsum(mat = dis, con = rep(1, nrow(dis))))[3]
stprod_time[9] <- system.time(stprod(mat = dis, con = rep(0, nrow(dis))))[3]
set.seed(42)
dis <- as.matrix(dist(cbind(runif(10000), runif(10000))))
stsum_time[10] <- system.time(stsum(mat = dis, con = rep(1, nrow(dis))))[3]
stprod_time[10] <- system.time(stprod(mat = dis, con = rep(0, nrow(dis))))[3]
table1 <- rbind(stsum_time, stprod_time)
colnames(table1) <- c("N = 1000", "N = 2000", "N = 3000", "N = 4000", "N = 5000",
                      "N = 6000", "N = 7000", "N = 8000", "N = 9000", "N = 10000")
rownames(table1) <- c("stsum()", "stprod()")
table1
# ------------- #
# ------------------------------------------------ #
# -- Section 5.2 - The role of beta -- #
library("sp")
data("meuse.grid", package = "sp")
dis_m <- as.matrix(dist(cbind(meuse.grid$x, meuse.grid$y)))
# --------- #
# -- PWD -- #
con <- rep(0, nrow(meuse.grid))
stand_dist_m_pwd <- stprod(mat = dis_m, con = con)$mat # it takes some minutes
set.seed(42)
s_pwd1 <- pwd(dis = stand_dist_m_pwd, n = 20, beta = 1)$s
set.seed(42)
s_pwd5 <- pwd(dis = stand_dist_m_pwd, n = 20, beta = 5)$s
set.seed(42)
s_pwd10 <- pwd(dis = stand_dist_m_pwd, n = 20, beta = 10)$s
pi <- rep(20 / nrow(meuse.grid), nrow(meuse.grid))
sbi(dis = dis_m, pi = pi, s = s_pwd1[1, ])
sbi(dis = dis_m, pi = pi, s = s_pwd5[1, ])
sbi(dis = dis_m, pi = pi, s = s_pwd10[1, ])
# ---------- #
# -- HPWD -- #
set.seed(42)
s_hpwd1 <- hpwd(dis = stand_dist_m_pwd, n = 20, beta = 1)
set.seed(42)
s_hpwd5 <- hpwd(dis = stand_dist_m_pwd, n = 20, beta = 5)
set.seed(42)
s_hpwd10 <- hpwd(dis = stand_dist_m_pwd, n = 20, beta = 10)
pi <- rep(20 / nrow(meuse.grid), nrow(meuse.grid))
sbi(dis = dis_m, pi = pi, s = s_hpwd1[1, ])
sbi(dis = dis_m, pi = pi, s = s_hpwd5[1, ])
sbi(dis = dis_m, pi = pi, s = s_hpwd10[1, ])
# --------- #
# -- SWD -- #
con <- rep(1, nrow(meuse.grid))
stand_dist_m_swd <- stsum(mat = dis_m, con = con)$mat
set.seed(42)
s_swd1 <- swd(dis = stand_dist_m_swd, n = 20, beta = 1)$s
set.seed(42)
s_swd5 <- swd(dis = stand_dist_m_swd, n = 20, beta = 5)$s
set.seed(42)
s_swd10 <- swd(dis = stand_dist_m_swd, n = 20, beta = 10)$s
pi <- rep(20 / nrow(meuse.grid), nrow(meuse.grid))
sbi(dis = dis_m, pi = pi, s = s_swd1[1, ])
sbi(dis = dis_m, pi = pi, s = s_swd5[1, ])
sbi(dis = dis_m, pi = pi, s = s_swd10[1, ])
# ---------- #
# -- plot -- #
data("meuse.area", package = "sp")
data("meuse.riv", package = "sp")
data("meuse.grid", package = "sp")
par(mfrow = c(3, 3), mar = c(2, 1, 2, 1))
plot(meuse.area, axes = F, type = "l", asp = 1, xlab = "", ylab = "", main = "PWD beta = 1")
lines(meuse.riv, type = "l")
polygon(meuse.riv, col = "gray")
points(meuse.grid[s_pwd1[1, ], 1:2], pch = 20)
plot(meuse.area, axes = F, asp = 1, type = "l", xlab = "", ylab = "", main = "PWD beta = 5")
lines(meuse.riv, type = "l")
polygon(meuse.riv, col = "gray")
points(meuse.grid[s_pwd5[1, ], 1:2], pch = 20)
plot(meuse.area, axes = F, asp = 1, type = "l", xlab = "", ylab = "", main = "PWD beta = 10")
lines(meuse.riv, type = "l")
polygon(meuse.riv, col = "gray")
points(meuse.grid[s_pwd10[1, ], 1:2], pch = 20)
plot(meuse.area, axes = F, asp = 1, type = "l", xlab = "", ylab = "", main = "HPWD beta = 1")
lines(meuse.riv, type = "l")
polygon(meuse.riv, col = "gray")
points(meuse.grid[s_hpwd1[1, ], 1:2], pch = 20)
plot(meuse.area,axes=F, asp = 1, type = "l", xlab = "", ylab = "", main = "HPWD beta = 5")
lines(meuse.riv, type = "l")
polygon(meuse.riv, col = "gray")
points(meuse.grid[s_hpwd5[1, ], 1:2], pch = 20)
plot(meuse.area,axes = F, asp = 1, type = "l", xlab = "", ylab = "", main = "HPWD beta = 10")
lines(meuse.riv, type = "l", asp = 1)
polygon(meuse.riv, col = "gray")
points(meuse.grid[s_hpwd10[1, ], 1:2], pch = 20)
plot(meuse.area,axes = F, asp = 1, type = "l", xlab = "", ylab = "", main = "SWD beta = 1")
lines(meuse.riv, type = "l", asp = 1)
polygon(meuse.riv, col = "gray")
points(meuse.grid[s_swd1[1, ], 1:2], pch = 20)
plot(meuse.area,axes=F, asp = 1, type = "l", xlab = "", ylab = "", main = "SWD beta = 5")
lines(meuse.riv, type = "l", asp = 1)
polygon(meuse.riv, col = "gray")
points(meuse.grid[s_swd5[1, ], 1:2], pch = 20)
plot(meuse.area,axes = F, asp = 1, type = "l", xlab = "", ylab = "", main = "SWD beta = 10")
lines(meuse.riv, type = "l")
polygon(meuse.riv, col = "gray")
points(meuse.grid[s_swd10[1, ], 1:2], pch = 20)
par(mfrow = c(1, 1))
# ---------- #
# ------------------------------------ #
# -- Section 5.3 - Spread on the covariate space -- #
data("meuse", package = "sp")
dis_m_3 <- as.matrix(dist(cbind(meuse$copper, meuse$lead, meuse$zinc)))
con <- rep(0, nrow(meuse))
stand_dis_m_3 <- stprod(mat = dis_m_3, con = con)$mat
set.seed(42)
s <- pwd(dis = stand_dis_m_3, n = 10)$s
sbi(dis = dis_m_3, pi = rep(10 / nrow(meuse), nrow(meuse)), s = s)
scatter3D(meuse$copper, meuse$lead, meuse$zinc, pch = 20, colvar = NULL)
points3D(meuse[s, 4], meuse[s, 5], meuse[s, 6], pch = 20, col = "red", add = TRUE)
# ------------------------------------------------- #