[R-sig-Geo] anisotropy modeling [SEC=UNCLASSIFIED]

Jon Olav Skoien jon.skoien at jrc.ec.europa.eu
Tue Aug 16 12:29:00 CEST 2011


Hi Jin,

some comments below:

On 16-Aug-11 3:38, Jin.Li at ga.gov.au wrote:
> Hi All,
>
> Following Edzer's examples, I have tried to test the effects of the specification of alpha in variogram on the results of fit.variogram as below:
>
>>   v =  variogram(zinc~1,meuse)
>>   v0 = variogram(zinc~1,meuse, alpha = 0, tol.hor = 45)
>>   v90 =  variogram(zinc~1,meuse, alpha = 90, tol.hor = 45)
>>   v090 =  variogram(zinc~1,meuse, alpha=c(0, 90), tol.hor=45)
>>   v45135 =  variogram(zinc~1,meuse, alpha=c(45, 135), tol.hor=45)
>>   v2 =  variogram(zinc~1,meuse, alpha=c(0, 45, 90, 135), tol.hor=45)
>
>>    fit.variogram(v, vgm(1, "Exp", 300,  anis=c(90,1)))
>    model  psill    range
> 1   Exp 166238 323.6674
>>    fit.variogram(v, vgm(1, "Exp", 300, anis=c(90,.1)))
>    model  psill    range ang1 anis1
> 1   Exp 166238 3236.674   90   0.1
>>    fit.variogram(v090, vgm(1, "Exp", 300, anis=c(90,.1)))
>    model    psill    range ang1 anis1
> 1   Exp 166262.4 3237.814   90   0.1
>>    fit.variogram(v45135, vgm(1, "Exp", 300, anis=c(90,.1)))
>    model    psill    range ang1 anis1
> 1   Exp 165919.4 3224.622   90   0.1
>>   fit.variogram(v2, vgm(1, "Exp", 300, anis=c(90,.1)))
>    model  psill    range ang1 anis1
> 1   Exp 166091 3231.219   90   0.1
>
>
>>    fit.variogram(v0, vgm(1, "Exp", 300, anis=c(90,.1)))
>    model    psill    range ang1 anis1
> 1   Exp 146068.7 2772.751   90   0.1
>>    fit.variogram(v90, vgm(1, "Exp", 300, anis=c(90,.1)))
>    model    psill    range ang1 anis1
> 1   Exp 200767.2 4027.282   90   0.1
>
> Above results show that using v, v090, v45135 and v2 would produce similar results for anis=c(90, 0.1), suggesting that the specification of alpha in variogram is largely redundant, while using v0 and v90 produced very different results, suggesting a contrary. Could anyone shed a light on this? Thanks.
> Jin
>

I would not say that the last results are very different, if you plot 
e.g. v2 with the three last models:
m1 = fit.variogram(v2, vgm(1, "Exp", 300, anis=c(90,.1)))
m2 = fit.variogram(v0, vgm(1, "Exp", 300, anis=c(90,.1)))
m3 = fit.variogram(v90, vgm(1, "Exp", 300, anis=c(90,.1)))
plot(v2, m1)
plot(v2, m2)
plot(v2, m3)

you can see that all of them give quite similar fits to the variogram. 
Small fitted range comes together with small fitted psill, which means 
that their shapes will be rather similar. The reason for v0 and v90 
giving different fitted variograms is probably that these just use a 
part of the data set (the angles 0 and 90), whereas the other variograms 
use the complete data set.

For anisotropy modelling, you can also use the 
estimateAnistropy-function in the intamap-package:
 > estimateAnisotropy(meuse, "zinc")
$ratio
[1] 1.618686

$direction
[1] 48.54651

$Q
            Q11       Q22        Q12
[1,] 0.9195639 0.8232523 -0.3870018

$doRotation
[1] TRUE

where the ratio is the inverse (1/ratio) of the one used in gstat. 
Through interpolate with methodName = "automap" (or with 
estimateParameters and spatialPredict if you want to do the steps 
separately), you get a variogram model fitted to the anisotropic sample 
variogram with direction and ratio from estimateAnisotropy.

Cheers,
Jon


