本篇是介绍克里金插值的第三篇推文,也是最后一篇。
因为前面两篇使用的数据中已知点的样本太少,本篇使用gstat
工具包说明文档中的数据集。该数据集来自sp
工具包。
实际上,gstat
工具包的方法对sp和sf两种格式的空间矢量对象同样适用,且使用方法也一致。gstat
的说明文档中是以sp为例的,而本系列的推文是以sf对象为例。
加载数据集:
library(sp)
data("meuse")
data("meuse.area")
class(meuse)
## [1] "data.frame"
names(meuse)
## [1] "x" "y" "cadmium" "copper" "lead" "zinc" "elev"
## [8] "dist" "om" "ffreq" "soil" "lime" "landuse" "dist.m"
meuse
是数据框结构,记录的是墨兹河附近区域的信息:
x和y:点的空间位置;
cadmium-zinc:镉、铜、铅、锌四种金属元素的含量;
elev:高度;
dist:距离墨兹河的距离。
前面已有推文介绍过通过数据框创建sf对象的方法了:
library(sf)
data <- st_as_sf(meuse, coords = c(x = "x",
y = "y"))
class(meuse.area)
## [1] "matrix"
colnames(meuse.area)
## [1] "x" "y"
meuse.area
是矩阵结构,记录的是区域边界的点坐标。这里将其转换为sf线对象,作为插值边界:
border <- st_multilinestring(list(meuse.area))
创建插值网格、可视化锌含量的空间分布:
grid <- st_make_grid(border, n = c(100, 100))
plot(data["zinc"], reset = F)
plot(border, add = T)
![](https://i-blog.csdnimg.cn/blog_migrate/ab4bee2203475fa1cefbabb711b7d78e.png)
加载相关工具包:
library(gstat)
library(tidyverse)
5 协同克里金
协同克里金(Cokriging)假设多种属性值相互之间存在相关性,因此可以同时对多个属性变量进行插值。协同克里金使用的函数是gstat
工具包的同名函数:
gstat(g, id, formula, locations, data, model = NULL, beta,
nmax = Inf, nmin = 0, omax = 0, maxdist = Inf, force = FALSE,
dummy = FALSE, set, fill.all = FALSE,
fill.cross = TRUE, variance = "identity", weights = NULL, merge,
degree = 0, vdist = FALSE, lambda = 1.0)
考虑到克里金插值要求变量满足空间平稳性,这里以四种金属元素含量的对数作为插值对象。计算已知点的协同半变异函数值的方法如下:
g <- gstat(NULL, "logCd", log(cadmium) ~ 1, data)
g <- gstat(g, "logCu", log(copper) ~ 1, data)
g <- gstat(g, "logPb", log(lead) ~ 1, data)
g <- gstat(g, "logZn", log(zinc) ~ 1, data)
vg <- variogram(g)
plot(vg$dist, vg$gamma)
![](https://i-blog.csdnimg.cn/blog_migrate/dbfcdaaa65abbbedc0ea08eb2bb253d6.png)
拟合协同半变异函数的方法如下:
vgm <- vgm(psill = 1, model = "Sph",
range = 800, nugget = 1)
fit <- fit.lmc(vg, g, vgm)
plot(vg, fit)
![](https://i-blog.csdnimg.cn/blog_migrate/4b7894367eab020180acf8196a1c9e8b.png)
插值过程使用predict
函数:
cokrige <- predict(fit, grid)
save(cokrige, file = "32-1.cokrige.rdata")
## Linear Model of Coregionalization found. Good.
## [using ordinary cokriging]
插值结果:
load("32-1.cokrige.rdata")
cokrige <- st_interp(cokrige, st_polygonize(border))
names(cokrige)
## [1] "x" "y" "logCd.pred" "logCd.var"
## [5] "logCu.pred" "logCu.var" "logPb.pred" "logPb.var"
## [9] "logZn.pred" "logZn.var" "cov.logCd.logCu" "cov.logCd.logPb"
## [13] "cov.logCu.logPb" "cov.logCd.logZn" "cov.logCu.logZn" "cov.logPb.logZn"
## [17] "geometry"
logCd.pred
、logCu.pred
等即为插值结果,logCd.var
、logCu.var
等为相应的插值结果的方差。
以锌为例对插值结果进行可视化:
ggplot(data = cokrige) +
geom_sf(aes(fill = exp(logZn.pred)), col = NA) +
geom_sf(data = border) +
theme_bw() +
theme(text = element_text(family = "mono")) +
scale_fill_continuous(low = "white", high = "red", name = "Zn")
![](https://i-blog.csdnimg.cn/blog_migrate/9bcf731190c855846e6b82ab3d5c661a.png)
6 交叉验证
可以使用交叉验证(Cross Validation,CV)对插值结果进行评价。常使用的k折交叉验证(K-fold cross-validation)将样本平均分为k份,其中k-1份作为训练样本,剩余1份作为验证样本,重复k次。
单变量克里金插值(即非协同克里金)的交叉验证函数为krige.cv
:
krige.cv(formula, locations, model = NULL, ..., beta = NULL,
nmax = Inf, nmin = 0, maxdist = Inf,
nfold = nrow(locations),
verbose = interactive(), debug.level = 0)
nfold:k参数;
verbose:控制是否显示进度条的参数。
v.fit <- vgm(0.59, "Sph", 874, 0.04)
cv <- krige.cv(log(zinc) ~ 1, data, v.fit,
nfold = 5, verbose = T)
cor(cv$var1.pred, cv$observed)
## [1] 0.8510259
协同克里金的交叉验证函数为gstat.cv
:
gstat.cv(object, nfold, remove.all = FALSE,
verbose = interactive(),
all.residuals = FALSE, ...)
gstat.cv
函数默认只对第一个属性变量进行交叉验证:
object:
gstat
函数的输出结果;all.residuals:若为TRUE,则生成所有变量的预测残差。
cv2 <- gstat.cv(g, verbose = T)
cor(cv2@data$logCd.pred, cv2@data$observed)
## [1] 0.5563733