0

I have a series of samples, which I wish to construct a variogram model of, and eventual kriging model. First, I neeed to detrend the data, as shown below:

samples
     x    y         z
1   180 1180  2.763441
2   180  240 -2.000000
3   380 1840  1.720087
4   720   80  4.056754
5   860  800  4.361503
6  620 1360  4.737717
7   980 1920  4.352956
8  1680  260  4.568255
9  1520  800  5.025272
10 1100 1220  4.693432
11  800 1460  2.470927
12  360 1900  1.455169
13  700 1760  2.894159
14  720 1540  2.115742
15  660 1480  1.749017
16  540 1680  3.291592
17  260 1280  2.962401
18  440 1640  2.422442
19  280 1260  2.966076
20  580 1580  3.178913
21  600 1220  3.752786
22  240 1700  1.748011
23  480 1440  3.106302
24  740 1880  4.827699
25  760 1320  3.603621
26 1560 1640  5.410076
27 1960 1980  6.145778
28 1520 1620  5.499064
29 1900 1820  5.316121
30 1780 1580  5.318344
31  100  740  2.019103
32  180  760  2.353693
33  140  200  1.714856
34  380  720  3.526107
35  240  580  3.075283
36  260  600  3.329397
37  340  360  3.188613
38  280  680  2.626241
39  420  700  3.211163
40  500  240  2.960805
41  460  280  3.171664
42  480  300  2.828883
43  400  640  3.227938
44  440  480  2.420358
45  300  560  4.021187
46 1380  220  5.364264
47 1500  740  5.344526
48 1240  380  4.632060
49 1420  360  4.012537
50 1280  800  4.122139
51 1400  600  5.033020
52 1300  640  4.215308
53 1460  200  5.116025
54 1220  440  4.550290
55 1200  520  3.788613
56 1540  340  5.772432
57 1520  660  5.656598
58 1480  260  5.423685
59 1360  780  4.728220
60 1260  240  3.683696


print(mean(samples$z))

h <- gstat(formula=z~1, locations=~x+y, data=samples)
samples.vgm <- variogram(h) 
plot(samples.vgm,main='Variogram of Samples NOT detrended') 


z = samples$z
x = samples$x
y = samples$y 
trend <- lm(z~x+y)

c = trend$coefficients[[1]]
a = trend$coefficients[[2]]
b = trend$coefficients[[3]]



Xs <- c()  
Ys <- c()  
Zs <- c()  

print('started the loop')
for (i in 1:nrow(samples)){
  i = samples[i,]
  x=i$x
  y=i$y
  z=i$z
  z_prime = z - (a*x+b*y+c)
  Xs <- c(Xs,x)
  Ys <- c(Ys,y)
  Zs <- c(Zs,z_prime) 
}
sampled <- data.frame(Xs,Ys,Zs)
print(sampled)
print('the length of sampled is')
print(length(sampled[[1]]))
print(levelplot(Zs~Xs+Ys,sampled))

x <- seq(0,2000,by=20)
y <- seq(0,2000,by=20)

pred.grid <- data.frame(x=rep(x,times=length(y)),y=rep(y,each=length(x)))

g <- gstat(formula=Zs~1, locations=~Xs+Ys, data=sampled)
sampled.vgm <- variogram(g) 
plot(sampled.vgm,main='Variogram of Samples hopefully detrended')

The problem is that the plot of detrended variogram (i.e. the variogram g above) looks exactly the same as the variogram h also above, which is NOT detrended. any reason why this happens??

The data is clearly different, The mean of the values in the detrended version is 0, as expected, but the non-detrended version the mean is around 3.556, also as expected.

Is there something I'm not catching here?

makansij
  • 9,303
  • 37
  • 105
  • 183
  • From which package comes `gstat`? Also, `object 'g' not found` in `samples.vgm <- variogram(g)` (third line). –  Apr 03 '15 at 02:03
  • Perhaps this should have been mentioned earlier: when I run the following **immediately** after the above code, it errors: `vg.sph <- vgm(psill=1.0,model='Sph', range = 500) fit.sph <- fit.variogram(sampled.vgm, model = vg.sph) sk1 <- krige(formula=Zs~1, locations=~Xs+Ys, data=sampled, newdata=pred.grid, model=fit.sph, beta=0)` And the error message is `Error in [.data.frame (object, , -coord.numbers, drop = FALSE) : undefined columns selected` – makansij Apr 06 '15 at 23:22

1 Answers1

0

Not sure this question belongs here, since I think the issue is conceptual, not related to your code. I'm new, though, so I'll just go ahead and give you some quick feedback.

The variogram plots the variance (or semi-variance technically) of your data within a given spatial lag. When you apply a linear transformation to your data, I don't believe you'll alter the variance, and so you shouldn't see different parameters come out of the variogram model. Instead, your kriged surface just takes on a different mean value.

p.s. It would be helpful to make the code something that anyone can copy and paste -- e.g., include coded test data.

Allison
  • 11
  • 2
  • My understanding is that it does alter the variogram. Here's my logic: A stationary variable would have a bounded variogram because the variance is well defined, i.e. bounded at some maximum. A non-stationary variable would not have a bounded variogram because there would be a geospatial trend in the variable. Thus, to transform a non-stationary variable to a stationary variable, we must detrend, and to detrend we must linearly transform. – makansij Apr 06 '15 at 21:37