Универсальный кригинг R с autoKrige()

Я пытаюсь использовать функцию autoKrige() в пакете automap для простого применения универсального кригинга. У меня есть неравномерная сетка измерений, и я хочу интерполировать между ними в точном пространственном масштабе. Пример кода:

    library('automap')

    # create an irregularly spaced grid
    y <-x <-c(-5,-4,-2,-1,-0.5,0,0.5,1,2,4,5)
    grid <-expand.grid(x,y)
    names(grid) <-c('x', 'y')

    # create some measurements, greatest in the centre, with some noise
    vals <-apply(grid,1, function(x) {12/(0.1+sqrt(x[1]^2 + x[2]^2))+rnorm(1,2,1.5)})

    # get data into sp format
    s <-SpatialPointsDataFrame(grid, data.frame(vals))

    # make some prediction locations and get them into sp format
    pred <-expand.grid(seq(-5,5,by=0.5), seq(-5,5,by=0.5))
    pred <-cbind(pred[,1], pred[,2])    # this seems to be needed, not sure why
    pred <-SpatialPoints(pred)

    # try universal kriging
    surf <-autoKrige(vals~x+y, s, new_data=pred)    

Это приводит к ошибке:

    Error in gstat.formula.predict(d$formula, newdata, na.action = na.action,  : 
      NROW(locs) != NROW(X): this should not occur

Я пытался сделать так, чтобы в new_data было то же количество строк, что и в исходных данных, и даже пытался сделать координаты в new_data точно такими же, как в исходных данных, но я все еще получаю эту ошибку. Я новичок в геостатистических методах, так что извините, если я делаю основную ошибку. Может ли кто-нибудь посоветовать, где я ошибаюсь? Спасибо.


person Jonathan Denniss    schedule 07.09.2015    source источник
comment
Шаг pred <-cbind(pred[,1], pred[,2]) для меня можно опустить, autoKrige работает и без него.   -  person Paul Hiemstra    schedule 08.09.2015


Ответы (1)


Проблема в том, что у вас неправильный синтаксис функции autoKrige. Ввод формулы в autoKrige определяет линейную модель, которую вы хотите использовать, например:

log(zinc) ~ dist

из набора данных meuse. В этом случае вы моделируете log(zinc) по сравнению с dist с помощью линейной модели, а остатки этой модели интерполируются с помощью вариограммы. По сути, универсальный кригинг представляет собой линейную регрессию с пространственно коррелированными остатками.

В вашем случае вы указываете:

val ~ x+y

поэтому autoKrige (фактически gstat) попытается сначала смоделировать линейную модель vals по сравнению с x и y (многомерная регрессия) и интерполировать остатки, используя модель вариограммы. Однако переменные x и y отсутствуют в файле SpatialPointsDataFrame.

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

vals ~ 1

где определяется среднее значение vals, а остатки интерполируются с использованием модели вариограммы. На самом деле это известно как обычный кригинг. Ваш вызов autoKrige будет выглядеть примерно так:

surf <-autoKrige(vals ~ 1, s, new_data=pred) 
person Paul Hiemstra    schedule 08.09.2015