Thanks, a really interesting post. May I ask a stupid question? Where are the following function documented?
getPolygonCoordsSlot(getPolygonsPolygonsSlot(getSpPpolygonsSlot(poly.shp )[[1]])[[1]]) They are in the sp package but I can't find references but in some posts on this ML. Thanks, Giovanni 2008/1/27, Roger Bivand <[EMAIL PROTECTED]>: > > 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 > [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
