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

[Lyndon Estes]

Hi Edzer,

Many thanks for the pointer.

My rough function is an order of magnitude slower than the 'over' one
you showed.

Sys.time() difference for my approach = 3.65 seconds, for yours =
0.323 seconds.

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

Thanks again, Lyndon


On Thu, Feb 17, 2011 at 2:07 AM, Edzer Pebesma
<edzer.pebesma at uni-muenster.de> wrote:
> 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>



-- 
Lyndon Estes
Research Associate
Woodrow Wilson School
Princeton University
+1-609-258-2392 (o)
+1-609-258-6082 (f)
+1-202-431-0496 (m)
lestes at princeton.edu