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@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/
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
