[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