krige {sgeostat} | R Documentation |
Kriging
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.obj |
a point object, generated by |
at |
the variable, contained in |
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 |
border |
optional polygon (list with two components |
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
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)