Мой код вариограммы отличается от результата variog () - PullRequest
Новая
       73

Мой код вариограммы отличается от результата variog ()

0 голосов
/ 30 января 2019

Я пишу код для создания вариограммы.Для проверки своего результата я проверил с помощью geoR :: variog (), но обе вариограммы разные.

Я попытался понять код variog (), чтобы увидеть, что происходит под капотом, но происходит так много вещейчто я не могу понять это.Я в своем коде использую параметры X-координата, Y-координата, значение данных, количество лагов, минимальное значение лага, интервал отставания, азимут (угол в градусах; 90 соответствует вертикальному направлению), допуск угла (в градусах) и максимальная пропускная способность. 

variogram = function(xcor, ycor, data, nlag, minlag, laginv, azm, atol, maxbandw){
dl <- length(data)
lowangle <- azm - atol
upangle <- azm + atol
gamlag <- integer(nlag)
n <- integer(nlag)

dist <- pairdist(xcor, ycor)
maxd <- max(dist)
llag <- seq(minlag, minlag + (nlag-1) * laginv, by = laginv)
hlag <- llag + laginv

for(i in 1:dl){
    for(j in i:dl){
        if(i != j){

            if(xcor[j]- xcor[i] == 0)
                theta <- 90
            else
                theta <- 180/pi * atan((ycor[j] - ycor[i])/(xcor[j] - xcor[i]))
            for(k in 1:nlag){
                d <- dist[j, i]
                b <- abs(d * sin(theta - azm))
                if((llag[k] <= d & d < hlag[k]) & (lowangle <= theta & theta < upangle) & (b <= maxbandw)){
                    gamlag[k] <- gamlag[k] + (data[i] - data[j])^2;
                    n[k] <- n[k] + 1
                }
            }
        }
    }   
}

gamlag <- ifelse(n == 0, NA, gamlag/(2*n))
tmp <- data.frame("lag" = llag, "gamma" = gamlag)
return(tmp)
}

вызов функции для приведенного выше кода 

ideal_variogram_2 <- variogram(data3[,1], data3[,2], data3[,3], 18, 0, 0.025, 90, 45, 1000000)
ideal_variogram_2 <-  na.omit(ideal_variogram_2)
plot(ideal_variogram_2$lag, ideal_variogram_2$gamma, main = "Using my code")

вызов функции для variog () 

geodata1 <- as.geodata(data3, coords.col = 1:2, data.col = 3)
ideal_variogram_1 <- variog(geodata1, coords = geodata1$coords, data = geodata1$data, option = "bin", uvec = seq(0, 0.45, by = 0.025), direction = pi/2, tolerance = pi/4)
df <- data.frame(u = ideal_variogram_1$u, v = ideal_variogram_1$v)
plot(df$u, df$v, main = "Using variog()")

2 полученных мне вариограммынаходятся по следующей ссылке: Вариограмма

...