1

我正在编写生成变异函数的代码。为了验证我的结果，我使用 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 个变异函数位于以下链接：变异 函数

4

0 回答 0