>
> -----Original Message-----
> From: r-sig-geo-bounces at r-project.org [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of Edzer Pebesma
> Sent: Tuesday, 19 July 2011 5:59 PM
> To: shweta jp
> Cc: r-sig-geo at r-project.org
> Subject: Re: [R-sig-Geo] anisotropy modeling
>
> The description says:
>
>    The automatic fitting is done through fit.variogram.
>
> estimateAnisotropy {intamap} does fit anistropy coefficients.
>
>
>
> On 07/18/2011 08:27 PM, shweta jp wrote:
>> Is the same true for autofitVariogram {automap} ?
>>
>> Shweta
>>
>> On Wed, Jul 13, 2011 at 1:13 PM, Edzer Pebesma<
>> edzer.pebesma at uni-muenster.de>  wrote:
>>
>>> help(fit.variogram) has the following description:
>>>
>>>      Fit ranges and/or sills from a simple or nested variogram model to
>>>      a sample variogram
>>>
>>> I suspected it would completely ignore anisotropy settings. However,
>>> this seems not to be the case:
>>>
>>>> library(gstat)
>>> Loading required package: sp
>>> Loading required package: spacetime
>>> Loading required package: zoo
>>> Loading required package: xts
>>>> loadMeuse()
>>>> v0 = variogram(zinc~1,meuse, alpha = 0, tol.hor = 45)
>>>> v90 = variogram(zinc~1,meuse, alpha = 90, tol.hor = 45)
>>>> v = variogram(zinc~1,meuse)
>>>> fit.variogram(v, vgm(1, "Exp", 300, anis=c(90,1)))
>>>   model  psill    range
>>> 1   Exp 166238 323.6674
>>>> fit.variogram(v, vgm(1, "Exp", 300, anis=c(90,.1)))
>>>   model  psill    range ang1 anis1
>>> 1   Exp 166238 3236.674   90   0.1
>>>> fit.variogram(v0, vgm(1, "Exp", 300, anis=c(90,.1)))
>>>   model    psill    range ang1 anis1
>>> 1   Exp 146068.7 2772.751   90   0.1
>>>> fit.variogram(v90, vgm(1, "Exp", 300, anis=c(90,.1)))
>>>   model    psill    range ang1 anis1
>>> 1   Exp 200767.2 4027.282   90   0.1
>>>
>>> On 07/13/2011 11:34 AM, Matevž Pavlič wrote:
>>>> Hi,
>>>>
>>>> thanks for the reply. Can you explain this a little bit more :
>>>> " it assumes the variogram to fit to is in the major (correlation)
>>> direction (or averaged over all directions)."
>>>> I don't know exactly ehat it means....
>>>> m
>>>>
>>>> -----Original Message-----
>>>> From: r-sig-geo-bounces at r-project.org [mailto:
>>> r-sig-geo-bounces at r-project.org] On Behalf Of Edzer Pebesma
>>>> Sent: Wednesday, July 13, 2011 10:10 AM
>>>> To: r-sig-geo at r-project.org
>>>> Subject: Re: [R-sig-Geo] anisotropy modeling
>>>>
>>>>
>>>>
>>>> On 07/12/2011 10:52 PM, Matevž Pavlič wrote:
>>>>> Hi all,
>>>>>
>>>>>
>>>>>
>>>>> i 'm am not entirely sure I understand anisotropy kriging in R.
>>>>>
>>>>> I have a data set in which is (at least i think so) anisotropy is
>>> clearly visible. So I made directional variograms....
>>>>> So as i understand, the only thing that is different from>normal<
>>> kriging is the>anis<  property in which you define the diffrenece of
>>> ranges? Is that coreect?
>>>> Yes. Please note that fit.variogram ignores anything about anisotropy, it
>>> assumes the variogram to fit to is in the major (correlation) direction (or
>>> averaged over all directions).
>>>>> Bellow is the code i use :
>>>>>
>>>>>
>>>>>
>>>>> a<-5000/7000
>>>>>
>>>>> print(plot(variogram(Z ~ 1, DF, map = TRUE, cutoff = 15000, width =
>>>>> 100), main = "Variogram map, podlaga",col.regions =
>>>>> terrain.colors(64)))
>>>>>
>>>>> v1.a<-variogram(Z~1, DF, alpha=c(45, 135))
>>>>>
>>>>> (vmf.a<- fit.variogram(v1.a, vgm(2900, "Pen", 5500, 300, anis =
>>>>> c(135, a))))
>>>>>
>>>>> print(plot(v1.a, pl = F, pch = 20, col = "blue"))
>>>>>
>>>>> print(plot(v1.a, plot.numbers = F, pch = 20, col = "darkblue", model =
>>>>> vmf.a))
>>>>>
>>>>>
>>>>>
>>>>> podlaga.aniso<-krige(Z~1, DF[-zerodist(DF)[,1],], grd, vmf.a )
>>>>>
>>>>> print(spplot(podlaga.aniso, zcol="var1.pred",
>>>>> col.regions=terrain.colors(64), contour=T, pretty=T, cuts=15,
>>>>> key.space="right"))
>>>>>
>>>>> print(spplot(podlaga.aniso,  zcol="var1.var",
>>>>> col.regions=terrain.colors(64), contour=T, pretty=T, cuts=15,
>>>>> key.space="right"))
>>>>>
>>>>> writeGDAL(podlaga.aniso, "podlaga_aniso.tif")
>>>>>
>>>>>
>>>>>
>>>>> thanks for info,
>>>>>
>>>>>
>>>>>
>>>>> m
>>>>>
>>>>>
>>>>>       [[alternative HTML version deleted]]
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> R-sig-Geo at r-project.org
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>> --
>>>> Edzer Pebesma
>>>> Institute for Geoinformatics (ifgi), University of Münster Weseler Straße
>>> 253, 48151 Münster, Germany. Phone: +49 251 8333081, Fax: +49 251 8339763
>>> http://ifgi.uni-muenster.de
>>>> http://www.52north.org/geostatistics      e.pebesma at wwu.de
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo at r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>> --
>>> Edzer Pebesma
>>> Institute for Geoinformatics (ifgi), University of Münster
>>> Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
>>> 8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
>>> http://www.52north.org/geostatistics      e.pebesma at wwu.de
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>



More information about the R-sig-Geo mailing list