我正在嘗試在R中使用Gstat軟件包進行通用剋剋裏格處理。我有一個腳本,我被幫助,但現在我卡住了,不能要求原始來源的援助。 問題是我無法更改cokriged數據的輸出分辨率。我想將插入的地圖導入ArcMap,點對柵格以非常低的分辨率離開我。R:Gstat通用剋剋裏格分辨率
我的腳本如下:
library(raster)
library(gstat)
library(sp)
library(rgdal)
library(FitAR)
加載我的數據集,即containes座標和採樣值:
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中導出的Shape文件:
border <- shapefile("area_2014.shp")
projection(kova)=projection(border)
據說這是爲了創造協同克里金法網格和RES =應該讓我指定我所要的輸出是,但不管我用什麼號碼什麼的分辨率輸出不會改變。
grid <- spsample(border,type="regular",res=25)
我刪除overlaping點:
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
我刪除了na!從疊加式的深度值的值(這些值應該是相同的,在各自座標上以前加載的深度協表,但是如果我離開的那部分,腳本停止運行)
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和腳本整體將不勝感激!