## AI-GEOSTATS: gstat and anisotropy parameters.

Expand Messages
• 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
• 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]
• ... 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.