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

[Tomislav Hengl]

# 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: file7b0208e.kml
Type: application/vnd.google-earth.kml+xml
Size: 423 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20080127/852e8171/attachment.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: file7b0208e.png
Type: image/png
Size: 2964 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20080127/852e8171/attachment.png>