| krige {sgeostat} | R Documentation |
Carry out spatial prediction (or kriging).
krige(s, point.obj, at, var.mod.obj, maxdist=NULL, extrap=FALSE, border)
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. |
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.
http://www.gis.iastate.edu/SGeoStat/homepage.html
# 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)