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

AI-GEOSTATS: gstat and anisotropy parameters.

Expand Messages
  • Nick Hamm
    Hi I have a set of two dimensional data( measured on a grid). Initial analysis strongly suggests that there is directional anisotropy, so I ve been trying to
    Message 1 of 3 , Sep 4, 2003
      Hi

      I have a set of two dimensional data( measured on a grid). Initial
      analysis strongly suggests that there is directional anisotropy, so I've
      been trying to deal with this.

      One way I've found is to fit the variogram model by maximum likelihood
      (or REML) using geoR. geoR seems to account for anisotropy by deforming
      the locations, according to the anisotropy parameters (direction of
      maximum range, and the anisotropy ratio). Cross--validation suggests
      that this performs well, relative to fitting various models to
      isotropic sample variogram.

      The thing is, I want to use gstat for kriging. The reason for this is
      that I want to do block kriging. However, gstat seems to account for
      anisotropy by adjusting the variogram range according to the anisotropy
      parameters.

      So it seems that geoR deforms the coordinates, whereas gstat
      "deforms" the variogram....

      The upshot of this is that, when I use gstat to do punctual kriging
      using the parameters estimated in geoR, I get different results to what
      I get when I use geoR for kriging. A possible alternative, is to use
      geoR to predict the deformed locations, and then use gstat to do the
      kriging WITHOUT specifying the anisotropy parameter. If I do this then
      gstat and geoR give the same results for kriging. This suggests that
      there is no difference in the algorithms used for kriging between the
      two packages. This approach seems to be OK for doing punctual
      kriging... However, I want to do block kriging (square blocks). The
      approach I've adopted for punctual kriging seems to fall over
      here. This is because I would specify the block size (in gstat) in the
      deformed coordinate space, whereas I actually want the blocks to be
      specified in the original geographic coordinate space. Am I making
      sense?

      What a mess!!! Can any one suggest an alternative, or is this all a bit
      silly?

      cheers

      Nick



      --
      * 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
    • Edzer J. Pebesma
      Nick, thanks for sharing this with us. The good news of course is that both gstat and geoR give the same results for isotropic kriging. The trouble is that
      Message 2 of 3 , Sep 5, 2003
        Nick, thanks for sharing this with us. The good news of course is that both
        gstat and geoR give the same results for isotropic kriging. The trouble is
        that software developers like Paulo and me are mainly bothered with
        developing our own thing, and not too much with compatibility to the
        others, although last March we agreed to make an effort for exactly this
        reason. A similar thing that came up a few years ago on ai-geostats was
        the specification of the parameters for the exponential and Gaussian
        variogram models -- see
        http://ai-geostats.org/Geostats_Faq/FAQ_software_conventions.htm

        and we need such a page for modelling (and fitting!) anisotropy
        parameters.

        I did a little experiment with two generated anisotropic
        random fields, x and y, R file and png images attached.
        The geoR output does estimate the ratio fairly well:
        a value of 0.3 for gstat corresponds to 3.333.. in geoR;
        however I do not understand the psiA value; geoR docs
        say:

        Anisotropy angle: defined here as the azimuth angle of the direction
        with greater spatial continuity, i.e. the angle between the y-axis and
        the direction with the maximum range.

        another difference is that gstat uses degrees (0-180),
        whereas I think geoR uses radial degrees (?), (0-pi).
        Still, in that case geoR should give values close to 0.11 (20/180)
        and 0.611 (110/180). Paulo?
        Also, the sigmasq values worry me a bit:
        > var(x[,3])
        [1] 0.8750844
        > var(y[,3])
        [1] 0.9214435


        generating variogram(x): 1sph(10, 20, 0.3);
        >
        likfit(as.geodata(x),ini=c(1,10),fix.nugget=T,cov.model="spherical",fix.psiA=F,
        psiA=0.111, fix.psiR=F, psiR=3)
        ....
        likfit: estimated model parameters:
        beta sigmasq phi psiA psiR
        "-0.0713" " 2.5342" " 9.5813" " 0.3248" " 3.7673"


        variogram(x): 1sph(10, 110, 0.3);
        >
        likfit(as.geodata(y),ini=c(1,10),fix.nugget=T,cov.model="spherical",fix.psiA=F,
        psiA=0.05*c(1:20)*pi, fix.psiR=F, psiR=3)
        likfit: searching for best initial value ... selected values:
        sigmasq phi tausq kappa lambda psiR psiA
        initial.value " 1.0" "10.0" " 0.0" " 0.5" " 1.0" " 3.0" " 2.0"
        status "est" "est" "fix" "fix" "fix" "est" "est"
        likelihood value: -539.58582544843
        ....
        beta sigmasq phi psiA psiR
        " 0.1186" " 2.3785" "10.0816" " 1.9426" " 3.1296"

        likfit: maximised log-likelihood = -337.4

        I hope we can resolve this soon.

        As a short-term answer: you can obtain block mean values
        by point kriging on a fine grid, and then averaging point kriged
        values within a block. It does not give you block kriging
        standard errors, though.
        --
        Edzer

        Nick Hamm wrote:

        >Hi
        >
        >I have a set of two dimensional data( measured on a grid). Initial
        >analysis strongly suggests that there is directional anisotropy, so I've
        >been trying to deal with this.
        >
        >One way I've found is to fit the variogram model by maximum likelihood
        >(or REML) using geoR. geoR seems to account for anisotropy by deforming
        >the locations, according to the anisotropy parameters (direction of
        >maximum range, and the anisotropy ratio). Cross--validation suggests
        >that this performs well, relative to fitting various models to
        >isotropic sample variogram.
        >
        >The thing is, I want to use gstat for kriging. The reason for this is
        >that I want to do block kriging. However, gstat seems to account for
        >anisotropy by adjusting the variogram range according to the anisotropy
        >parameters.
        >
        >So it seems that geoR deforms the coordinates, whereas gstat
        >"deforms" the variogram....
        >
        >The upshot of this is that, when I use gstat to do punctual kriging
        >using the parameters estimated in geoR, I get different results to what
        >I get when I use geoR for kriging. A possible alternative, is to use
        >geoR to predict the deformed locations, and then use gstat to do the
        >kriging WITHOUT specifying the anisotropy parameter. If I do this then
        >gstat and geoR give the same results for kriging. This suggests that
        >there is no difference in the algorithms used for kriging between the
        >two packages. This approach seems to be OK for doing punctual
        >kriging... However, I want to do block kriging (square blocks). The
        >approach I've adopted for punctual kriging seems to fall over
        >here. This is because I would specify the block size (in gstat) in the
        >deformed coordinate space, whereas I actually want the blocks to be
        >specified in the original geographic coordinate space. Am I making
        >sense?
        >
        >What a mess!!! Can any one suggest an alternative, or is this all a bit
        >silly?
        >
        >cheers
        >
        >Nick
        >
        >
        >
        >

        ----------

        "x" <-
        structure(list(x = c(-9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5,
        -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5,
        9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5,
        0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5,
        -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5,
        3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5,
        -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5,
        7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5,
        -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
        -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5,
        1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5,
        -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5,
        4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5,
        -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5,
        8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5,
        -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5,
        -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5,
        2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5,
        -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5,
        5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5,
        -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5,
        8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5,
        -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5,
        -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5,
        2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5,
        -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5,
        5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5,
        -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5,
        8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5,
        -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5,
        -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5,
        2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5,
        -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5,
        5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5,
        -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5,
        8.5, 9.5), y = c(9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5,
        9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 8.5, 8.5,
        8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5,
        8.5, 8.5, 8.5, 8.5, 8.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5,
        7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 6.5,
        6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5,
        6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5,
        5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5,
        4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5,
        4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5,
        3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5,
        3.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5,
        2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 1.5, 1.5, 1.5, 1.5, 1.5,
        1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5,
        1.5, 1.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
        0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, -0.5, -0.5, -0.5,
        -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
        -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -1.5, -1.5, -1.5, -1.5, -1.5,
        -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5,
        -1.5, -1.5, -1.5, -1.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5,
        -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5,
        -2.5, -2.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5,
        -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5,
        -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5,
        -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -5.5, -5.5,
        -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5,
        -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -6.5, -6.5, -6.5, -6.5,
        -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5,
        -6.5, -6.5, -6.5, -6.5, -6.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5,
        -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5,
        -7.5, -7.5, -7.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5,
        -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5,
        -8.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5,
        -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5),
        attr = c(-0.976396, 0.699642, 1.02706, 1.8274, -0.0842941,
        0.730785, 1.17241, 0.17138, -0.499801, 0.0290917, 0.971627,
        0.159505, 0.0978222, -0.0628971, 1.10981, 1.60993, 0.157566,
        -1.31656, -0.814423, -1.08538, -0.664614, 1.3288, 1.7715,
        0.648596, 0.512488, 0.125008, 1.15626, -0.293581, -1.40186,
        0.345669, 1.36306, 0.828358, -0.342381, 0.00905169, 1.24502,
        0.576963, -0.0811922, -1.44835, -0.319535, -0.341521, 0.515442,
        1.15643, 1.48897, 0.418224, 0.129235, 0.0275096, -0.505079,
        -1.1466, -0.922303, 0.501884, 1.3434, 0.690577, -0.0654742,
        0.240356, 0.833001, 0.465874, -1.25548, -0.841794, 0.228543,
        0.924724, 0.461841, 2.01067, 0.684979, -0.194078, -0.311548,
        0.73425, -0.682227, -1.32422, -1.31838, 0.762056, 1.92055,
        0.664635, 0.169215, 1.13904, -0.0118303, -0.42341, -2.09248,
        -0.443041, -0.521656, 1.20924, 1.60895, 0.541853, -1.32658,
        -0.34973, 0.320017, 0.0198023, -1.85797, -0.837399, -0.861881,
        0.752807, 1.19442, 0.55206, -0.337885, -0.268511, -0.0354799,
        -0.577635, -1.1632, -0.0418927, -0.190198, 1.38387, 0.864981,
        -0.230297, -0.908201, -0.317284, 0.540593, 0.588873, -2.15811,
        -1.29086, 0.564227, 0.45145, 0.00254214, 0.223103, -0.72224,
        -0.691557, 0.293142, -0.229129, -0.0911508, -0.446708, 0.0914658,
        1.91125, -0.0655018, -0.409656, -0.785249, 0.0342057, 1.0731,
        0.226422, -1.59694, -1.20873, 0.506341, 0.979477, 0.164993,
        -0.676357, -0.608324, -0.989886, -0.74416, 0.844582, -0.828441,
        0.627724, 1.75291, 1.97964, -0.147929, -0.929661, -1.73665,
        -0.349646, 1.48683, -0.289314, -2.0301, -0.418626, 0.404249,
        1.43853, 0.0772116, -1.2291, -2.07862, -1.02916, -0.3847,
        0.233113, -0.484702, 1.29267, 2.53204, 1.27968, -0.598517,
        -0.369162, -0.592886, 0.937435, 1.3464, -0.48181, -0.404962,
        0.930906, 1.7464, 1.18944, -0.177596, -1.65704, -1.86563,
        -0.59391, 0.396574, 0.285294, 0.150976, 1.0885, 3.13262,
        1.42033, -1.07894, -0.422835, 0.340313, 1.04254, 1.02172,
        -0.930621, -0.130838, 0.891712, 1.83336, 0.903686, -0.121988,
        -1.51042, -0.680568, 0.262939, 0.345589, -0.0345113, 1.29514,
        1.91254, 1.89058, 1.40523, -0.352911, -0.471958, 1.17492,
        1.06129, 0.0261414, -0.00657254, 0.270249, 0.831697, 1.82552,
        0.873856, -0.684648, -1.86114, -0.119203, 0.311865, 0.167739,
        -0.0225983, 1.37158, 1.93348, 0.282233, 1.81257, -0.0534109,
        -0.178551, 0.559363, 0.706015, -0.252069, -0.175302, -0.720996,
        -0.166321, 0.743316, 0.812269, -1.17577, -0.932122, -0.566561,
        0.718219, 0.188124, 0.735218, 1.6418, 0.178471, 0.407502,
        1.36858, -0.000146809, 0.291759, -0.468605, 0.0206043, -0.944232,
        -0.61307, 0.442481, -0.968002, 1.81035, 0.181125, -1.23142,
        -0.763612, 0.0737858, 1.12902, -0.254041, 0.927865, 0.886834,
        -0.571706, -0.13527, 1.27622, -0.907889, -0.657191, -0.758956,
        0.561848, -0.681267, -0.222691, 0.463942, 0.224039, 1.25724,
        0.25348, -1.35065, -0.709735, 1.73603, 0.84702, -0.958807,
        -0.203773, -0.749048, -0.215812, 0.948336, 0.8805, -0.953084,
        -1.49239, -0.0665188, -0.533943, -0.0828005, -0.870736, 0.65511,
        -0.137514, 0.552696, 0.0139825, -1.84643, -0.360674, 2.06584,
        0.0382007, -0.198137, -0.187399, -1.5025, -0.560836, 0.790284,
        0.477532, -0.86596, -0.68644, -0.690822, 0.627588, -0.635185,
        0.685255, 0.0380479, 0.681636, 0.503216, -0.0772367, -0.858596,
        0.890215, 0.857633, -0.705095, -0.759286, -1.27039, -0.652695,
        -0.123117, 1.34745, -0.0751269, -0.418802, -0.0673295, -0.265125,
        1.11624, -1.15422, 0.19013, -0.219525, 0.53884, -0.0382443,
        -0.948078, -0.477328, -0.634133, -0.369468, -1.50007, -1.08964,
        -1.31225, 0.270723, 0.442437, 1.9753, 0.434132, -1.16969,
        -0.198799, 0.606902, 0.711261, -0.719003, 0.473304, -0.649792,
        0.472475, -0.314922, -1.41962, -0.428929, -0.653978, -0.777673,
        -1.13312, -1.23852, -0.990636, 0.0964782, 0.26112, 1.37736,
        0.977292, 0.0433025, 0.878182, 1.06444, 0.89754, -0.130647,
        -0.538042, 0.0716454, 0.191813, -0.884949, -2.00202, -0.0712418,
        -1.35596, -0.373052, -1.25609, -1.68115, -1.31557, -0.605708,
        0.638957, 0.841501, 0.701605, -0.322548, -0.374385, 0.692801,
        0.963631, -0.0116963, -0.122099, 0.125823, -1.04663, -0.769083,
        -1.01459, -0.524525, -1.90085, -0.898233, -1.45772, -1.68847,
        -1.82645, -0.429825, 0.0917008, 0.948507, 0.305031)), .Names = c("x",
        "y", "attr"), class = "data.frame", row.names = c("1", "2", "3",
        "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15",
        "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26",
        "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37",
        "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "48",
        "49", "50", "51", "52", "53", "54", "55", "56", "57", "58", "59",
        "60", "61", "62", "63", "64", "65", "66", "67", "68", "69", "70",
        "71", "72", "73", "74", "75", "76", "77", "78", "79", "80", "81",
        "82", "83", "84", "85", "86", "87", "88", "89", "90", "91", "92",
        "93", "94", "95", "96", "97", "98", "99", "100", "101", "102",
        "103", "104", "105", "106", "107", "108", "109", "110", "111",
        "112", "113", "114", "115", "116", "117", "118", "119", "120",
        "121", "122", "123", "124", "125", "126", "127", "128", "129",
        "130", "131", "132", "133", "134", "135", "136", "137", "138",
        "139", "140", "141", "142", "143", "144", "145", "146", "147",
        "148", "149", "150", "151", "152", "153", "154", "155", "156",
        "157", "158", "159", "160", "161", "162", "163", "164", "165",
        "166", "167", "168", "169", "170", "171", "172", "173", "174",
        "175", "176", "177", "178", "179", "180", "181", "182", "183",
        "184", "185", "186", "187", "188", "189", "190", "191", "192",
        "193", "194", "195", "196", "197", "198", "199", "200", "201",
        "202", "203", "204", "205", "206", "207", "208", "209", "210",
        "211", "212", "213", "214", "215", "216", "217", "218", "219",
        "220", "221", "222", "223", "224", "225", "226", "227", "228",
        "229", "230", "231", "232", "233", "234", "235", "236", "237",
        "238", "239", "240", "241", "242", "243", "244", "245", "246",
        "247", "248", "249", "250", "251", "252", "253", "254", "255",
        "256", "257", "258", "259", "260", "261", "262", "263", "264",
        "265", "266", "267", "268", "269", "270", "271", "272", "273",
        "274", "275", "276", "277", "278", "279", "280", "281", "282",
        "283", "284", "285", "286", "287", "288", "289", "290", "291",
        "292", "293", "294", "295", "296", "297", "298", "299", "300",
        "301", "302", "303", "304", "305", "306", "307", "308", "309",
        "310", "311", "312", "313", "314", "315", "316", "317", "318",
        "319", "320", "321", "322", "323", "324", "325", "326", "327",
        "328", "329", "330", "331", "332", "333", "334", "335", "336",
        "337", "338", "339", "340", "341", "342", "343", "344", "345",
        "346", "347", "348", "349", "350", "351", "352", "353", "354",
        "355", "356", "357", "358", "359", "360", "361", "362", "363",
        "364", "365", "366", "367", "368", "369", "370", "371", "372",
        "373", "374", "375", "376", "377", "378", "379", "380", "381",
        "382", "383", "384", "385", "386", "387", "388", "389", "390",
        "391", "392", "393", "394", "395", "396", "397", "398", "399",
        "400"))
        "y" <-
        structure(list(V1 = c(-9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5,
        -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5,
        9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5,
        0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5,
        -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5,
        3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5,
        -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5,
        7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5,
        -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
        -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5,
        1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5,
        -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5,
        4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5,
        -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5,
        8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5,
        -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5,
        -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5,
        2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5,
        -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5,
        5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5,
        -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5,
        8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5,
        -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5,
        -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5,
        2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5,
        -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5,
        5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5,
        -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5,
        8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5,
        -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5,
        -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5,
        2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5,
        -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5,
        5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5,
        -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5,
        8.5, 9.5), V2 = c(9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5,
        9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 8.5, 8.5,
        8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5,
        8.5, 8.5, 8.5, 8.5, 8.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5,
        7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 6.5,
        6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5,
        6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5,
        5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5,
        4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5,
        4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5,
        3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5,
        3.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5,
        2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 1.5, 1.5, 1.5, 1.5, 1.5,
        1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5,
        1.5, 1.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
        0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, -0.5, -0.5, -0.5,
        -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
        -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -1.5, -1.5, -1.5, -1.5, -1.5,
        -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5,
        -1.5, -1.5, -1.5, -1.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5,
        -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5,
        -2.5, -2.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5,
        -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5,
        -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5,
        -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -5.5, -5.5,
        -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5,
        -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -6.5, -6.5, -6.5, -6.5,
        -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5,
        -6.5, -6.5, -6.5, -6.5, -6.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5,
        -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5,
        -7.5, -7.5, -7.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5,
        -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5,
        -8.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5,
        -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5),
        V3 = c(-0.358227, -0.883175, 0.428496, -0.213165, -0.121203,
        0.504585, 1.29542, 0.644304, 0.350139, 0.114496, 0.174196,
        0.328178, -0.0868328, 0.226939, 0.346281, 0.401713, 1.30887,
        1.76147, 0.853352, 0.121347, -0.337769, -0.953408, -1.69177,
        -1.00548, -0.837798, -0.422826, -0.0688711, 0.323329, 0.995217,
        0.638197, 0.762255, 0.870313, 0.711054, 1.04473, 1.34499,
        1.17326, 0.218925, 0.936686, 0.861577, 0.383784, 0.938214,
        1.08439, -0.0492151, -0.251227, -1.49285, -1.03058, -2.03424,
        -1.24326, -0.634366, -0.193798, 0.390262, 1.66108, 1.92824,
        1.5671, 0.191612, 0.338176, 1.50281, 1.0535, 0.551203, 0.0413718,
        0.377746, 0.723457, 0.142913, -0.89396, -0.782161, -0.981936,
        -0.973828, -1.30591, -1.73934, -1.84269, -0.0832526, 1.38839,
        1.51343, 2.0416, 2.24402, 1.90511, 1.97194, 0.324144, -0.23269,
        -0.692648, 1.05821, 1.03786, 0.784226, -0.339956, -0.676398,
        -1.01696, -0.95513, -0.487205, -1.22426, -1.56292, -1.04767,
        -0.265129, 0.585775, 1.53192, 1.06778, 2.4384, 1.06568, 1.29809,
        0.924595, -0.606835, 1.10346, 1.56817, 1.02077, 1.24847,
        0.248581, 0.224182, -0.717934, -0.188762, -0.3741, -1.05474,
        -0.597772, -0.564637, -0.921164, -0.271499, -0.906596, 0.492014,
        0.630646, 0.854272, 1.96541, 0.999921, 0.496202, 0.351973,
        1.60986, 2.33434, 1.51661, 0.167276, -0.146389, -0.225696,
        -0.892759, -0.0899573, -1.39838, -0.981012, -0.58096, -0.694041,
        -1.8745, -1.89785, -1.75994, -1.4497, 0.201812, 0.619787,
        -0.485159, 0.273018, 1.14061, 1.97519, 2.11462, 0.839004,
        0.717381, 0.824517, 0.259545, 0.817707, -1.27184, -0.457396,
        0.338396, -0.228584, 0.194588, -0.077346, -0.926435, -1.50319,
        -1.55432, -1.17051, -1.09185, -0.917399, -0.23146, -0.375039,
        0.155894, 0.220507, 0.231177, 0.97191, -0.0836007, -1.33339,
        -0.622234, -0.201292, -0.854076, -0.367267, 0.384303, 1.13566,
        1.03534, 0.516219, -0.48998, -1.18473, 0.196995, 0.262778,
        -0.0359566, 0.47255, -0.0905125, 0.357794, -0.285613, -0.00798318,
        -0.11457, -0.761042, -0.439379, -0.835278, -0.639411, -0.785074,
        -0.00523192, 0.1666, 1.05992, 1.06818, 0.833633, 0.385857,
        0.921412, 0.276106, 0.463932, 0.358063, -0.0233177, -0.238753,
        -0.42071, -0.251407, -0.282355, -0.0486227, -0.130755, -0.595462,
        0.490839, 0.339338, -0.0560451, 0.489945, 0.517167, 0.40683,
        1.00145, 0.589293, 1.53035, 0.337886, 0.536092, 0.748265,
        0.985178, 0.618256, -0.744364, -1.40118, -0.225276, -0.133966,
        -0.741826, -0.662154, -1.13363, -0.199749, 0.218788, -0.381617,
        0.579232, 1.01703, 0.765987, 0.745737, 1.26144, 0.722743,
        0.232609, 0.296584, 0.867954, 2.10431, 1.52535, 1.06925,
        0.314842, 0.356067, 0.0942154, -0.108979, -0.988771, -1.90368,
        -1.91997, -1.34157, -0.217298, 1.22345, 1.15753, 0.441481,
        1.55639, 0.882917, -0.792265, -0.0559937, 0.201641, 1.13765,
        1.73539, 1.28896, 0.854913, 0.901594, -0.0555984, -0.7484,
        -0.195218, -0.602409, -1.71683, -2.07703, -1.95228, -0.548744,
        0.0371549, 0.595531, 1.21788, 1.17501, -0.147261, -0.510002,
        -1.14883, -0.700317, -0.134978, 0.268485, 0.327767, 0.1389,
        0.243102, -0.289124, 0.309648, -0.290261, -0.564847, -0.641719,
        -1.75654, -2.05613, -1.75529, -0.0565981, 1.88894, 1.48677,
        0.646488, 0.0356698, -0.113009, -0.45765, -0.720445, -0.786172,
        -0.724996, 0.0962278, 0.0526066, -2.20774, -1.21379, -0.291403,
        -0.665068, -0.912585, -0.64879, -1.8832, -1.095, -0.90682,
        0.646028, 2.11661, 2.45733, 1.26387, 1.07263, -0.416093,
        -0.0959732, -0.327427, -0.167697, 0.422089, 0.144705, -0.401572,
        -0.348742, -0.650393, -0.136084, -1.01573, -0.546468, -1.49028,
        -0.893513, -0.51415, 0.0848532, 1.01496, 2.22971, 2.30453,
        2.29134, 0.976251, 0.176245, -0.561158, -0.930977, -0.628719,
        -0.172885, 0.518523, -0.825708, -1.04278, -0.0333577, -1.53102,
        -0.990126, -0.942557, -1.55075, -0.831715, 0.198945, 0.367643,
        1.02887, 1.86658, 0.936428, 1.05558, 0.426798, 0.176225,
        0.113876, -0.509315, -0.00407159, -0.0213874, 0.444478, 0.545452,
        -0.359192, -0.562048, -0.0992335, -0.869162, -1.66234, -1.9709,
        -0.449777, -0.263866, -0.270827, 0.194169, 1.67132, 1.44206,
        1.94036, 0.553896, -0.106475, 0.304114, 0.0289134, 0.0534254,
        -0.0438635, -0.445235, -0.455341, 0.202635, -0.232922, -0.798149,
        0.756377, -0.150466)), .Names = c("V1", "V2", "V3"), class = "data.frame", row.names = c("1",
        "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13",
        "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24",
        "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35",
        "36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46",
        "47", "48", "49", "50", "51", "52", "53", "54", "55", "56", "57",
        "58", "59", "60", "61", "62", "63", "64", "65", "66", "67", "68",
        "69", "70", "71", "72", "73", "74", "75", "76", "77", "78", "79",
        "80", "81", "82", "83", "84", "85", "86", "87", "88", "89", "90",
        "91", "92", "93", "94", "95", "96", "97", "98", "99", "100",
        "101", "102", "103", "104", "105", "106", "107", "108", "109",
        "110", "111", "112", "113", "114", "115", "116", "117", "118",
        "119", "120", "121", "122", "123", "124", "125", "126", "127",
        "128", "129", "130", "131", "132", "133", "134", "135", "136",
        "137", "138", "139", "140", "141", "142", "143", "144", "145",
        "146", "147", "148", "149", "150", "151", "152", "153", "154",
        "155", "156", "157", "158", "159", "160", "161", "162", "163",
        "164", "165", "166", "167", "168", "169", "170", "171", "172",
        "173", "174", "175", "176", "177", "178", "179", "180", "181",
        "182", "183", "184", "185", "186", "187", "188", "189", "190",
        "191", "192", "193", "194", "195", "196", "197", "198", "199",
        "200", "201", "202", "203", "204", "205", "206", "207", "208",
        "209", "210", "211", "212", "213", "214", "215", "216", "217",
        "218", "219", "220", "221", "222", "223", "224", "225", "226",
        "227", "228", "229", "230", "231", "232", "233", "234", "235",
        "236", "237", "238", "239", "240", "241", "242", "243", "244",
        "245", "246", "247", "248", "249", "250", "251", "252", "253",
        "254", "255", "256", "257", "258", "259", "260", "261", "262",
        "263", "264", "265", "266", "267", "268", "269", "270", "271",
        "272", "273", "274", "275", "276", "277", "278", "279", "280",
        "281", "282", "283", "284", "285", "286", "287", "288", "289",
        "290", "291", "292", "293", "294", "295", "296", "297", "298",
        "299", "300", "301", "302", "303", "304", "305", "306", "307",
        "308", "309", "310", "311", "312", "313", "314", "315", "316",
        "317", "318", "319", "320", "321", "322", "323", "324", "325",
        "326", "327", "328", "329", "330", "331", "332", "333", "334",
        "335", "336", "337", "338", "339", "340", "341", "342", "343",
        "344", "345", "346", "347", "348", "349", "350", "351", "352",
        "353", "354", "355", "356", "357", "358", "359", "360", "361",
        "362", "363", "364", "365", "366", "367", "368", "369", "370",
        "371", "372", "373", "374", "375", "376", "377", "378", "379",
        "380", "381", "382", "383", "384", "385", "386", "387", "388",
        "389", "390", "391", "392", "393", "394", "395", "396", "397",
        "398", "399", "400"))


        [Non-text portions of this message have been removed]
      • Nick Hamm
        ... This doesn t seem to work..... For another data set (no anisotropy) I have defined a 5 by 5 block and discretized it on a fine grid (400 points). Using
        Message 3 of 3 , Sep 11, 2003
          On Fri, 5 Sep 2003, Paulo Justiniano Ribeiro Jr wrote:

          > Edzer Pebesma wrote:
          >
          > > As a short-term answer: you can obtain block mean values
          > > by point kriging on a fine grid, and then averaging point kriged
          > > values within a block. It does not give you block kriging
          > > standard errors, though.
          > > --
          >
          > I think we can:
          > If we generate simulations on the fine grid, for every simulation compute
          > the mean at specified block (block kriging estimate) and then the variance
          > of those means you will have the block kriging std. errors



          This doesn't seem to work.....

          For another data set (no anisotropy) I have defined a 5 by 5 block and
          discretized it on a fine grid (400 points). Using gstat I use SGS to
          simulate 1000 realisations of each point. for each realisation I
          compute the mean of the 400 points (to give the value of the simulation
          on the block). I then have 1000 simulations of the mean for the block
          -- I comupte the meand and variance of these and get:

          > "Mean: " 8.748847
          > "Var : " 0.000704

          If I do block kriging straight off I get:

          > "Mean: " 8.746475 (ie block kriging prediction)
          > "Var: " 0.001332 (ie block kriging variance)

          If I do 1000 simulations of the block (ie SGS for the block) I get:

          > "Mean: " 8.748569
          > "Var: " 0.001390

          So the mean and variance from the simulations over the block seems to be
          basically the same as the block kriging prediction and
          variance. Using Paulo's approach the Mean looks ok, but the Variance
          seems to be underestimated...... I've checked my code and it looks
          ok...

          Nick



          --
          * 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
        Your message has been successfully submitted and would be delivered to recipients shortly.