[R-sig-Geo] gstat variogram & great circle distance

[Timothy W. Hilton]

Hello,

I am trying to use gstat to compute a semivariogram for data whose coordinates are latitude/longitude pairs.  I would like to use the great circle distance between pairs.  The documentation implies that gstat can do this, but I am not having any success.  If anyone could suggest the correct syntax, I would greatly appreciate it.

Here is a sample of my data (see output from dump below):


> foo
            z        lon      lat
1  -1.9582483 -125.29228 49.87217
2  -1.9158902  -82.15560 48.21670
3   4.2221176  -98.52472 55.90583
4   3.2335693  -99.94833 56.63583
5   1.1203839 -104.69174 53.91626
6   0.3461385  -79.42083 39.06333
7   1.1258993 -105.10053 48.30788
8  23.5179123  -88.29187 40.00610
9   3.0519159  -72.17148 42.53776
10  3.2026143 -121.55694 44.44889
11 -2.1094711  -89.34765 46.24202

I can calculate a variogram:

>coordinates(foo) <- ~lon+lat
>variogram(z~1, foo)
   np      dist       gamma dir.hor dir.ver   id
1   1  1.599865   0.4886139       0       0 var1
2   2  5.545490   1.1163957       0       0 var1
3   4  6.712018  86.6319381       0       0 var1
4   1  8.038953   3.6606156       0       0 var1
5   3  9.422337  91.0816908       0       0 var1
6   2 10.149322 164.1162183       0       0 var1
7   2 11.868366   7.6772788       0       0 var1
8   1 13.326965  20.0445076       0       0 var1
9   1 14.846073  14.2740402       0       0 var1
10  1 15.887767   5.2338108       0       0 var1
11  3 16.792331  72.2669527       0       0 var1
12  2 17.828085  16.0787636       0       0 var1

The distances are clearly not great circle distances, though.  Setting the "projected" flag to "false" gives me this error:

> variogram(z~1, foo, projected=FALSE)
Error in variogram.default(y, locations, X, trend.beta = beta, grid = grid,  :
  formal argument "projected" matched by multiple actual arguments

Thanks in advance for any help,
Tim

==================
`foo` <-
structure(list(z = c(-1.95824831109744, -1.91589016435630, 4.22211761150161,
3.23356929459598, 1.12038389231868, 0.34613850821113, 1.12589932643631,
23.5179122516170, 3.05191586902680, 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")