krige: Kriging

In sgeostat: An Object-Oriented Framework for Geostatistical Modeling in S+

Description

Carry out spatial prediction (or kriging).

Usage

krige(s, point.obj, at, var.mod.obj, maxdist=NULL, extrap=FALSE, border) 

Arguments

s	a point object, generated by point(), at which prediction is carried out
point.obj	a point object, generated by point(), containing the sample points and data
at	the variable, contained in point.obj, for which prediction will be carried out
var.mod.obj	variogram object
maxdist	an optional maximum distance. If entered, then only sample points (i.e, in point.obj) within maxdist of each prediction point will be used to do the prediction at that point. If not entered, then all n sample points will be used to make the prediction at each point.
extrap	logical, indicates if prediction outside the convex hull of data points should be done, default FALSE
border	optional polygon (list with two components x and y of same length) representing a (possibly non convex) region of interest to be used instead of the convex hull. Needs extrap=TRUE.

Value

A point object which is a copy of the s object with two new variables, zhat and sigma2hat, which are, repspectively, the predicted value and the kriging variance.

References

http://www.gis.iastate.edu/SGeoStat/homepage.html

See Also

est.variogram,fit.variogram

Examples

# a single point:
prdpnt <- point(data.frame(list(x=180000,y=331000)))
prdpnt <- krige(prdpnt, maas.point, 'zinc', maas.vmod)
prdpnt

# kriging on a grid (slow!)
grid <- list(x=seq(min(maas$x),max(maas$x),by=100),
             y=seq(min(maas$y),max(maas$y),by=100))
grid$xr <- range(grid$x)
grid$xs <- grid$xr[2] - grid$xr[1]
grid$yr <- range(grid$y)
grid$ys <- grid$yr[2] - grid$yr[1]
grid$max <- max(grid$xs, grid$ys)
grid$xy <- data.frame(cbind(c(matrix(grid$x, length(grid$x), length(grid$y))),
             c(matrix(grid$y, length(grid$x), length(grid$y), byrow=TRUE))))
colnames(grid$xy) <- c("x", "y")
grid$point <- point(grid$xy)
data(maas.bank)
grid$krige <- krige(grid$point,maas.point,'zinc',maas.vmod,
                    maxdist=1000,extrap=FALSE,border=maas.bank)
op <- par(no.readonly = TRUE)
par(pty="s")
plot(grid$xy, type="n", xlim=c(grid$xr[1], grid$xr[1]+grid$max),
                    ylim=c(grid$yr[1], grid$yr[1]+grid$max))
image(grid$x,grid$y,
      matrix(grid$krige$zhat,length(grid$x),length(grid$y)),
      add=TRUE)
contour(grid$x,grid$y,
        matrix(grid$krige$zhat,length(grid$x),length(grid$y)),
        add=TRUE)
data(maas.bank)
lines(maas.bank$x,maas.bank$y,col="blue")
par(op)

sgeostat documentation built on May 1, 2019, 6:31 p.m.