局部塊克里格與gstat本地變差函數

Question

我一直無法找到任何信息，具體到本地塊克里格與使用gstat包河還有就是免費從澳大利亞中心精準農業所謂VESPER，能夠局部變差要做到這一點，從我讀過的應該可以在R中，我可以使用一些幫助，把一個for循環，使gstat函數在本地工作。

使用默茲數據集作爲一個例子，我已經能夠計算和適應全球變差函數的數據的集合：

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包中的默認數據集一樣，我的網格和數據數據幀的尺寸也不相同

非常感謝您提供任何能夠解決此問題的信息。很高興發佈我目前正在使用的代碼，如果它會有用。

Answer 1

我做了一個循環，我認爲完成了你的要求。我不認爲塊克里格是必需的，因爲循環預測每個網格單元。

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) 
    } 
    } 
}