Я пытаюсь сделать универсальный кокригинг в R с пакетом Gstat. У меня есть сценарий, с которым мне помогли, но теперь я застрял и не могу обратиться за помощью к первоисточнику. Проблема в том, что я не могу изменить выходное разрешение данных 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:
border <- shapefile("area_2014.shp")
projection(kova)=projection(border)
Это должно создать сетку для кокригинга, и 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 и сценария в целом будет принята с благодарностью!
