我正在尝试使用 Gstat 包在 R 中进行通用协同克里金。我有一个脚本,我得到了帮助,但现在我被卡住了,无法从原始来源寻求帮助。问题是我无法更改 cokriged 数据的输出分辨率。我想将插值地图导入 ArcMap,而点对点栅格的分辨率非常低。
我的脚本如下:
library(raster)
library(gstat)
library(sp)
library(rgdal)
library(FitAR)
加载我的数据集,其中包含坐标和采样值:
kova<-read.table("katvus_point_modif3.txt",sep=" ",header=T)
coordinates(kova)=~POINT_X+POINT_Y
在与前一个相同的坐标处加载深度值,这是我的协变量:
Sygavus<-read.table("sygavus_point_cokrig.txt",sep=" ",header=T)
coordinates(Sygavus)=~POINT_X+POINT_Y
overlay <- over(kova,Sygavus)
kova$Sygavus <- overlay$Sygavus
这应该为我的插值设置边界,该文件是从 ArcMap 导出的 shapefile:
border <- shapefile("area_2014.shp")
projection(kova)=projection(border)
这应该为 cokriging 创建一个网格, res= 应该让我指定我希望输出的分辨率,但无论我使用什么数字,输出都不会改变。
grid <- spsample(border,type="regular",res=25)
我删除重叠点:
zero <- zerodist(kova)
kova <- kova[-zero[,2],]
我加载深度协变量光栅文件。这是从 ArcMap 到 ascii 格式的深度栅格导出:
depth <- raster("htp_depth_covar.asc")
projection(depth)=projection(border)
overlay <- extract(depth,kova)
kova$depth <- overlay
我去掉呐!来自叠加深度值的值(这些值应与先前在相应坐标处加载的深度协变量表相同,但如果我不考虑该部分,脚本将停止运行)
kova <- kova[!is.na(kova$depth),]
kova.gstat <- gstat(id="Kova",formula=kova~depth,data=kova)
kova.gstat <- gstat(kova.gstat,id="Sygavus",formula=Sygavus~depth,data=kova)
var.kova <- variogram(kova.gstat)
plot(var.kova)
kova.gstat <- gstat(kova.gstat,id=c("Kova","Sygavus"),model=vgm(psill=cov(kova$kova,kova$Sygavus),model="Mat",range=12000,nugget=0))
kova.gstat <- fit.lmc(var.kova,kova.gstat,model=vgm(psill=cov(kova$kova,kova$Sygavus),model="Mat",range=12000,nugget=0))
plot(var.kova,kova.gstat$model)
overlay <- extract(depth,grid)
grid <- as.data.frame(grid)
grid$depth <- overlay
coordinates(grid)=~x1+x2
projection(grid)=projection(border)
krige <- predict.gstat(kova.gstat,grid)
spplot(krige,c("Kova.pred"))
write.table(krige, "kova.raster1.ck.csv", sep=";", dec=",", row.names=F)
对于理解 gstat cokriging 和整体脚本的任何帮助将不胜感激!