[R-sig-Geo] Efficient way to obtain gridded count of overlapping polygons?

[Edzer Pebesma]

Yes, the (relatively) new method "over" does this for you, as in

pts$n = sapply(over(pts, geometry(sppX3), returnList = TRUE), length)

In the context of a running script:

require(sp)
con <- url(
"http://sites.google.com/site/ldemisc/polygon-intersect/spp.shp.Rdata")
load(con)
close(con)
proj4string(pts) = proj4string(sppX3) # otherwise over will fail here;
pts$n = sapply(over(pts, geometry(sppX3), returnList = TRUE), length)
gridded(pts)=TRUE
spplot(pts["n"], sp.layout=list("sp.polygons", sppX3, first=F))


On 02/17/2011 04:58 AM, Lyndon Estes wrote:
> Hello,
> 
> I am interested in creating a raster that summarizes the number of
> overlapping polygons that underlie each cell (in order to display how
> many species occur at that point).
> 
> I have come up with a function that seems to work, but I was wondering
> if I had missed one that already exists, or if someone can recommend a
> better way?
> 
> Herewith the code:
> 
> 
> CountIntersections <- function(x, y, CID) {
> # Counts number of polygons intersecting spatial points, specifically
> designed for situations in which
> # SpatialPolygons object contains overlapping polygons
> # Args:
> #   x: Input SpatialPointsDataFrame
> #   y: SpatialPolygons* on which to count intersects
> #   CID: Name of column containing unique polygon IDs
> # Returns:
> #   SpatialPointsDataFrame with column providing summary of number of
> polygons with unique IDs intersected
> 
>   # Vector of unique IDs
>   ids <- as.character(unique(y at data[, CID]))
> 
>   # lapply function in which each unique polygons is extracted and
> overlaid by points, with 1 assigned for sp
>   # presence, 0 for absence
>   p.int <- lapply(ids, function(z) {
>     p.ex <- y[y at data[, CID] == z, ]  # Extract polygon(s)
> corresponding to unique ID extracted by lapply
>     ov <- overlay(x, p.ex)  # Intersect points with sp polygon
>     ov2 <- ov / ov  # Reduce intersected polygon IDs to 1
>     ov2[is.na(ov2)] <- 0  # Convert NAs to 0
>     return(ov2)
>   })
>   x$P.cnt <- Reduce("+", p.int)  # Sum list of intersect for each
> species, and add value to dataframe of x
>   return(x)
> }
> 
> # Load spatial points (pts) and polygons of 3 species' ranges (sppX3)
> extracted from IUCN mammals database
> # downloaded here: http://www.iucnredlist.org/technical-documents/spatial-data
> con <- url("http://sites.google.com/site/ldemisc/polygon-intersect/spp.shp.Rdata")
> load(con)
> close(con)
> 
> # Run function and convert to raster
> pts.ct <- CountIntersections(pts, sppX3, "BINOMIAL")
> pts.grd <- SpatialPixelsDataFrame(pts.ct, pts.ct at data)
> pts.rst <- raster(pts.grd, layer = "P.cnt")
> 
> # Seems to work based on plots
> plot(pts.rst)
> plot(sppX3, add = T)
> 
> Thanks in advance for any advice.
> 
> Cheers, Lyndon
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      e.pebesma at wwu.de