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

[Edzer Pebesma]

Timothy, for some reason the projected argument was not meant to be set 
by users at this level of abstraction; I'll look into it. The following 
seems to work:

 > proj4string(foo)=CRS("+longlat")
 > proj4string(foo)
[1] "+longlat"
 > variogram(z~1,foo)
  np      dist        gamma dir.hor dir.ver   id
1  1  177.7815   0.48861387       0       0 var1
2  2  614.0040   1.11639574       0       0 var1
3  3  715.8402 115.50300578       0       0 var1
4  1  829.6047   0.01873678       0       0 var1
5  1  893.0483   3.66061567       0       0 var1
6  1  992.9436 268.46555052       0       0 var1
7  4 1095.7359  83.25299050       0       0 var1
8  1 1274.4497  12.33954872       0       0 var1
9  1 1357.4796   3.01500925       0       0 var1

--
Edzer


Timothy W. Hilton wrote:
> 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")
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>