[52N Geostatistics] Gstat package. Problem with 3d anisotropic variogram

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

[52N Geostatistics] Gstat package. Problem with 3d anisotropic variogram

Julien Viarre
Dear all,

Sorry for the disturbance but I have a problem calculating anisotropic 3D-variograms with gstat package.

I simulated a rainfall squall line propagating from West to East with a velocity of 1km/min (figure 1 attached).
position.x represents longitude values
position.y represents latitude values 
position.t represents rescaled timesteps (time is considered as the third dimension)
For each time step, only one longitude band receives rainfall (other position have no rain) (positive rainfall values don't change)

Then I calculated anisotropic variogram in 3 dimensions, specifically with alpha = 0° (meridian direction) and alpha = 90° (zonal direction) ; beta = 0° (no velocity) and beta = 45° (velocity of the squall line simulated) 

vario = variogram(estim.d[,4]~1location= ~position.x+position.y+position.tdata = estim.d , alpha=angle.p[kp] , beta=angle.q[kq] , tol.hor=0 , tol.ver=0boundaries=classe_vario)

For alpha = 90° and beta = 45°, estimations should be correlated so I should expect a variogram with gamma=0 for all the classes...
Nevertheless, I obtained : 

    np       dist        gamma dir.hor dir.ver   id

1  30150  0.07071068  0.5773396664      90      45 var1

2  26550  0.14142136  0.5196323354      90      45 var1

3  47400  0.24584641  0.5358185764      90      45 var1

4   7550          0.35355339  0.5773182290      90      45 var1

5  32150  0.44746773  0.4290722330      90      45 var1

6  16300  0.56568542  0.4675393932      90      45 var1

7  16650  0.63639610  0.3920274991      90      45 var1

8  20050  0.75806784  0.3445372217      90      45 var1

9  12500  0.84852814  0.2323875687      90      45 var1

10 14050       0.93936997  0.0779285188      90      45 var1

11  8450  1.06066017  0.0003281226      90      45 var1

12  5000  1.13137085 0.0003727332      90      45 var1

13 10850       1.23336367  0.0004393653      90      45 var1

14  3150  1.34350288  0.0005170953      90      45 var1

15  3600  1.45644355  0.0005944342      90      45 var1

16  1800  1.55563492 0.0006728895      90      45 var1

17  1950  1.65172892  0.0007428695      90      45 var1

18   300          1.76776695  0.0007899247      90      45 var1

19   100          1.83847763 0.0007689443      90      45 var1

20    50          1.90918831 0.0005528209      90      45 var1


That is very surprising because for the first classes we have positive semivariance values.


And the most surprising is when I removed zero values to keep only the meridian positive band (figure 2 attached), I obtained consistent variogram estimation : 



Variogram associated : 

  

   np       dist gamma dir.hor dir.ver   id

1   600 0.07071068     0      90      45 var1

2   525 0.14142136     0      90      45 var1

3  1050 0.24917096     0      90      45 var1

4   150 0.35355339     0      90      45 var1

5   725 0.44620876     0      90      45 var1

6   350 0.56568542     0      90      45 var1

7   450 0.63639610     0      90      45 var1

8   575 0.75937120     0      90      45 var1

9   400 0.84852814     0      90      45 var1

10  475 0.93784689     0      90      45 var1

11  325 1.06066017     0      90      45 var1

12  200 1.13137085     0      90      45 var1

13  525 1.23575328     0      90      45 var1

14  175 1.34350288     0      90      45 var1

15  275 1.45921127     0      90      45 var1

16  150 1.55563492     0      90      45 var1

17  200 1.65286210     0      90      45 var1

18   50 1.76776695     0      90      45 var1

19   25 1.83847763     0      90      45 var1

20   25 1.90918831     0      90      45 var1




Any idea?

Thanks in advance,

With best regards,

Julien Viarre


_______________________________________________
Geostatistics mailing list
[hidden email]
http://list.52north.org/mailman/listinfo/geostatistics
http://geostatistics.forum.52north.org

figure1.png (349K) Download Attachment
figure2.png (126K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: [52N Geostatistics] Gstat package. Problem with 3d anisotropic variogram

Edzer Pebesma
Administrator
Julien, thanks for you detailed email. I hope to get back to you soon on
this issue, but it will take some time. Putting my brains in the 3D
anisotropy "mode" has always been quite an exercise, in the past.

In the meantime: you are aware of the errors made in the original GSLIB,
which were copied into gstat?

On 10/08/2010 05:59 PM, Julien Viarre wrote:

> Dear all,
>
> Sorry for the disturbance but I have a problem calculating anisotropic
> 3D-variograms with gstat package.
>
> I simulated a rainfall squall line propagating from West to East with a
> velocity of 1km/min (figure 1 attached).
> position.x represents longitude values
> position.y represents latitude values
> position.t represents rescaled timesteps (time is considered as the third
> dimension)
> For each time step, only one longitude band receives rainfall (other
> position have no rain) (positive rainfall values don't change)
>
> Then I calculated anisotropic variogram in 3 dimensions, specifically with
> alpha = 0° (meridian direction) and alpha = 90° (zonal direction) ; beta =
> 0° (no velocity) and beta = 45° (velocity of the squall line simulated)
>
> vario = variogram(estim.d[,4]~1, location= ~position.x+position.y+position.t
> , data = estim.d , alpha=angle.p[kp] , beta=angle.q[kq] , tol.hor=0 ,
> tol.ver=0, boundaries=classe_vario)
>
> For alpha = 90° and beta = 45°, estimations should be correlated so I should
> expect a variogram with gamma=0 for all the classes...
> Nevertheless, I obtained :
>
>     np       dist        gamma dir.hor dir.ver   id
>
> 1  30150  0.07071068  0.5773396664      90      45 var1
>
> 2  26550  0.14142136  0.5196323354      90      45 var1
>
> 3  47400  0.24584641  0.5358185764      90      45 var1
>
> 4   7550          0.35355339  0.5773182290      90      45 var1
>
> 5  32150  0.44746773  0.4290722330      90      45 var1
>
> 6  16300  0.56568542  0.4675393932      90      45 var1
>
> 7  16650  0.63639610  0.3920274991      90      45 var1
>
> 8  20050  0.75806784  0.3445372217      90      45 var1
>
> 9  12500  0.84852814  0.2323875687      90      45 var1
>
> 10 14050       0.93936997  0.0779285188      90      45 var1
>
> 11  8450  1.06066017  0.0003281226      90      45 var1
>
> 12  5000  1.13137085 0.0003727332      90      45 var1
>
> 13 10850       1.23336367  0.0004393653      90      45 var1
>
> 14  3150  1.34350288  0.0005170953      90      45 var1
>
> 15  3600  1.45644355  0.0005944342      90      45 var1
>
> 16  1800  1.55563492 0.0006728895      90      45 var1
>
> 17  1950  1.65172892  0.0007428695      90      45 var1
>
> 18   300          1.76776695  0.0007899247      90      45 var1
>
> 19   100          1.83847763 0.0007689443      90      45 var1
>
> 20    50          1.90918831 0.0005528209      90      45 var1
>
>
> That is very surprising because for the first classes we have positive
> semivariance values.
>
>
> And the most surprising is when I removed zero values to keep only the
> meridian positive band (figure 2 attached), I obtained consistent variogram
> estimation :
>
>
>
> Variogram associated :
>
>
>
>    np       dist gamma dir.hor dir.ver   id
>
> 1   600 0.07071068     0      90      45 var1
>
> 2   525 0.14142136     0      90      45 var1
>
> 3  1050 0.24917096     0      90      45 var1
>
> 4   150 0.35355339     0      90      45 var1
>
> 5   725 0.44620876     0      90      45 var1
>
> 6   350 0.56568542     0      90      45 var1
>
> 7   450 0.63639610     0      90      45 var1
>
> 8   575 0.75937120     0      90      45 var1
>
> 9   400 0.84852814     0      90      45 var1
>
> 10  475 0.93784689     0      90      45 var1
>
> 11  325 1.06066017     0      90      45 var1
>
> 12  200 1.13137085     0      90      45 var1
>
> 13  525 1.23575328     0      90      45 var1
>
> 14  175 1.34350288     0      90      45 var1
>
> 15  275 1.45921127     0      90      45 var1
>
> 16  150 1.55563492     0      90      45 var1
>
> 17  200 1.65286210     0      90      45 var1
>
> 18   50 1.76776695     0      90      45 var1
>
> 19   25 1.83847763     0      90      45 var1
>
> 20   25 1.90918831     0      90      45 var1
>
>
>
> Any idea?
>
> Thanks in advance,
>
> With best regards,
>
> Julien Viarre
> [hidden email]
>
>
>
>
> _______________________________________________
> Geostatistics mailing list
> [hidden email]
> http://list.52north.org/mailman/listinfo/geostatistics
> http://geostatistics.forum.52north.org

--
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      [hidden email]
_______________________________________________
Geostatistics mailing list
[hidden email]
http://list.52north.org/mailman/listinfo/geostatistics
http://geostatistics.forum.52north.org
Reply | Threaded
Open this post in threaded view
|

Re: [52N Geostatistics] Gstat package. Problem with 3d anisotropic variogram

Julien Viarre
Thank you very much for spending time to answer me!
I'm not aware of this problem.. Do you have some details? (I'm going  
to check on the net)
++

Le 22 oct. 10 à 12:23, Edzer Pebesma a écrit :

> Julien, thanks for you detailed email. I hope to get back to you  
> soon on
> this issue, but it will take some time. Putting my brains in the 3D
> anisotropy "mode" has always been quite an exercise, in the past.
>
> In the meantime: you are aware of the errors made in the original  
> GSLIB,
> which were copied into gstat?
>
> On 10/08/2010 05:59 PM, Julien Viarre wrote:
>> Dear all,
>>
>> Sorry for the disturbance but I have a problem calculating  
>> anisotropic
>> 3D-variograms with gstat package.
>>
>> I simulated a rainfall squall line propagating from West to East  
>> with a
>> velocity of 1km/min (figure 1 attached).
>> position.x represents longitude values
>> position.y represents latitude values
>> position.t represents rescaled timesteps (time is considered as the  
>> third
>> dimension)
>> For each time step, only one longitude band receives rainfall (other
>> position have no rain) (positive rainfall values don't change)
>>
>> Then I calculated anisotropic variogram in 3 dimensions,  
>> specifically with
>> alpha = 0° (meridian direction) and alpha = 90° (zonal direction) ;  
>> beta =
>> 0° (no velocity) and beta = 45° (velocity of the squall line  
>> simulated)
>>
>> vario = variogram(estim.d[,4]~1, location= ~position.x+position.y
>> +position.t
>> , data = estim.d , alpha=angle.p[kp] , beta=angle.q[kq] , tol.hor=0 ,
>> tol.ver=0, boundaries=classe_vario)
>>
>> For alpha = 90° and beta = 45°, estimations should be correlated so  
>> I should
>> expect a variogram with gamma=0 for all the classes...
>> Nevertheless, I obtained :
>>
>>    np       dist        gamma dir.hor dir.ver   id
>>
>> 1  30150  0.07071068  0.5773396664      90      45 var1
>>
>> 2  26550  0.14142136  0.5196323354      90      45 var1
>>
>> 3  47400  0.24584641  0.5358185764      90      45 var1
>>
>> 4   7550          0.35355339  0.5773182290      90      45 var1
>>
>> 5  32150  0.44746773  0.4290722330      90      45 var1
>>
>> 6  16300  0.56568542  0.4675393932      90      45 var1
>>
>> 7  16650  0.63639610  0.3920274991      90      45 var1
>>
>> 8  20050  0.75806784  0.3445372217      90      45 var1
>>
>> 9  12500  0.84852814  0.2323875687      90      45 var1
>>
>> 10 14050       0.93936997  0.0779285188      90      45 var1
>>
>> 11  8450  1.06066017  0.0003281226      90      45 var1
>>
>> 12  5000  1.13137085 0.0003727332      90      45 var1
>>
>> 13 10850       1.23336367  0.0004393653      90      45 var1
>>
>> 14  3150  1.34350288  0.0005170953      90      45 var1
>>
>> 15  3600  1.45644355  0.0005944342      90      45 var1
>>
>> 16  1800  1.55563492 0.0006728895      90      45 var1
>>
>> 17  1950  1.65172892  0.0007428695      90      45 var1
>>
>> 18   300          1.76776695  0.0007899247      90      45 var1
>>
>> 19   100          1.83847763 0.0007689443      90      45 var1
>>
>> 20    50          1.90918831 0.0005528209      90      45 var1
>>
>>
>> That is very surprising because for the first classes we have  
>> positive
>> semivariance values.
>>
>>
>> And the most surprising is when I removed zero values to keep only  
>> the
>> meridian positive band (figure 2 attached), I obtained consistent  
>> variogram
>> estimation :
>>
>>
>>
>> Variogram associated :
>>
>>
>>
>>   np       dist gamma dir.hor dir.ver   id
>>
>> 1   600 0.07071068     0      90      45 var1
>>
>> 2   525 0.14142136     0      90      45 var1
>>
>> 3  1050 0.24917096     0      90      45 var1
>>
>> 4   150 0.35355339     0      90      45 var1
>>
>> 5   725 0.44620876     0      90      45 var1
>>
>> 6   350 0.56568542     0      90      45 var1
>>
>> 7   450 0.63639610     0      90      45 var1
>>
>> 8   575 0.75937120     0      90      45 var1
>>
>> 9   400 0.84852814     0      90      45 var1
>>
>> 10  475 0.93784689     0      90      45 var1
>>
>> 11  325 1.06066017     0      90      45 var1
>>
>> 12  200 1.13137085     0      90      45 var1
>>
>> 13  525 1.23575328     0      90      45 var1
>>
>> 14  175 1.34350288     0      90      45 var1
>>
>> 15  275 1.45921127     0      90      45 var1
>>
>> 16  150 1.55563492     0      90      45 var1
>>
>> 17  200 1.65286210     0      90      45 var1
>>
>> 18   50 1.76776695     0      90      45 var1
>>
>> 19   25 1.83847763     0      90      45 var1
>>
>> 20   25 1.90918831     0      90      45 var1
>>
>>
>>
>> Any idea?
>>
>> Thanks in advance,
>>
>> With best regards,
>>
>> Julien Viarre
>> [hidden email]
>>
>>
>>
>>
>> _______________________________________________
>> Geostatistics mailing list
>> [hidden email]
>> http://list.52north.org/mailman/listinfo/geostatistics
>> http://geostatistics.forum.52north.org
>
> --
> 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      [hidden email]

_______________________________________________
Geostatistics mailing list
[hidden email]
http://list.52north.org/mailman/listinfo/geostatistics
http://geostatistics.forum.52north.org
Reply | Threaded
Open this post in threaded view
|

Re: [52N Geostatistics] Gstat package. Problem with 3d anisotropic variogram

Edzer Pebesma
Administrator
If you read the documentation of ?vgm, or look into its source at:

https://svn.52north.org/svn/geostatistics/main/gstat/S/man/vgm.Rd

you'll find the relevant information. The link it refers to is not
active anymore, regrettably. Googling "gslib errors" gave me:

http://www.gslib.com/sec_gb.html

which contains similar information. Let us know if this helped.

On 10/22/2010 06:42 PM, Julien Viarre wrote:

> Thank you very much for spending time to answer me!
> I'm not aware of this problem.. Do you have some details? (I'm going to
> check on the net)
> ++
>
> Le 22 oct. 10 à 12:23, Edzer Pebesma a écrit :
>
>> Julien, thanks for you detailed email. I hope to get back to you soon on
>> this issue, but it will take some time. Putting my brains in the 3D
>> anisotropy "mode" has always been quite an exercise, in the past.
>>
>> In the meantime: you are aware of the errors made in the original GSLIB,
>> which were copied into gstat?
>>
>> On 10/08/2010 05:59 PM, Julien Viarre wrote:
>>> Dear all,
>>>
>>> Sorry for the disturbance but I have a problem calculating anisotropic
>>> 3D-variograms with gstat package.
>>>
>>> I simulated a rainfall squall line propagating from West to East with a
>>> velocity of 1km/min (figure 1 attached).
>>> position.x represents longitude values
>>> position.y represents latitude values
>>> position.t represents rescaled timesteps (time is considered as the
>>> third
>>> dimension)
>>> For each time step, only one longitude band receives rainfall (other
>>> position have no rain) (positive rainfall values don't change)
>>>
>>> Then I calculated anisotropic variogram in 3 dimensions, specifically
>>> with
>>> alpha = 0° (meridian direction) and alpha = 90° (zonal direction) ;
>>> beta =
>>> 0° (no velocity) and beta = 45° (velocity of the squall line simulated)
>>>
>>> vario = variogram(estim.d[,4]~1, location=
>>> ~position.x+position.y+position.t
>>> , data = estim.d , alpha=angle.p[kp] , beta=angle.q[kq] , tol.hor=0 ,
>>> tol.ver=0, boundaries=classe_vario)
>>>
>>> For alpha = 90° and beta = 45°, estimations should be correlated so I
>>> should
>>> expect a variogram with gamma=0 for all the classes...
>>> Nevertheless, I obtained :
>>>
>>>    np       dist        gamma dir.hor dir.ver   id
>>>
>>> 1  30150  0.07071068  0.5773396664      90      45 var1
>>>
>>> 2  26550  0.14142136  0.5196323354      90      45 var1
>>>
>>> 3  47400  0.24584641  0.5358185764      90      45 var1
>>>
>>> 4   7550          0.35355339  0.5773182290      90      45 var1
>>>
>>> 5  32150  0.44746773  0.4290722330      90      45 var1
>>>
>>> 6  16300  0.56568542  0.4675393932      90      45 var1
>>>
>>> 7  16650  0.63639610  0.3920274991      90      45 var1
>>>
>>> 8  20050  0.75806784  0.3445372217      90      45 var1
>>>
>>> 9  12500  0.84852814  0.2323875687      90      45 var1
>>>
>>> 10 14050       0.93936997  0.0779285188      90      45 var1
>>>
>>> 11  8450  1.06066017  0.0003281226      90      45 var1
>>>
>>> 12  5000  1.13137085 0.0003727332      90      45 var1
>>>
>>> 13 10850       1.23336367  0.0004393653      90      45 var1
>>>
>>> 14  3150  1.34350288  0.0005170953      90      45 var1
>>>
>>> 15  3600  1.45644355  0.0005944342      90      45 var1
>>>
>>> 16  1800  1.55563492 0.0006728895      90      45 var1
>>>
>>> 17  1950  1.65172892  0.0007428695      90      45 var1
>>>
>>> 18   300          1.76776695  0.0007899247      90      45 var1
>>>
>>> 19   100          1.83847763 0.0007689443      90      45 var1
>>>
>>> 20    50          1.90918831 0.0005528209      90      45 var1
>>>
>>>
>>> That is very surprising because for the first classes we have positive
>>> semivariance values.
>>>
>>>
>>> And the most surprising is when I removed zero values to keep only the
>>> meridian positive band (figure 2 attached), I obtained consistent
>>> variogram
>>> estimation :
>>>
>>>
>>>
>>> Variogram associated :
>>>
>>>
>>>
>>>   np       dist gamma dir.hor dir.ver   id
>>>
>>> 1   600 0.07071068     0      90      45 var1
>>>
>>> 2   525 0.14142136     0      90      45 var1
>>>
>>> 3  1050 0.24917096     0      90      45 var1
>>>
>>> 4   150 0.35355339     0      90      45 var1
>>>
>>> 5   725 0.44620876     0      90      45 var1
>>>
>>> 6   350 0.56568542     0      90      45 var1
>>>
>>> 7   450 0.63639610     0      90      45 var1
>>>
>>> 8   575 0.75937120     0      90      45 var1
>>>
>>> 9   400 0.84852814     0      90      45 var1
>>>
>>> 10  475 0.93784689     0      90      45 var1
>>>
>>> 11  325 1.06066017     0      90      45 var1
>>>
>>> 12  200 1.13137085     0      90      45 var1
>>>
>>> 13  525 1.23575328     0      90      45 var1
>>>
>>> 14  175 1.34350288     0      90      45 var1
>>>
>>> 15  275 1.45921127     0      90      45 var1
>>>
>>> 16  150 1.55563492     0      90      45 var1
>>>
>>> 17  200 1.65286210     0      90      45 var1
>>>
>>> 18   50 1.76776695     0      90      45 var1
>>>
>>> 19   25 1.83847763     0      90      45 var1
>>>
>>> 20   25 1.90918831     0      90      45 var1
>>>
>>>
>>>
>>> Any idea?
>>>
>>> Thanks in advance,
>>>
>>> With best regards,
>>>
>>> Julien Viarre
>>> [hidden email]
>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> Geostatistics mailing list
>>> [hidden email]
>>> http://list.52north.org/mailman/listinfo/geostatistics
>>> http://geostatistics.forum.52north.org
>>
>> --
>> 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      [hidden email]
>

--
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      [hidden email]
_______________________________________________
Geostatistics mailing list
[hidden email]
http://list.52north.org/mailman/listinfo/geostatistics
http://geostatistics.forum.52north.org