同时试图从克里格苯值
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, coords.col = 2:3, data.col = 4, )
ggplot(a, aes(x= X, y= Y, colour="green", label=WELL.ID)) + geom_point() + geom_text(aes(label=WELL.ID),hjust=0, vjust=0)
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))
看看变异函数..(不是很漂亮,但正在努力)
### variogram ###
## on geo-object
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
plot(v1, type = "b", main = "Variogram: BENZENE at CRAIG BP")
使用 ksline 生成克里格值...
## variance of benzene readings = sd^2
sd <- sd(a$BENZENE)
var = sd^2
fitted_model <- variofit(vario=v1, ini.cov.pars=c(var, .29), cov.model='exp')
q <- ksline(c, cov.model=fitted_model$cov.model, cov.pars=fitted_model$cov.pars,
nugget=fitted_model$nugget, locations=in_mat)
但后来它拿着手机,当我尝试对结果进行成像时出错!!!!
> image(q, val = q$predict)
Error in eval(x$call$geodata, envir = attr(x, "parent.env"))$borders :
object of type 'builtin' is not subsettable
这似乎完全超出了左侧字段,因为我已经多次检查过......我用谷歌搜索了错误,似乎我正在尝试调用函数的子集,而 90% 的答案是我的语法某处是错误的,但我已经检查了所有内容,但我无法弄清楚......任何帮助将不胜感激。
谢谢
ZR