Dear List, my problem is to aggregate the outcome frome kernel density estimations to some polygons using saga gis and/or the overlay command. In other words, applying a point-in-polygon analysis and calculating the mean of the points inside every polygon.
I hope the following (cumbersome) example-code helps to reproduce my problem. Thank you very much for your help! Regards Marco ---- Marco Helbich Institute for Urban and Regional Research Austrian Academy of Sciences Postgasse 7/4/2, A-1010 Vienna, Austria e-mail: marco.helbich(at)oeaw.ac.at e-mail priv.: marco.helbich(at)gmx.at *** Keep on rockin' in the free world *** ####### CODE ############# library(rgdal); library(splancs); library(RSAGA); library(sp) # test data sids <- readOGR(system.file("shapes/sids.shp", package = "maptools"), "sids") crs <- CRS("+proj=longlat +datum=NAD27") proj4string(sids) <- crs plot(sids) ncrs <- CRS("+proj=utm +zone=18 +datum=NAD27") tsids <- spTransform(sids, ncrs) win <- unionSpatialPolygons(tsids, rep("x", length(slot(tsids, "polygons")))) proj4string(win) <- ncrs tcoord <- cbind(coordinates(tsids)[,1], coordinates(tsids)[,2]) tpoi <- SpatialPoints(tcoord) proj4string(tpoi) <- ncrs plot(win) plot(tpoi, add=T) tpoi.pts <- as.points(coordinates(tpoi)) win.pol = as.points(list(x=...@polygons[[1]]@polygons[[...@coords[,1], y=...@polygons[[1]]@polygons[[...@coords[,2])) # define grid, bandwith for kernel density estimation pixs <- 15000 dimx <- ceiling(abs((bbox(win)[1,1] - bbox(win)[1,2])/pixs)) dimy <- ceiling(abs((bbox(win)[2,1] - bbox(win)[2,2])/pixs)) grd <- GridTopology(cellcentre.offset=c(bbox(win)[1,1], bbox(win)[2,1]), cellsize=c(pixs, pixs), cells.dim=c(dimx,dimy)) summary(grd) # kernel density estimation kde <- spkernel2d(tpoi.pts, win.pol, 2000, grd) kdf <- data.frame(k0=kde) k <- SpatialGridDataFrame(grd, data=kdf) spplot(k, sp.layout = list(list("sp.polygons", tsids))) k1 <- as(k, "SpatialPointsDataFrame") k1[["x"]] <- coordinates(k1)[,1] k1[["y"]] <- coordinates(k1)[,2] # the following part shows my problem: how can i aggregate the estimations # (using the mean) for every polygon in tsids? # version 1: using saga... # i can not find any command to import the shape (have i missed somthing ?), # i would like to calculated the point in polygon analysis there and export it # back to R... writeOGR(k1, dsn=".", layer="test", driver="ESRI Shapefile") #write.dbf(slot(k1, "data"), file = "testdbf") rsaga.env(path="C:/Progra~1/saga_vc") rsaga.get.modules("shapes_polygons") rsaga.get.usage("shapes_polygons", 4) # version 2: using overlay xx <- overlay(k1, tsids, fn = mean) # how can i combine the outcome with my transformed sids to plot the result? -- für nur 19,99 Euro/mtl.!* http://portal.gmx.net/de/go/dsl02 _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo