3D variogram model

---------- Forwarded message ----------
From: *Edzer Pebesma* <edzer.pebesma at uni-muenster.de 
<mailto:edzer.pebesma at uni-muenster.de>>
Date: Sun, Oct 25, 2009 at 3:15 PM
Subject: Re: [R-sig-Geo] 3D variogram model
To: matzke at berkeley.edu <mailto:matzke at berkeley.edu>
Cc: r-sig-geo at stat.math.ethz.ch <mailto:r-sig-geo at stat.math.ethz.ch>, 
gstat-info at geo.uu.nl <mailto:gstat-info at geo.uu.nl>


Nick, thanks for bringing this up! I uploaded a script and the (large)
resulting pdf to

http://ifgi.uni-muenster.de/~epebe_01/pcb/ 
<http://ifgi.uni-muenster.de/%7Eepebe_01/pcb/>

I used it for a workshop about a year ago, just before geoenv 2008. It's
supposed to more or less do what you had in mind. The variography in
continuous space-time is hard to infer for this particular data set, I
think, and what I did below is rather quick and dirty, and builds upon
the earlier 2004 work which was already heavy on the side of
assumptions. I'll add this script to the demo section of the next gstat
R package.

A few more comments to your script are inline, below.

Nick Matzke wrote:

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

after your:
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

coordinates(pcb)=~x+y+rescaled_year
variogram(PCB138~1,pcb,beta=90,tol.ver=1) # or any small value

However you won't see much, as this data set doesn't time-replicated
locations.

 >
 > Thanks!
 > Nick
 >
 > ====================
 >
 > # 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))

This wouldn't be needed for 3D kriging.

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

You don't factor year correlation out here; you just exclude (dX=.1)
point pairs from different years. So it's spatial correlation only that
ends up in this v3 variogram.

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

This variogram computes distances in 3D, which is correct provided that
you took care (and knew) the appropriate xy vs t anisotropy in advance.
I don't think the dX makes sense in this case.

 >
 >
 > # 3-dimensional model with nugget component and sill component
 >
 > #
 > v4gm = vgm(0.3, "Sph", 2000, anis=c(0, 90, 0, 1, 1), add.to 

<http://add.to>=v3gm)
as both anisotropy ratios are 1, this model is isotropic after all, and
the 0,90,0 could be anything.

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

--
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of M?nster
Weseler Stra?e 253, 48151 M?nster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/
http://www.springer.com/978-0-387-78170-9 e.pebesma at wwu.de 
<mailto:e.pebesma at wwu.de>