Локальный блочный кригинг с локальной вариограммой с помощью gstat

Question

Мне не удалось найти какую-либо информацию, относящуюся к локальному блочному кригингу, с использованием локальной вариограммы с использованием пакета gstat в R. Есть бесплатная программа под названием 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 для прогнозирования значений в каждой точке сетки на основе локальной вариограммы вокруг каждой точки сетки

примечание: как и в случае набора данных meuse, встроенного в пакет gstat, размеры моей сетки и фреймов данных отличаются

Большое спасибо за помощь, если кто-то может решить этот вопрос. Рад опубликовать код, с которым я работаю до сих пор, если это будет полезно.

Answer 1

Я сделал цикл for, который, я думаю, выполняет то, что вы просите. Я не думаю, что для этого необходим блочный кригинг, потому что цикл предсказывает в каждой ячейке сетки.

rad Параметр - это радиус поиска, который может быть установлен на другие величины, но в настоящее время ссылается на глобальный диапазон вариограммы (с эффектом самородка). Я думаю, что было бы лучше поискать немного дальше для точек, потому что, если вы будете искать только до глобального диапазона вариограммы, локальное соответствие вариограммы может не сходиться (т.е. нет наблюдаемого диапазона).

k параметр для минимального количества ближайших соседей в пределах rad, Это важно, потому что некоторые места могут не иметь точек внутри rad, что приведет к ошибке.

Вы должны отметить, что так, как вы указали model=vgm("Sph", "Exp") кажется, взять первый из перечисленных методов. Итак, я использовал сферическую модель в цикле for, но вы можете перейти к тому, что вы хотите использовать. 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)
    }
  }
}