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

Question

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

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

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

Answer 1

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

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

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

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