[R-sig-Geo] Fwd: 3D variogram model
Nick Matzke
matzke at berkeley.edu
Tue Oct 27 05:58:25 CET 2009
Dr. Pebesma,
Thanks so much for your help!! I understand the pcb dataset is not a
great one for this, but it is extremely useful to see an example of the
basic procedure implemented. Thanks again!
Nick
Nick Matzke wrote:
>
>
> ---------- 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>
>
>
--
====================================================
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
More information about the R-sig-Geo
mailing list