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

[Lyndon Estes]

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