[Gstat-info] retrieving gammas from a variogram model

[Paul Hiemstra]

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
>
>
>   


-- 
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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: s_custom.c
Type: text/x-csrc
Size: 31149 bytes
Desc: not available
Url : http://mailman.geo.uu.nl/pipermail/gstat-info/attachments/20080715/7c7774ac/s_custom-0001.bin