1

我一直无法使用 R 中的 gstat 包使用本地变异函数找到任何特定于本地块克里金法的信息。澳大利亚精准农业中心有一个名为 VESPER 的免费软件能够做到这一点，并且从我读过的内容来看在 R 中应该是可能的，我可以使用一些帮助来组合一个 for 循环以使 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_fit[2,3]）） 2. 在这个数据子集上，我想计算变异函数（如上）并拟合模型（如上） 3. 基于这个模型，我想执行块克里金法以获得该网格点的预测值和方差 4. 将上述三个步骤构建为一个 for 循环，以根据每个网格点周围的局部变异函数预测每个网格点的值

注意：与 gstat 包中内置的 meuse 数据集一样，我的网格和数据数据框的维度是不同的

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

4

1 回答 1

1

我做了一个 for 循环，我认为它可以满足您的要求。我认为不需要块克里金法，因为循环预测每个网格单元。

参数为搜索半径，rad可设置为其他量，但目前引用全局变差函数范围（带金块效果）。我认为最好进一步搜索点，因为如果您只搜索全局变异函数范围，局部变异函数拟合可能不会收敛（即没有观察到的范围）。

该k参数用于 内的最近邻居的最小数量rad。这很重要，因为某些位置可能在 内没有点rad，这会导致错误。

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

#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)
    }
  }
}
于 2017-02-27T01:41:54.087 回答