[R-sig-Geo] 3D variogram model

[Nick Matzke]

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



-- 
====================================================
Nicholas J. Matzke
Ph.D. Student, Graduate Student Researcher
Huelsenbeck Lab
Center for Theoretical Evolutionary Genomics
4151 VLSB (Valley Life Sciences Building)
Department of Integrative Biology
University of California, Berkeley

Lab websites:
http://ib.berkeley.edu/people/lab_detail.php?lab=54
http://fisher.berkeley.edu/cteg/hlab.html
Dept. personal page: 
http://ib.berkeley.edu/people/students/person_detail.php?person=370
Lab personal page: http://fisher.berkeley.edu/cteg/members/matzke.html
Lab phone: 510-643-6299
Dept. fax: 510-643-6264
Cell phone: 510-301-0179
Email: matzke at berkeley.edu

Mailing address:
Department of Integrative Biology
3060 VLSB #3140
Berkeley, CA 94720-3140

-----------------------------------------------------
"[W]hen people thought the earth was flat, they were wrong. When people 
thought the earth was spherical, they were wrong. But if you think that 
thinking the earth is spherical is just as wrong as thinking the earth 
is flat, then your view is wronger than both of them put together."

Isaac Asimov (1989). "The Relativity of Wrong." The Skeptical Inquirer, 
14(1), 35-44. Fall 1989.
http://chem.tufts.edu/AnswersInScience/RelativityofWrong.htm