[R-sig-Geo] variogram on 3D points with gstat ?

[Edzer J. Pebesma]

Matthew, try

library(gstat)
demo(gstat3D)

It doesn't use the sp classes, and in the next version it will be like this:

# $Id: gstat3D.R,v 1.4 2006-02-10 19:05:02 edzer Exp $
# simple demo of 3D interpolation of 50 points with random normal values,
# randomly located in the unit cube
n <- 50

data3D <- data.frame(x = runif(n), y = runif(n), z = runif(n), v = rnorm(n))
coordinates(data3D) = ~x+y+z

range1D <- seq(from = 0, to = 1, length = 20)
grid3D <- expand.grid(x = range1D, y = range1D, z = range1D)
gridded(grid3D) = ~x+y+z

res3D <- krige(formula = v ~ 1, data3D, grid3D, model = vgm(1, "Exp", .2))

library(lattice)

levelplot(var1.pred ~ x + y | z, as.data.frame(res3D))

--
Edzer

Matthew Perry wrote:
> Any one know if kriging data in 3 dimensions is possible with R/gstat? Any
> examples?
>
> Thanks,
> Matt
>
>
> On 2/11/07, Matthew Perry <perrygeo at gmail.com> wrote:
>   
>> Hey folks,
>>
>>  I'm trying to use R/gstat to krige data in 3 dimensions. Is this
>> possible? I haven't had much luck...
>>
>>     
>>> library(gstat)
>>>       
>> Loading required package: sp
>>     
>>> d <- read.table("C:\\Workspace\\temp\\points.txt",sep="\t",header=TRUE)
>>> coordinates(d) = ~x+y+z
>>> summary(d)
>>>       
>> Object of class SpatialPointsDataFrame
>> Coordinates:
>>   min max
>> x   1   3
>> y   1   3
>> z   1   3
>> Is projected: NA
>> proj4string : [NA]
>> Number of points: 27
>> Data attributes:
>>       conc
>>  Min.   :   2.0
>>  1st Qu.:   5.0
>>  Median :  20.0
>>  Mean   : 653.2
>>  3rd Qu.:1000.0
>>  Max.   :2000.0
>>     
>>> d
>>>       
>>    coordinates conc
>> 1    (1, 1, 1)  500
>> 2    (1, 1, 2) 1000
>> 3    (1, 1, 3) 2000
>> 4    (2, 1, 1)  500
>> 5    (2, 1, 2) 1000
>> 6    (2, 1, 3) 2000
>> 7    (3, 1, 1)  500
>> 8    (3, 1, 2) 1000
>> 9    (3, 1, 3) 2000
>> 10   (1, 2, 1)    5
>> 11   (1, 2, 2) 1000
>> 12   (1, 2, 3) 2000
>> 13   (2, 2, 1)    5
>> 14   (2, 2, 2)   20
>> 15   (2, 2, 3) 2000
>> 16   (3, 2, 1)    5
>> 17   (3, 2, 2)   20
>> 18   (3, 2, 3) 2000
>> 19   (1, 3, 1)    2
>> 20   (1, 3, 2)    5
>> 21   (1, 3, 3)   20
>> 22   (2, 3, 1)    2
>> 23   (2, 3, 2)    5
>> 24   (2, 3, 3)   20
>> 25   (3, 3, 1)    2
>> 26   (3, 3, 2)    5
>> 27   (3, 3, 3)   20
>>     
>>> cvgm <- variogram(log(conc)~1,data=d)
>>>       
>> Error in as.vector(x, mode) : invalid argument 'mode'
>>
>>
>> When I set "coordinates(d) = ~x+y", the variogram works fine. Does the R
>> interface to gstat allow for 3D data?
>>
>> --
>> Matthew T. Perry
>> http://www.perrygeo.net
>>     
>
>
>
>
>