kernel density estimation to google earth

# Modified by T. Hengl (http://spatial-analyst.net)
# 27.01.2008

library(rgdal)
library(maptools)
library(splancs)

# Import the points and study area:

data.shp <- readOGR("C:/", layer="events")
str(data.shp)

poly.shp <- readOGR("C:/", layer="hull")
str(poly.shp)

poly <-
getPolygonCoordsSlot(getPolygonsPolygonsSlot(getSpPpolygonsSlot(poly.shp)[[1]])[[1]])
grd <-
GridTopology(cellcentre.offset=c(round(poly.shp at bbox["r1","min"],0),
round(poly.shp at bbox["r2","min"],0)), cellsize=c(2000,2000),
cells.dim=c(30,25))

# Run the 2D kernel smoother:

kbw2000 <- spkernel2d(data.shp, poly, h0=2000, grd)
hist(kbw2000)

# Pack and plot a SpatialGridDataFrame:

kbw2000.grd <- SpatialGridDataFrame(grd, data=data.frame(kbw2000))
proj4string(kbw2000.grd) <- data.shp at proj4string
spplot(kbw2000.grd, col.regions=terrain.colors(16),
sp.layout=list("sp.points",pch="+",cex=1.2,col="black",data.shp))

# Export to KML
# First, reproject the grid to longlat:

kbw2000.ll <- spTransform(kbw2000.grd, CRS("+proj=longlat +datum=WGS84"))
str(kbw2000.ll)

# The cell size you need to determine yourself!!

width =
(kbw2000.grd at bbox["coords.x1","max"]-kbw2000.grd at bbox["coords.x1","min"])/2000
height =
(kbw2000.grd at bbox["coords.x2","max"]-kbw2000.grd at bbox["coords.x2","min"])/2000
geogrd.cell = (kbw2000.ll at bbox["s1","max"]-kbw2000.ll at bbox["s1","min"])/width

# Define a new grid:

geogrd = spsample(kbw2000.ll, type="regular",
cellsize=c(geogrd.cell,geogrd.cell))
gridded(geogrd) = TRUE

gridparameters(geogrd)
# cellcentre.offset   cellsize cells.dim
# x1          15.90165 0.02636685        30
# x2          47.95541 0.02636685        16

# This is an empty grid without any topology (only grid nodes are defined)
and coordinate
# system definition. To create topology, we coerce a dummy variable (1s),
then
# specify that the layer has a full topology:

nogrids = geogrd at grid@cells.dim["x1"]*geogrd at grid@cells.dim["x2"]
geogrd = SpatialGridDataFrame(geogrd at grid, data=data.frame(rep(1,
nogrids)), proj4string=kbw2000.ll at proj4string)

# and estimate the values of the reprojected map at new grid locations
using the bilinear resampling:
# this can be time-consuming for large grids!!!

library(gstat)
kbw2000.llgrd = krige(kbw2000~1, kbw2000.ll, geogrd, nmax=4)

# Optional, convert the original shape to latlong coordinates:

data.ll <- spTransform(data.shp, CRS("+proj=longlat +datum=WGS84"))
spplot(kbw2000.llgrd["var1.pred"], col.regions=terrain.colors(16),
scales=list(draw=TRUE),
sp.layout=list("sp.points",pch="+",cex=1.2,col="black",data.ll))

# The final grid map can be exported to KML format using the maptools
package and kmlOverlay method:

kbw2000.kml = GE_SpatialGrid(kbw2000.llgrd)

tf <- tempfile()
png(file=paste(tf, ".png", sep=""), width=kbw2000.kml$width,
height=kbw2000.kml$height, bg="transparent")

par(mar=c(0,0,0,0), xaxs="i", yaxs="i")
image(as.image.SpatialGridDataFrame(kbw2000.llgrd[1]), col=bpy.colors(),
xlim=kbw2000.kml$xlim, ylim=kbw2000.kml$ylim)
plot(data.ll, pch="+", cex=1.2, add=TRUE, bg="transparent")

kmlOverlay(kbw2000.kml, paste(tf, ".kml", sep=""), paste(tf, ".png", sep=""))
dev.off()

# see also:
# Hengl, T., 2007. A Practical Guide to Geostatistical Mapping of
Environmental Variables. EUR 22904 EN Scientific and Technical Research
series, Office for Official Publications of the European Communities,
Luxemburg, 143 pp.

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))
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))

# 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)
dev.off()
kmlOverlay(SGxx, paste(tf, ".kml", sep=""), paste(tf, ".png", sep=""))
###########

I appreciate every hint! Thanks.
Marco

--
Marco Helbich
Institute for Urban and Regional Research
Austrian Academy of Sciences
Postgasse 7/4/2, A-1010 Vienna, Austria (EU)
e-mail: marco.helbich(at)oeaw.ac.at