[R-sig-Geo] gstat::variogram - distance calculation is incorrect; kriging on sphere

[Edzer Pebesma]

Dear all,

In continuation of this thread, I've spent some time looking at kriging 
on the sphere, corrected some bugs, and need further help.

First of all, distances for covariances on the sphere were computed 
incorrectly as well in gstat, so I hope not too many people have been 
relying on this--kriging seemed to happen still in some Euclidian mode. 
The good news is that it seems to work now (gstat 0.9-53, accepted on 
CRAN). Inverse distance interpolation for spherical data seemed to work 
already.

Covariances on the sphere now work, but the models present do not 
include those specially deviced for spherical data. Can anyone provide 
me with or point me to useful, preferably simple covariance functions 
that are positive definite on the sphere? The example below seems to 
work but in certain cases without nugget the interpolation may go crazy. 
You'll need to update your sp to 0.9-27 to run it (also accepted on CRAN).

Below is an example script. It also needs the new sp and gstat versions.

library(gstat)
library(rgdal)
world = expand.grid(long=seq(-177.5,177.5,5),lat=seq(-87.5,87.5,5))
world.sp = SpatialPixels(SpatialPoints(world,CRS("+proj=longlat")))
plot(world.sp,axes=T)
pts=data.frame(long=runif(100,-180,180),lat=runif(100,-90,90),val=rnorm(100))
coordinates(pts)=~long+lat
proj4string(pts)=CRS("+proj=longlat")
points(pts,col='red')

# inverse distance interpolation on the sphere:
idw.out = idw(val~1,pts,world.sp)
image(idw.out, axes = TRUE, ylim = c(-90,90))
points(pts, pch=3)
idw.spdf = as(idw.out, "SpatialPolygonsDataFrame")
newproj = CRS("+proj=moll")
idw.spdf.moll = spTransform(idw.spdf, newproj)
spplot(idw.spdf.moll, "var1.pred",col.regions=bpy.colors(),col=0,
    sp.layout = list(sp.points, spTransform(pts, newproj), col = 'black'))

# kriging on the sphere
kr.out = krige(val~1,pts,world.sp,vgm(1, "Exp", 3000))
idw.spdf.moll$kr = kr.out[[1]]
spplot(idw.spdf.moll, "kr", col.regions=bpy.colors(), col=0,
    sp.layout = list(sp.points, spTransform(pts, newproj), col = 'black'))

-- 
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/
http://www.springer.com/978-0-387-78170-9 e.pebesma at wwu.de

-------------- next part --------------
A non-text attachment was scrubbed...
Name: world.png
Type: image/png
Size: 22219 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20081030/fbb4efe5/attachment.png>