# 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([EMAIL PROTECTED]"r1","min"],0), round([EMAIL PROTECTED]"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) <- [EMAIL PROTECTED] 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 = ([EMAIL PROTECTED]"coords.x1","max"[EMAIL PROTECTED]"coords.x1","min"])/2000 height = ([EMAIL PROTECTED]"coords.x2","max"[EMAIL PROTECTED]"coords.x2","min"])/2000 geogrd.cell = ([EMAIL PROTECTED]"s1","max"[EMAIL PROTECTED]"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 = [EMAIL PROTECTED]@cells.dim["x1"[EMAIL PROTECTED]@cells.dim["x2"] geogrd = SpatialGridDataFrame([EMAIL PROTECTED], data=data.frame(rep(1, nogrids)), [EMAIL PROTECTED]) # 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([EMAIL PROTECTED],1], [EMAIL PROTECTED],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@stat.math.ethz.ch > https://stat.ethz.ch/mailman/listinfo/r-sig-geo >
file7b0208e.kml
Description: application/vnd.google-earth.kml
<<attachment: file7b0208e.png>>
_______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo