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
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
