Hello Edzer,
I was curious to run the example you provided (using the latest CRAN
versions of all packages), but as written the line
> idw.spdf = as(idw.out, "SpatialPolygonsDataFrame")
produces
Error in as(idw.out, "SpatialPolygonsDataFrame") : no method or default for
coercing "SpatialPixelsDataFrame" to "SpatialPolygonsDataFrame"
Regards,
Greg.
2008/10/30 Edzer Pebesma <[EMAIL PROTECTED]>
> 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 [EMAIL PROTECTED]
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
