用gstat局部变异函数进行局部块克里格

Question

我一直无法找到任何特定于本地块克里格的信息，使用R中的gstat包进行局部变异函数，澳大利亚精确农业中心提供了一种名为VESPER的免费软件，能够做到这一点。据我所读到，在R中应该是可能的，我只需要一些帮助就可以将gstat函数放在本地工作。

以meuse数据集为例，我能够计算出一个全局变异函数并将其拟合到一个数据集中：

    library(gstat)
    data(meuse)
    coordinates(meuse) = ~x+y
    data(meuse.grid)
    gridded(meuse.grid) = ~x+y

    logzinc_vgm<- variogram(log(zinc)~1, meuse)
    logzinc_vgm_fit <- fit.variogram(logzinc_vgm, model=vgm("Sph", "Exp"))
    logzinc_vgm_fit

    plot(logzinc_vgm, logzinc_vgm_fit)

这给出了一个很好的变异图的整个数据集与拟合模型。然后，我可以使用它对整个数据集执行块克里格：

    logzinc_blkkrig <- krige(log(zinc)~1, meuse, meuse.grid, model = logzinc_vgm_fit, block=c(100,100))
    spplot(logzinc_blkkrig["var1.pred"], main = "ordinary kriging predictions")
    spplot(logzinc_blkkrig["var1.var"],  main = "ordinary kriging variance")

这产生了插值数据的图以及每个预测点的方差图。所以如果我想让这些函数为我的整个数据集工作一次，那就太完美了.

但是我一直无法生成一个for循环来处理本地级别的这些函数。

我的目标是: 1.对于我的网格文件中的每一点(我尝试作为数据框架和SpatialPointsDataFrame)，我想从我的数据文件点在全局变差图中给定的范围对角线上的距离内子集(即logzinc_vgm_fit2，3) 2。在这个数据子集上，我想计算变异图(如上)，并将模型拟合到它(如上面) 3。我希望执行块kriging来获得网格点4的预测值和方差。

注意:与内置在gstat包中的meuse数据集一样，我的网格和数据帧的尺寸是不同的。

如果有人能解决这个问题的话，非常感谢你的参与。很高兴发布我目前正在使用的代码，如果它有用的话。

Answer 1

我做了一个for循环，我认为它完成了您的要求。我不认为这是需要块克里格，因为循环预测在每个网格单元。

rad参数是搜索半径，可以设置为其他数量，但目前引用全局变异函数范围(具有金块效应)。我认为最好进一步搜索点，因为如果只搜索到全局变异函数范围，则局部方差函数拟合可能不会收敛(即没有观测范围)。

k参数用于rad中最近邻的最小数目。这一点很重要，因为某些位置可能在rad中没有点，这将导致错误。

您应该注意到，您指定的model=vgm("Sph", "Exp")方式似乎采用了第一个列出的方法。所以，我在for循环中使用了球形模型，但是您可以更改为您想要使用的。母亲可能是一个很好的选择，如果你认为形状会随着位置的变化。

#Specify the search radius for the local variogram
rad = logzinc_vgm_fit[2,3]
#Specify minimum number of points for prediction
k = 25
#Index to indicate if any result has been stored yet
stored = 0
for (i in 1:nrow(meuse.grid)){
  #Calculate the Euclidian distance to all points from the currect grid cell
  dists = spDistsN1(pts = meuse, pt = meuse.grid[i,], longlat = FALSE)

  #Find indices of the points within rad of this grid point
  IndsInRad = which(dists < rad)

  if (length(IndsInRad) < k){
    print('Not enough nearest neighbors')
  }else{
    #Calculate the local variogram with these points
    locVario = variogram(log(zinc)~1, meuse[IndsInRad,])

    #Fit the local variogram
    locVarioFit = fit.variogram(logzinc_vgm, model=vgm("Sph"))

    #Use kriging to predict at grid cell i. Supress printed output.
    loc_krig <- krige(log(zinc)~1, meuse[IndsInRad,], meuse.grid[i,], model = locVarioFit, debug.level = 0)

    #Add result to database
    if (stored == 0){
      FinalResult = loc_krig
      stored = 1
    }else{
      FinalResult = rbind(FinalResult, loc_krig)
    }
  }
}