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([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 >> > -- 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: [EMAIL PROTECTED] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo