【问题标题】:What KInd of variogram should I plot before doing Kriging in R?在 R 中进行克里金之前,我应该绘制哪种变异函数?
【发布时间】:2015-10-19 10:23:57
【问题描述】:

我有一个名为 seoul3112 的 csv 文件,其中包含 PM10 浓度。 please, download..我尝试绘制样本变异函数并在其上拟合模型。

library(sp)
    library(gstat)
    library(rgdal)
    seoul3112<-read.csv("seoul3112.csv",row.name=1)
    seoul3112<-na.omit(seoul3112)

#分配一个 CRS 并重新投影

coordinates(seoul3112)=~LON+LAT
proj4string(seoul3112) =  "+proj=longlat +datum=WGS84" 
seoul3112<-spTransform(seoul3112, CRS("+proj=utm +north +zone=52 +datum=WGS84"))

#plot 半变异函数

g<-gstat(id="PM10",formula=PM10~LON+LAT, data=seoul3112)
seoul3112.var<-variogram(g, cutoff=70000, width=6000)
seoul3112.var
plot(seoul3112.var, col="black", pch=16,cex=1.3,
     xlab="Distance",ylab="Semivariance",
     main="Omnidirectional Variogram for seoul 3112")

#模型拟合

model.3112<- fit.variogram(seoul3112.var,vgm(700,"Gau",40000,400),fit.method = 2)
                       
plot(seoul3112.var,model=model.3112, col="black", pch=16,cex=1.3,
     xlab="Distance",ylab="Semivariance",
     main="Omnidirectional Variogram for seoul 3112")

写完这些代码后,我得到了一个像这样的半变异函数。

由于我是地统计学的新手,所以我很困惑我上面的变异函数是否适合我的数据集。因为,在典型的变异函数中,半方差值在门槛处变成水平的。但是这个变异函数是向上的!我应该在我的代码中做一些更正吗?

另一件事是,实际上我的最终目标是对我的数据集(seoul3112)进行克里金插值。我不明白,对于进行克里金,这个样本变异函数就足够了,还是我应该绘制方向变异函数或其他任何东西? 谁能详细解释一下?

【问题讨论】:

  • 因为您的问题是关于选择和解释统计方法而不是关于编程,我相信stats.stackexchange.com 是一个更合适的站点。我已经投票决定转让它。

标签: r spatial kriging gstat


【解决方案1】:

如果您查看您的拟合模型,它将有一个sill 参数(即nugget + psill),但它超出了您的样本范围。

 sill = sum(model.3112$psill)

您可能没有足够远的点对距离到达rangesill。我认为这不是问题,只要您使用此半变异函数在数据的空间域 + 60000 m(您有数据支持的距离)内进行预测。我采取这个立场是因为要拟合的半变异函数最重要的部分是开始,并且在数据的范围内拟合是好的。

需要检查的是您有多少点对 (np) 支持您绘制的 bins(点对越多越好)。

另一个建议是使用水平图寻找各向异性

seoul3112.var_map<-variogram(g, cutoff=70000, width=6000, map=TRUE)
plot(seoul3112.var_map, col.regions=terrain.colors(20))

或通过设置alpha 并绘制模型拟合来查看多个方向的拟合。您只需要检查 180 度,因为它是对称的(您可以在下图中看到,其中 0 和 180 相同)。

seoul3112.var_angles<-variogram(g, cutoff=70000, width=6000, alpha=seq(0,180,30))
plot(seoul3112.var_angles, model=model.3112, pch=16, ylim=c(0,3000))

如果在方向上存在差异,您可以使用 fit.variogram 和为 vgm() 模型定义的 anis 对各向异性的长轴和短轴进行建模。 示例:

vgm(700,"Gau",40000,400, anis=c(someangle, someratio))

【讨论】:

  • 从以上7个方向变差函数如何理解是否存在各向异性?我知道几何和纬向各向异性。但是我可以从上面的数字中说,窗台在各个方向上几乎相同吗?而且我也无法理解这些数字中不同或不同的范围!您能否根据此图评论一下各向异性?
  • 您真正要寻找的是rangesill 的哪个方向最短,哪个方向最长。从水平图和方向变异函数看来,30 度可能是短轴(半方差有一些增加,可能是门槛),这将使 120 度成为长轴。我建议分别绘制这些角度并查看这些变异函数的样子。
  • 谢谢。但是我看到了一些示例,它们在表格中显示了每个方向的窗台和范围的数值。通过它我可以清楚地了解主轴和次轴。链接:techmat.vgtu.lt/~art/proc/file/BudrLi.pdf(第 365 页)。 R中是否有任何内置函数,我可以通过它查看每个方向的范围和窗台的数值?
  • 如果sills 在各个方向都到达,那么是的,您可以通过在为vgm() 定义的@987654346 上使用fit.variogram 来制作一个类似于您引用的表格@ 与感兴趣的alpha(在alpha 上有合理的tol.hor)。然后使用 sill = sum(model.3112$psill)range = model.3112$range 之类的东西来提取它们。
  • 你是这个意思吗? seoul3112.var1&lt;-variogram(g,width=8000,cutoff=80000, alpha=seq(0,180,45),tol.hor=15) // model&lt;- fit.variogram(seoul3112.var1,vgm(800,"Sph",40000,200),fit.method = 2) // sill = sum(model$psill) // range = model$range // sill // range 我试过了。但没有得到我想要的结果。 (这里。我在问题中创建的对象 g)
猜你喜欢
  • 2016-12-08
  • 2015-11-28
  • 2019-02-09
  • 2017-04-12
  • 2017-10-08
  • 1970-01-01
  • 1970-01-01
  • 2013-04-16
  • 1970-01-01
相关资源
最近更新 更多