variogram on 3D points with gstat ?
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