[R-sig-Geo] kernel density estimation to google earth

[Roger Bivand]

On Sat, 26 Jan 2008, Marco Helbich wrote:

> Dear List,
>
> I've calculated a kernel density estimation (splancs function) and now I
> want to export the result as a kml-file to open it in the google earth
> viewer, but I get stuck at converting the SpatialGridDataFrame to a
> SpatialPixelDataFrame...
>
> here (http://www.zshare.net/download/689742375ef5fb/) you can download some
> testdata and my R-code so far:
>
> ###########
> data.shp <- readOGR("C:/", layer="events")
> prj <- data.shp@ proj4string@ projargs
> dat <- data.shp
> str(dat)
> poly.shp <- readOGR("C:/", layer="hull")
> str(poly.shp)
>
> dat.SP <- as(dat, "SpatialPoints")
> str(dat.SP)
> pp_poi <- as.points(dat.SP at coords[,1], dat.SP at coords[,2])
> poly <-
> getPolygonCoordsSlot(getPolygonsPolygonsSlot(getSpPpolygonsSlot(poly.shp)[[1]])[[1]])
> polymap(poly)
> points(pp_poi)
>
> grd <- GridTopology(cellcentre.offset=c(590511, 396191), cellsize=c(2000,
> 2000),
>  cells.dim=c(30,25))

Here cellsize= is in metres, cells.dim= the numbers of columns and rows.

> kbw2000 <- spkernel2d(pp_poi, poly, h0=2000, grd)
> spplot(SpatialGridDataFrame(grd, data=data.frame(kbw2000)),
> col.regions=terrain.colors(16))
>
> test <- SpatialGridDataFrame(grd, data=data.frame(kbw2000))
> proj4string(test) <- CRS(prj)
> str(test)
> test1 <- spTransform(test, CRS("+proj=longlat +datum=WGS84"))
>
> ## here I got following error message:
> # validityMethod(as(object, superClass)): Geographical CRS given to
> non-conformant data
> test2 <- spsample(test1, type="regular", cellsize=c(2000,2000))

But test1 is now in geographical coordinates, so 2K degrees is a bit 
extreme, and the validation mechanism traps you. I'm not sure what your 
cellsize and cell.dims should be here in longlat.

>
> # export the SpatialPixelDataFrame (Code (not testet) from the help-file)
> tf <- tempfile()
> SGxx <- GE_SpatialGrid(test2)
> png(file=paste(tf, ".png", sep=""), width=SGxx$width, height=SGxx$height,
>  bg="transparent")
> par(mar=c(0,0,0,0), xaxs="i", yaxs="i")
> plot(x, xlim=SGxx$xlim, ylim=SGxx$ylim, setParUsrBB=TRUE)

There is no x defined here.

Roger

> dev.off()
> kmlOverlay(SGxx, paste(tf, ".kml", sep=""), paste(tf, ".png", sep=""))
> ###########
>
> I appreciate every hint! Thanks.
> Marco
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no