Count overlapping polygons in single shapefile with R

I have a shapefile of species ranges for ~ 346 species. Some species have multiple polygon ranges so there are 1,078 polygons and many polygons are multi-part. It looks something like this:

I would like to get the max count of overlapping polygons anywhere within the layer bounds. There are two ways I have considered accomplishing this, but I wanted to see if there is a simple way before I undertake the effort to try either. I would prefer to work in R, but could potentially use ArcGIS 10.2.

I tried using the 'Count Overlapping Polygons' toolbox from http://www.arcgis.com/home/item.html?id=1dd4a6832b3d40b494dbf8521cc5134c which seems like it would allow me to extract the information I need. However, it failed with an 'Out of memory' error. I have a Winx64 computer with 8GB of RAM and only had a few other processes running.

I came across this post by @Roger Bivand but it appears to only function between two shapefiles, which will not work for me.

The other approach that I thought might work is to rasterize each polygon and then sum all the rasters. I found a function here that allows the conversion to raster. I thought this might be a very time-consuming and computationally expensive, so I have not yet tried that approach.

Answer

In R, you can used the sp package and over function to do this. I adapted an example data set and the solution from this post by Roger Bivand.

Set up the example data:

library(sp)
library(rgeos)
library(raster)
library(rworldmap)


box <- readWKT("POLYGON((-180 90, 180 90, 180 -90, -180 -90, -180 90))")
proj4string(box) <- CRS("+proj=cea +datum=WGS84")
set.seed(1)
pts <- spsample(box, n=2000, type="random")
pols <- gBuffer(pts, byid=TRUE, width=50) # create circle polys around each point
merge = sample(1:40, 100, replace = T) # create vector of 100 rand #s between 0-40 to merge pols on

Sp.df1 <- gUnionCascaded(pols, id = merge) # combine polygons with the same 'merge' value
# create SPDF using polygons and randomly assigning 1 or 2 to each in the @data df

Sp.df <- SpatialPolygonsDataFrame(Sp.df1,
                    data.frame(z = factor(sample(1:2, length(Sp.df1), replace = TRUE)),
                                                    row.names= unique(merge)))
Sp.df <- crop(Sp.df, box)
colors <- c(rgb(r=0, g=0, blue=220, alpha=50, max=255), rgb(r=220, g=0, b=0, alpha=50, max=255))

land <- getMap()

overlay.map <- spplot(Sp.df, zcol = "z", col.regions = colors,
                   col = NA, alpha = 0.5, breaks=c(0,1),

                   sp.layout = list("sp.polygons", land, fill = "transparent",
                                    col = "grey50"))
overlay.map

Then to get the count (or other statistics)...

# find the count of polygons below each grid cell
GT <- GridTopology(c(-179.5, -89.5), c(1, 1), c(360, 180))
SG <- SpatialGrid(GT)
proj4string(SG) <- CRS("+proj=cea +datum=WGS84")



o <- over(SG, Sp.df1, returnList=TRUE)
ct <- sapply(o, length)
summary(ct)
SGDF <- SpatialGridDataFrame(SG, data=data.frame(ct=ct))
spplot(SGDF, "ct", col.regions=bpy.colors(20))