[R-sig-Geo] [Gstat-info] retrieving gammas from a variogram model

[Paul Hiemstra]

Thanks,

A good way of incorporating it.

Paul

Edzer Pebesma schreef:
> (I'm adding r-sig-geo to this thread)
>
> Paul/Ashton,
>
> I included your code, Paul; I also added the option (with the same 
> argument name) to variogramLine so you can now also do e.g. a
>
> variogramLine(vgm(1, "Sph", 1, anis=c(0, 0.5)), dir = 
> c(1/sqrt(c(2,2)),0), dist_vector = c(0, 0.1, 0.5, 1, 2))
>
> (that is, irregular distances at a given direction and anisotropic model)
> --
> Edzer
>
>
> Paul Hiemstra wrote:
>> Hi Ashton,
>>
>> I had the same problem and also tried your workaround. But the 
>> problem that that is really slow. To get past this I delved into te 
>> C-source code of gstat and reworked one of the C-functions to get 
>> what I wanted. This solution works much faster (few orders of 
>> magnitude).
>>
>> In order to get this working you need to reinstall gstat with one 
>> modified .c file (s.c, located in the src directory of the source 
>> version of the gstat package). I added the modified .c file as an 
>> attachment (hope that is possible on the mailing list). Once you 
>> replaced s.c you install gstat to get the new C-function available. 
>> The new C-function accepts a vector of distances as input and returns 
>> a vector of gammas. An example providing an R-wrapper function:
>>
>> getGammas = function(dist_vector, varModel, covariance = FALSE) {
>> # dist_vector is the input vector with distances
>> # varModel is the variogram model
>> # covariance is a logical to choose for semivariances (FALSE)
>> # or for covariances (TRUE)
>>   .Call("gstat_init", as.integer(0))
>>    gstat:::load.variogram.model(varModel)
>>    ret = .Call("get_covariance_list", as.integer(c(0,0)), 
>> as.integer(covariance), as.numeric(dist_vector))[[2]]
>>    .Call("gstat_exit", 0)
>>   ret   }
>> library(gstat)
>> data(meuse)
>> vgm1 <- variogram(log(zinc)~1, ~x+y, meuse)
>> m = fit.variogram(vgm1, vgm(1,"Sph",300,1))
>> getGammas(seq(1,1000, 100), varModel = m)
>>
>> It is experimental code but it works very well for me, much faster 
>> than the workaround with variogramLine. Please let me know what you 
>> think of this solution.
>>
>> cheers and hth,
>> Paul
>>
>> Ashton Shortridge wrote:
>>> I have a vector of lags for which I'd like to calculate covariances, 
>>> given a variogram (covariance) model. variogramLines() almost seems 
>>> to do this; however, it generates a table of distances and 
>>> covariances for a sequence of values between user-specified min and 
>>> max distances - the user can't pass a vector of h-values.
>>>
>>> I was hoping I could simply modify the R wrapper code for 
>>> variogramLines to do what I wished, but it seems that the relevant 
>>> code is called from elsewhere.
>>>
>>> My current workaround is to call variogramLines for each lag value, 
>>> specifying an npoints of two. For example, using a vgm object for an 
>>> omnidirectional variogram called testmod, and interested in the 
>>> covariance for a distance of 15, I run:
>>> gstat.covar <- variogramLine(testmod, 15, 2, covariance=TRUE)
>>>
>>> and then grab the second row in gstat.covar. This is tedious, and 
>>> I'm guessing a better way is out there!
>>>
>>> Thanks for any help,
>>>
>>> Ashton
>>>
>>>
>>>   
>>
>>
>> ------------------------------------------------------------------------
>>
>> _______________________________________________
>> Gstat-info mailing list
>> Gstat-info at geo.uu.nl
>> http://mailman.geo.uu.nl/mailman/listinfo/gstat-info
>
> -- 
> 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/
>   


-- 
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone:     +31302535773
Fax:    +31302531145
http://intamap.geo.uu.nl/~paul