Loading ...
Sorry, an error occurred while loading the content.

1006Re: AI-GEOSTATS: kriging estimates in S+

Expand Messages
  • Edzer J. Pebesma
    May 9 4:05 AM
    • 0 Attachment
      Vanessa,

      you are right. I was confused, as MASS has a free library
      called "spatial", but you meant the commercial add-on module
      "S+SpatialStats", which is also called "module spatial" by
      the vendor of S-Plus.

      It seems that S+SpatialStats, when predicting at point locations,
      does as if the locations are shifted a very small amount; I guess
      this is saying the same thing that Isobel meantioned earlier (see
      my sample below). When predicting without a nugget effect, the
      differences are also too large, IMO (see x2 example); I do not
      have an explanation for that.

      I also tried the examples with the (my) gstat S library, which can
      be downloaded from http://www.gstat.org/ , where both examples do
      what you'd expect; again see below.

      Thanks for sharing this with us.
      --
      Edzer

      Here's the output of the S-Plus session:

      S-PLUS : Copyright (c) 1988, 2002 Insightful Corp.
      S : Copyright Lucent Technologies, Inc.
      Version 6.1.2 Release 2 for Linux 2.2.12 : 2002
      Working data will be in .Data
      > module(spatial)
      > x<-as.numeric(1:10)
      > y<-as.numeric(1:10)
      > z<-as.numeric(1:10)
      > d<-data.frame(x=x,y=y,z=z)
      > d
      x y z
      1 1 1 1
      2 2 2 2
      3 3 3 3
      4 4 4 4
      5 5 5 5
      6 6 6 6
      7 7 7 7
      8 8 8 8
      9 9 9 9
      10 10 10 10
      > krige(z~loc(x,y),data=d,covfun=gauss.cov,range=3,sill=1,nugget=1)
      Call:
      krige(formula = z ~ loc(x, y), data = d, covfun = gauss.cov, range = 3, sill =
      1, nugget = 1)

      Coefficients:
      constant
      5.5

      Number of observations: 10
      > x<-krige(z~loc(x,y),data=d,covfun=gauss.cov,range=3,sill=1,nugget=1)
      > predict(x, d)
      x y fit se.fit
      1 1 1 2.827990 1.206120
      2 2 2 2.778034 1.155678
      3 3 3 3.432520 1.152567
      4 4 4 4.299386 1.153139
      5 5 5 5.113279 1.153538
      6 6 6 5.886721 1.153538
      7 7 7 6.700614 1.153139
      8 8 8 7.567480 1.152567
      9 9 9 8.221966 1.155678
      10 10 10 8.172010 1.206120
      > d2<-data.frame(x=y+1e-7,y=y+1e-7)
      > d2
      x y
      1 1 1
      2 2 2
      3 3 3
      4 4 4
      5 5 5
      6 6 6
      7 7 7
      8 8 8
      9 9 9
      10 10 10
      > predict(x, d2)
      x y fit se.fit
      1 1 1 2.827990 1.206120
      2 2 2 2.778034 1.155678
      3 3 3 3.432520 1.152567
      4 4 4 4.299386 1.153139
      5 5 5 5.113279 1.153538
      6 6 6 5.886721 1.153538
      7 7 7 6.700614 1.153139
      8 8 8 7.567480 1.152567
      9 9 9 8.221966 1.155678
      10 10 10 8.172010 1.206120
      > x2<-krige(z~loc(x,y),data=d,covfun=gauss.cov,range=3,sill=2,nugget=0)
      > predict(x2, d)
      x y fit se.fit
      1 1 1 1.000004 0.001272791
      2 2 2 1.999994 0.001272785
      3 3 3 3.000009 0.001272774
      4 4 4 3.999990 0.001272763
      5 5 5 5.000010 0.001272757
      6 6 6 5.999990 0.001272757
      7 7 7 7.000010 0.001272763
      8 8 8 7.999991 0.001272774
      9 9 9 9.000006 0.001272785
      10 10 10 9.999996 0.001272791
      >
      # Now using the gstat library:
      > library(gstat)
      > krige(z~1,~x+y,d,d,vgm(1,"Gau",3,1))
      [using ordinary kriging]
      x y var1.pred var1.var
      1 1 1 1 0.000000e+00
      2 2 2 2 0.000000e+00
      3 3 3 3 5.185705e-33
      4 4 4 4 2.074282e-32
      5 5 5 5 5.185705e-33
      6 6 6 6 2.074282e-32
      7 7 7 7 0.000000e+00
      8 8 8 8 5.185705e-33
      9 9 9 9 0.000000e+00
      10 10 10 10 0.000000e+00
      > krige(z~1,~x+y,d,d,vgm(2,"Gau",3,0))
      [using ordinary kriging]
      x y var1.pred var1.var
      1 1 1 1 7.106437e-33
      2 2 2 2 0.000000e+00
      3 3 3 3 4.440892e-16
      4 4 4 4 2.220446e-16
      5 5 5 5 2.220446e-16
      6 6 6 6 4.440892e-16
      7 7 7 7 2.220446e-16
      8 8 8 8 2.842575e-32
      9 9 9 9 2.842575e-32
      10 10 10 10 0.000000e+00
      >

      --
      * To post a message to the list, send it to ai-geostats@...
      * As a general service to the users, please remember to post a summary of any useful responses to your questions.
      * To unsubscribe, send an email to majordomo@... with no subject and "unsubscribe ai-geostats" followed by "end" on the next line in the message body. DO NOT SEND Subscribe/Unsubscribe requests to the list
      * Support to the list is provided at http://www.ai-geostats.org
    • Show all 7 messages in this topic