How do I estimate a variogram for data on a globe?

Question

I have data with associated longitudes and latitudes. How do I get a variogram for this data based on the great-circle distances between the points?

This simple example has all the data on the equator:

require(geoR)

long <- seq(-179, 180)
x <- sin(pi * long / 180) + rnorm(length(long))
V <- variog(data=x, coords=cbind(long, 0))
# variog: computing omnidirectional variogram
plot(V)

The first and last points are actually only 1 degree apart, but my naive attempt results in variog thinking they're separated by 359 degrees.

Answer 1

You should use the Semi-variogram from nmle. It allows you to specify a distance matrix, which you can trivially work out for yourself.

Answer 2

From a post on R-sig-geo (mailing list dedicated to spatial data in R) I seem to remember that there are no ready-to-go functions in R that support great circle distances:

http://r-sig-geo.2731867.n2.nabble.com/Great-Circle-distances-in-Automap-Gstat-td6863940.html

My suggestion would be to project your data and than perform interpolation on the projected data.