我有一个小的位置数据集和以 mg/kg 为单位的苯浓度
WELL.ID X Y BENZENE
1 MW-02 268.8155 282.83 0.00150
2 IW-06 271.6961 377.01 0.00050
3 IW-07 251.0236 300.41 0.01040
4 IW-08 278.9238 300.37 0.03190
5 MW-10 281.4008 414.15 2.04000
6 MW-12 391.3973 449.40 0.01350
7 MW-13 309.5307 335.55 0.01940
8 MW-15 372.8967 370.04 0.01620
9 MW-17 250.0000 428.04 0.01900
10 MW-24 424.4025 295.69 0.00780
11 MW-28 419.3205 250.00 0.00100
12 MW-29 352.9197 277.27 0.00031
13 MW-31 309.3174 370.92 0.17900
我正在尝试像这样对网格中的值(这些井所在的属性)进行计算
setwd("C:/.....")
getwd()
require(geoR)
require(ggplot2)
a <- read.table("krigbenz_loc.csv", sep = ",", header = TRUE)
b <- data.matrix(a)
c <- as.geodata(b)
x.range <- as.integer(range(a[,2]))
y.range <- as.integer(range(a[,3]))
x = seq(from=x.range[1], to=x.range[2], by=1)
y = seq(from=y.range[1], to=y.range[2], by=1)
length(x)
length(y)
xv <- rep(x,length(y))
yv <- rep(y, each=length(x))
in_mat <- as.matrix(cbind(xv, yv))
这是我开始克里格的时候
q <- ksline(c, cov.model="exp", cov.pars=c(10,3.33), nugget=0, locations=in_mat)
但是,当查看这个的输出时
cbind(q$predict[1:10], q$krige.var[1:10])
我明白
[,1] [,2]
[1,] 343.8958 10.91698
[2,] 343.8958 10.91698
[3,] 343.8958 10.91698
[4,] 343.8958 10.91698
[5,] 343.8958 10.91698
[6,] 343.8958 10.91698
[7,] 343.8958 10.91698
[8,] 343.8958 10.91698
[9,] 343.8958 10.91698
[10,] 343.8958 10.91698
前 5000 行的这些值不会改变...(无法查看更多信息,因为 max.print = 5000... 也不知道如何更改,但这是一个切线..)
我意识到我的
cov.pars = c(10,3.33)
范围和门槛,可能是问题所在。
geoR.pdf ,第 19 页描述了cov.pars 的预期,但是我不确定我应该如何决定这些协方差参数需要是什么。
有没有一种方法可以从我现有的数据中找到适当的值,或者我可以将它们设置为通用值,我的输出将类似于在 ESRI 的 ArcGIS 的空间分析包中执行的克里金法?
ZR
::::编辑:::
我的地理数据对象转换不正确......这是正确的方法
c <- as.geodata(b, coords.col = 2:3, data.col = 4, )
也......对于变异函数,
v1 <- variog(c)
length(v1$n)
v1.summary <- cbind(c(1:11), v1$v, v1$n)
colnames(v1.summary) <- c("lag", "semi-variance", "# of pairs")
v1.summary