Hi all,
I am trying to figure out if I can get 3-D kriging to work in the R
gstat package. The demo given here:
# =================================================
# Edzer J. Pebesma, Richard N.M. Duin (2005) Spatio-temporal mapping of
# sea floor sediment pollution in the North Sea. In: Ph. Renard, and
# R. Froidevaux, eds. Proceedings GeoENV 2004 -- Fifth European Conference
# on Geostatistics for Environmental Applications; Springer.
#
# Run the demo:
demo(pcb)
#==================================================
...doesn't really do full 3D kriging as far as I can tell, it just
models the cross-variograms between data from different years. I would
eventually like to do a kriging prediction map for e.g. any one of, say,
100 or 1000 different years, so I don't think the cross-kriging approach
will work.
Anyway, for now, I am just seeing if I can get a simple 3D kriging to
work with the pcb dataset. In an attempt to make it work, I rescaled
the "year" data to the approximate dimensions of the x and y data, and
then I added 1% variability in all the locations to see if that would
avoid the singularity problem, which I gather (?) can be caused by
points too close together.
But still, no luck. I basically just want to fit a variogram model
which captures the variability in space (xy) and time (rescaled_year).
Here's what I've got, below, any comments/help VERY appreciated.
(PS: Does anyone have a bit of code that will calculate an empirical
variogram in the third dimension? The gstat variogram() function
evidently won't do it, even when beta=90 is specified.
Thanks!
Nick
====================
data(pcb)
# rescale year to similar units as space
(horiz_range = max(pcb$x) - min(pcb$x))
(vert_range = max(pcb$year) - min(pcb$year))
(range_ratio = horiz_range / vert_range)
pcb$rescaled_year = (pcb$year - min(pcb$year)) * range_ratio
# add a little noise to the data locations in case there are overlapping
points
pcb$x = as.double(pcb$x + runif(length(pcb$x), -0.01*horiz_range,
0.01*horiz_range))
pcb$y = as.double(pcb$y + runif(length(pcb$y), -0.01*horiz_range,
0.01*horiz_range))
pcb$rescaled_year = as.double(pcb$rescaled_year +
runif(length(pcb$rescaled_year), -0.01*vert_range, 0.01*vert_range))
# do a 2D variogram for various years
v3gm = NULL
v4gm = NULL
# get residuals after factoring out depth
pcb$res=residuals(lm(log(PCB138)~rescaled_year+depth, pcb))
# Get a variogram of the residuals by location, after factoring out any
year correlation
v3 = variogram(res ~ rescaled_year, ~x+y, pcb, dX=.1,
bound=c(0,1000,3000,5000,(1:16)*10000))
v3gm = vgm(.224,"Sph",17247,.08)
print(plot(v3, model = v3gm, plot.numbers = TRUE))
(v3gm.f <- fit.variogram(v3, v3gm, fit.ranges=FALSE))
print(plot(v3, model = v3gm.f, plot.numbers = TRUE))
print.data.frame(v3gm)
print.data.frame(v3gm.f)
# do the 3D variogram
v4 = variogram(res ~ 1, ~x+y+rescaled_year, pcb, dX=.5)
print(plot(v4, model = vgm(.224,"Exp",17247,.08), plot.numbers = TRUE))
# 3-dimensional model with nugget component and sill component
#
v4gm = vgm(0.3, "Sph", 2000, anis=c(0, 90, 0, 1, 1), add.to=v3gm)
print.data.frame(v3gm)
print.data.frame(v4gm)
(v4gm.f <- fit.variogram(v4, v4gm, fit.sills=TRUE, fit.ranges=TRUE))
print.data.frame(v4gm.f)
print(plot(v4, model = v4gm.f, plot.numbers = TRUE))
====================