This works great, Edzer. I expected there was a simple solution. Many thanks!
-Tim On Wed, Apr 2008, 30 at 07:42:36AM +0200, Edzer Pebesma wrote: > 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@stat.math.ethz.ch > >https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > > _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo