首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >R: Gstat通用cokriging决议

R: Gstat通用cokriging决议
EN

Stack Overflow用户
提问于 2015-03-25 07:21:11
回答 1查看 595关注 0票数 1

我试着用Gstat软件包做R中的通用cokriging。我有一个脚本,我得到了帮助,但现在我被困住了,无法从原来的来源寻求帮助。问题是,我不能改变焦化数据的输出分辨率。我想导入插值的地图到ArcMap和点对点光栅留给我一个非常低的分辨率。

我的脚本如下:

代码语言:javascript
复制
  library(raster)
library(gstat)
library(sp)
library(rgdal)
library(FitAR)

加载数据集,其中包含坐标和采样值:

代码语言:javascript
复制
kova<-read.table("katvus_point_modif3.txt",sep="    ",header=T)
coordinates(kova)=~POINT_X+POINT_Y

在与前面相同的坐标下加载深度值,这是我的协变量:

代码语言:javascript
复制
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:

代码语言:javascript
复制
border <- shapefile("area_2014.shp")
projection(kova)=projection(border)

这应该是为cokriging创建一个网格,res=应该允许我指定我希望输出的分辨率,,但是不管我使用哪个数字,输出都不会改变

代码语言:javascript
复制
grid <- spsample(border,type="regular",res=25)

我移除重叠点:

代码语言:javascript
复制
zero <- zerodist(kova)
kova <- kova[-zero[,2],]              

我加载深度协变量光栅文件。这是一个从ArcMap到ascii表单的深度光栅导出:

代码语言:javascript
复制
depth <- raster("htp_depth_covar.asc")
projection(depth)=projection(border)

overlay <- extract(depth,kova)
kova$depth <- overlay

我移除娜!来自覆盖深度值的值(--这些值应该与先前在各自坐标下加载的深度协变量表相同,但如果我省略了该部分,脚本将停止运行)

代码语言:javascript
复制
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和脚本的整体将是非常感谢的!

EN

回答 1

Stack Overflow用户

发布于 2015-03-25 09:06:59

因为您没有提供一个可重复的示例,我只能猜测,但是我认为spsample忽略了res=25参数。尝试n=1000,然后增加这个值以获得更高的分辨率。

票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/29249681

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档