Many thanks for the pointer.

you showed.

0.323 seconds.

That will make a big difference when I run it for several thousand species.

> 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@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@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@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

>>

[hidden email]
>>

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 [hidden email]
>

> _______________________________________________

> R-sig-Geo mailing list

>

[hidden email]
>

https://stat.ethz.ch/mailman/listinfo/r-sig-geo>