[R-sig-Geo] gstat::variogram - distance calculation is incorrect

[Edzer Pebesma]

Thanks for the bug report; indeed, distance computation for longlat data 
variograms went wrong so far. I just uploaded a new gstat version to 
/incoming on CRAN that seems to repair this bug.

Next thing I'm not 100% sure of is the direction selections for 
directional variograms of longlat data.

I hope you don't mind I added your code below to one of the regression 
tests in the package (unproj.R).
--
Edzer

Timothy W. Hilton wrote:
> Hello,
>
> I am trying to use 'variogram' from the gstat package to compute an empirical semivariogram for a dataset that spans the North American continent.  I have been unable to get gstat to calculate the great circle distances correctly.  For example, in the example below, the first two lines of vg.foo suggest to me that the calculated distance between foo[1,] and foo[2,] is 4804, and between foo[1,] and foo[3,] is 3047.  I assume those units are km.  
>
> spDistsN1, however tells me that those distances should be 3115 km and 1907 km, respectively.
>
> Could someone suggest to me my error, or how I might use calculate an emprical variogram with correct distances?  I do not think UTM coordinates are an option for me, as my data span approximately 5000 km.
>
> Many thanks,
> Tim
>
> =====
>
> code to illustrate the problem:
>
> library(sp)
> library(rgdal)
> library(gstat)
>
> foo <-
> structure(list(z = c(-1.95824831109744, -1.9158901643563, 4.22211761150161, 
> 3.23356929459598, 1.12038389231868, 0.34613850821113, 1.12589932643631, 
> 23.517912251617, 3.0519158690268, 3.20261431141517, -2.10947106854739
> ), lon = c(-125.29228, -82.1556, -98.524722, -99.948333, -104.691741, 
> -79.420833, -105.100533, -88.291867, -72.171478, -121.556944, 
> -89.34765), lat = c(49.87217, 48.2167, 55.905833, 56.635833, 
> 53.916264, 39.063333, 48.307883, 40.0061, 42.537756, 44.448889, 
> 46.242017)), .Names = c("z", "lon", "lat"), row.names = c(NA, 
> -11L), class = "data.frame")
>
> coordinates(foo) <- ~lon+lat
>
> proj4string(foo) <- CRS('+proj=longlat')
>
> vg.foo <- variogram(z~1, foo, cloud=TRUE, cutoff=1e10)
>
> cat('==========\nvariogram:\n')
> print(head(vg.foo))
>
> cat('==========\nspDistsN1 Distances:\n')
> print(spDistsN1(coordinates(foo), coordinates(foo)[1,], longlat=T))
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> 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/