gstat variogram & great circle distance
Hi Edzer, I am still having some trouble with the great distance calculations in 'variogram'. Your suggestion below works, but the distances are not correct (at least, not in kilometers, meters, or miles). I do not have proj.4 or gdal libraries installed, nor do I have the R packages proj4 or rgdal. It seems like I should not need them, as I am not doing a projection. Given that I am setting the 'proj4string' attribute in order to achieve a great cicle distance calculation, though, I wonder if not having those packages installed is a problem. An example (dump output for foo is below):
foo
z lon lat 1 -8.2920352 -68.74028 45.20407 2 0.3962574 -157.40894 70.46961 3 -5.3976371 -89.97919 46.08268
coordinates(foo) <- ~lon+lat
proj4string(foo)=CRS("+longlat")
variogram(z~1, foo, cloud=T, cutoff=10000)
dist gamma dir.hor dir.ver id left right
1 9880.046 37.74321 0 0 var1 2 1
2 2366.200 4.18877 0 0 var1 3 1
The distance from 1 to 2 is 5299 km, and from 3 to 1 is 1650 km. I'm not sure what the 9880.046 and 2366.200 represent. The ratio of the variogram$dist values to the correct distances in km are not the same, so those values cannot both be correct distances in any units.
Perhaps I am missing a package that gstat needs?
Many thanks for your help.
Cheers,
Tim
=======
foo <-
structure(list(z = c(-8.29203519866722, 0.396257381218808, -5.39763713302683), lon = c(-68.740278, -157.408944, -89.97919), lat = c(45.20407,70.469611, 46.08268)), .Names = c("z", "lon", "lat"), class = "data.frame", row.names = c(NA,3L))
=======
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 at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo