# \donttest{
p = c("sf", "sp", "terra", "gstat")
if(any(!unlist(lapply(p, requireNamespace, quietly=TRUE)))) {
m = which(!unlist(lapply(p, requireNamespace, quietly=TRUE)))
message("Can't run examples, please install ", paste(p[m], collapse = " "))
} else {
invisible(lapply(p, require, character.only=TRUE))
data(meuse, package = "sp")
meuse <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992,
agr = "constant")
data(meuse.grid, package = "sp")
meuse.grid <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992,
agr = "constant")
ref <- rast(ext(meuse.grid), resolution = 40)
crs(ref) <- crs(meuse)
e <- ext(179407.8, 181087.9, 331134.4, 332332.1)
# GRID-1 log(copper):
v1 <- variogram(log(copper) ~ 1, meuse)
x1 <- fit.variogram(v1, vgm(1, "Sph", 800, 1))
G1 <- krige(zinc ~ 1, meuse, meuse.grid, x1, nmax = 30)
G1 <- crop(rasterize(G1, ref, "var1.pred"),e)
names(G1) <- "copper"
# GRID-2 log(elev):
v2 <- variogram(log(elev) ~ 1, meuse)
x2 <- fit.variogram(v2, vgm(1, "Sph", 800, 1))
G2 <- krige(zinc ~ 1, meuse, meuse.grid, x2, nmax = 30)
G2 <- crop(rasterize(G2, ref, "var1.pred"),e)
names(G2) <- "elev"
# Raster corrected correlation
acor <- raster.modified.ttest(G1, G2)
plot(acor)
# Sample-based corrected correlation
( cor.hex <- raster.modified.ttest(G1, G2, sample = "hexagonal") )
plot(cor.hex["corr"], pch=20)
}
# }
Run the code above in your browser using DataLab