本篇是介绍克里金插值的第三篇推文，也是最后一篇。
因为前面两篇使用的数据中已知点的样本太少，本篇使用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对象的方法了：
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)
加载相关工具包：
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)
拟合协同半变异函数的方法如下：
vgm<-vgm(psill=1,model="Sph",range=800,nugget=1)fit<-fit.lmc(vg,g,vgm)plot(vg,fit)
插值过程使用predict函数：
cokrige<-predict(fit,grid)save(cokrige,file="32-1.cokrige.rdata")##LinearModelofCoregionalizationfound.Good.##[usingordinarycokriging]
插值结果：
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")
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
如果觉得《gstat | 空间插值（四）——克里金插值之协同克里金和交叉验证》对你有帮助，请点赞、收藏，并留下你的观点哦！
