Covariance Parameters for Krig in geoR ksline

376 views Asked by At

I have a small data set of locations and benzene concentrations in 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

and I am trying to krig the values in a grid (the property these wells reside on) like so

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))

this is when I start the Krig with

q <- ksline(c, cov.model="exp", cov.pars=c(10,3.33), nugget=0, locations=in_mat)

however, when looking at the output of this with

cbind(q$predict[1:10], q$krige.var[1:10])

i see

         [,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

these values do not change for the first 5000 rows... (cant view more because max.print = 5000... not sure how to change this either but that is a tangent..)

I am realizing that my

cov.pars = c(10,3.33)

being range and sill, are probably the issue.

the geoR.pdf, pg 19 describes what is expected from cov.pars however I am not sure how I should decide what these covariance parameters need to be.

Is there a method to find the appropriate values from my existing data or can I set these to generic values where my output will be similar to a kriging performed in the spatial analyst package of ESRI's ArcGIS?

ZR

::::EDIT:::

my geodata object was improperly converted... here is the correct way to do this

c <- as.geodata(b, coords.col = 2:3, data.col = 4, )

also...for the variogram,

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
1

There are 1 answers

2
Jesse Anderson On

One way to do this is to use the variofit function (also in the geoR package) to estimate the covariance parameters. For example, using your data and initial values:

vario <- variog(c)  # See other options here for binning, etc
# Note that the order of the cov.pars is variance, then range, (see your question)
fitted_model <- variofit(vario=vario, ini.cov.pars=c(10, 3.33), 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)

It is worth your time to look at the variogram, by the way.