[R-sig-Geo] aggregate outcome kernel density estimation using saga gis and overlay

[Marco Helbich]

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=win at polygons[[1]]@Polygons[[1]]@coords[,1],
  y=win at polygons[[1]]@Polygons[[1]]@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