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

[Roger Bivand]

On Sun, 27 Jan 2008, Tomislav Hengl wrote:

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

Thanks, very clear - the kernel values at the grid of projected points are 
not an a grid in geographical coordinates, and need to be interpolated 
(warped) to a grid in geographical coordinates. The alternative is to make 
the kernel function use Great Circle distances.

Roger

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

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