5

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:

enter image description here

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.

  • Have you checked how much RAM you have allocated to ArcGIS? – I Heart Beats Oct 7 '14 at 13:44
  • With regards to a possible solution, I am confused as to how you have 2 dimensional overlapping polygons in a single shape file. If you had multiple shapefiles where polygons overlapped, you can use the Union Tool to obtain the overlaps and count those using something like excel, once you pulled the attribute file data. Can you provide an image of your shape file? Union tool: resources.arcgis.com/en/help/main/10.2/index.html#//… – I Heart Beats Oct 7 '14 at 13:52
  • 1
    There is also this solution, which perhaps can be adapted to your problem: "Shapefile with overlapping polygons: calculate average values" stackoverflow.com/questions/20949926/… – I Heart Beats Oct 7 '14 at 13:53
  • Thanks for the suggestions @IHeartBeats. I haven't checked the RAM allocation. Can you give instructions how? I'm on Windows 8.1 with arcGIS 10.2.0. Neither have I installed the 64-bit background processing, so that might be a way to solve the issue for using ArcGIS. I'm not sure what you mean regarding '2 dimensional overlapping polygons in a single shape file.` Indeed, I could split the polygons into 1078 individual shapefiles, but that seems like it would be better to avoid. I uploaded an image similar to my shape file above. – Cotton.Rockwood Oct 9 '14 at 6:38
  • Since ArcGIS was failing and I prefer to keep everything in R anyway, if possible, I looked for a solution in R. I was able to find what seems like a reasonable and flexible solution. See the answer below. – Cotton.Rockwood Oct 9 '14 at 6:45
4

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))
  • 2
    Another approach would be the "gIntersects" function in the rgeos package. – Jeffrey Evans Oct 9 '14 at 18:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.