I'm working on a project involving bringing together multiple subjects circling difficult to find parts of an image. The output from users is pretty awesome.

Four classifications

Cool, right?

However, what I want is to get the area of the image that has 1, 2, 3, 4, etc. users selecting it. And I'm using R for the analysis with the sp, raster, and maptools packages (with a little rgeos throw in when need be). Essentially, my workflow is to 1) create a single polygon from each user's input (I actually have the track of each individual selection per user, but it seemed easier to combine them into a single SpatialPolygon per user) 2) merge them all together into a single SpatialPolygonsDataFrame Object 3) Rasterize the object using the count function 4) Assess area of the raster with n number of users selecting it

Here's the output SpatialPolygonsDataFrame to play with

This works great, but is slooow due to the rasterization. I feel like I should be able to do something with the SpatialPolygons object. And I've tried a few things with union, gIntersects, gUnion, and intersect. However, I keep getting errors such as the following

Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, "rgeos_difference") : 
TopologyException: Input geom 1 is invalid: Self-intersection at or near point

Fair enough. Often a given users SpatialPolygons are invalid when I test them using gIsValid. Looking at, say, just the first two, this becomes obvious that there are many overlapping points.

Polygons one and two

You can see the intersecting points pretty clearly. Furthermore, in poly2 (the blue) there appears to be some self-intersection, causing gIsValid to return false.

I've also tried unionSpatialPolygons (with avoidRGEOS=T) which creates a unified object unionSpatialPolygons(polysData[1:2,], IDs=names(polysData[1:2,]@polygons)) But the borders still overlap like in the above image, rather than being one smooth polygon outer border - i.e., the polygons are not dissolved together. Then using any further functions on this new SpatialPolygons object, I have the same problems.

So, is there a way to use polygons instead of going to rasters? It sure would make my life faster, which would be excellent. I feel like this must be a common problem with a wheel invented that I have yet to find. Thoughts?

  • Interesting question ... are you open to other FOSS tools like GRASS or PostGIS? – Simbamangu Jul 11 '15 at 16:23
  • As long as it's scriptable (as I have a huge amount of datasets just like this one), then yes! – jebyrnes Jul 12 '15 at 15:05
  • If the answer I gave below doesn't help, can find a GRASS solution too. – Simbamangu Jul 13 '15 at 6:18
up vote 3 down vote accepted

Sounds like cleaning the geometry up will get you on the road to your non-raster solution ... herewith something of a kludge which does help with fixing bad geometry:

# Load the library and problematic data
library(rgeos)
load("oneImage2_spdf.Rdata")

>gIsValid(polysData, reason = T)
Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, "rgeos_union") : 
TopologyException: Input geom 0 is invalid: Self-intersection at or near point -119.84228271000001 34.349193950783807 at -119.84228271000001 34.349193950783807

That's the offending data; perhaps only a single point, but doing a buffer by 0 width should fix it.

# Let's buffer by 0:
polysData <- gBuffer(polysData, width=0, byid = T)

gBuffer doesn't like this - there's an unprojected CRS attached to it!

Warning message:
In gBuffer(polysData, width = 0, byid = T) :
  Spatial object is not projected; GEOS expects planar coordinates

Time to trick rgeos:

# Set to a projected CRS - could be anything
polysCRS <- polysData@proj4string
polysData@proj4string <- CRS("+init=epsg:32737")

# Now buffer it but preserve individual objects
polysData <- gBuffer(polysData, width = 0, byid = T)

# Reset the original CRS
polysData@proj4string <- polysCRS

Check if geometry's valid:

>gIsValid(polysData, reason = T)
"Valid Geometry"
  • @jebyrnes, did this work for you? – Simbamangu Jul 21 '15 at 5:59
  • Huh. It did! Thanks! Now I need to figure out a bit more about how to tell how many of the polygons overlapped from gUnionCascaded - that heatmap, if you will, but the gBuffer made the geometry valid! – jebyrnes Jul 27 '15 at 15:09
  • 1
    Just a word of caution, The gBuffer function 'fixes' the problem by removing the self-intersecting polygons. To visualize what I mean before running gBuffer, create a copy of polysData: polysData_original <- polysData then use gBuffer as above and then plot both sets of polygons: for(i in seq_along(polysData)){ plot(polysData_original[i,],col='red'); plot(polysData[i,],col='green',add=T); title(i); } – Remi Daigle Aug 3 '15 at 17:04

Without your original data, I can't be sure this will work, but I thought it might help you out. I didn't bring it all the way there, this solution still likely needs some level of automation, but might give you a general way forward

First, I create some spatial polygons

polypoints1 <- matrix(c(1,2,2,1,1,2,2,1,1,2),ncol=2)
polypoints2 <- matrix(c(1,3,3,1,1,3,3,1,1,3),ncol=2)
polypoints3 <- matrix(c(1,2,2,1,1,2,2,1,1,2)+1.1,ncol=2)
polypoints4 <- matrix(c(1,2,2,1,1,2,2,1,1,2)+0.5,ncol=2)

p1 <- Polygon(polypoints1)
ps1 <- Polygons(list(p1),1)
sps1 <- SpatialPolygons(list(ps1))

p2 <- Polygon(polypoints2)
ps2 <- Polygons(list(p2),2)
sps2 <- SpatialPolygons(list(ps2))

p3 <- Polygon(polypoints3)
ps3 <- Polygons(list(p3),3)
sps3 <- SpatialPolygons(list(ps3))

p4 <- Polygon(polypoints4)
ps4 <- Polygons(list(p4),4)
sps4 <- SpatialPolygons(list(ps4))

I plotted them just to see

plot(sps2,col='green')
plot(sps1,add=T,col='blue')
plot(sps3,add=T,col='yellow')
plot(sps4,add=T,col='purple')

I merged them into an spdf

data=data.frame(c(x=rep(1,4)),row.names=c(1:4))
sps <- SpatialPolygons(list(ps1,ps2,ps3,ps4))
spdf <- SpatialPolygonsDataFrame(sps,data)

You can identify which polygon overlaps which like so:

gIntersects(spdf,spdf,byid =T)

From the above command you could create some kind of loops to do the overlapping combinations below (I'm just ignoring sps4 for brevity at this point)

poly2a <- gIntersection(spdf[2,],spdf[1,],drop_lower_td=T)
poly2a <- SpatialPolygonsDataFrame(poly2a,data.frame(c(x=1),row.names=c(1)))
plot(poly2a,add=T,col='red')

This time we need to change the ID since we're going to rbind these later

poly2b <- gIntersection(spdf[2,],spdf[3,],drop_lower_td=T)
poly2b <- spChFIDs(poly2b,"2")
poly2b <- SpatialPolygonsDataFrame(poly2b,data.frame(c(x=1),row.names=c(2)))
plot(poly2b,add=T,col='red')

Merge the overlapping polygons into another spdf

spdf_overlaps <- rbind(poly2a,poly2b)
poly2 <- unionSpatialPolygons(spdf_overlaps,rep(1,2))
plot(poly2,add=T,col='blue')

Now we have poly2 which is where we have 2 layers overlapping (except combinations with sps4) then to figure out 3 layers, we just have to check out where poly2 and spdf overlap (if you make a more automated version of this, you'll need to make sure that 'poly2' does not include sps4 as in this example)

gIntersects(poly2,spdf[4,],byid =T)

poly3 <- gIntersection(poly2,spdf[4,],drop_lower_td=T)
plot(poly3,add=T,col="red") 

Check it out

gIsValid(poly2)
gIsValid(poly3)

Alternatively, you could always do a pseudo rasterization, much easier, but you loose some detail depending on your cell size:

First make the grid:

bb <- bbox(spdf)
cs <- c(0.1,0.1)  # cell size
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells per direction
grd <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)


sp_grd <- SpatialGridDataFrame(grd,
                           data=data.frame(id=1:prod(cd)))

Then, make grid into a polygon which used for overlap

library(Grid2Polygons)
grid <- Grid2Polygons(sp_grd)
plot(grid)

Then count the number of polygons that overlap each grid cell

count <- apply(gContains(spdf,grid,byid=T),1,sum)

Finally, plot it!

plot(grid)
for(i in 1:length(grid)){
    plot(grid[i,],col=rev(heat.colors(3))[count[i]],add=T)
}
  • FYI - github.com/jebyrnes/floatingForests/blob/master/data/… is the spdf from the above example. – jebyrnes Jul 10 '15 at 12:19
  • The gIntersects method still fails with the above data. I like the pseudo-rasterization idea. Indeed, the original polygons can be made on the scale of the image - so, 1 pixel = 1 point, and then the ultimate overlapping polygon can be projected for area calculation. This would take a wee bit of doing with respect to my workflow, but, if it's the only solution... – jebyrnes Jul 10 '15 at 12:31
  • Ahhh, I see, the original data has self intersecting polygons, I think the only purely R solution will be to use the pseudo-rasterized method. As far as I know, R can't handle cleaning those up. But they could be cleaned up in GRASS, or PostGIS, or begrudgingly in ArcGIS Check out: gis.stackexchange.com/questions/141179/… stat.ethz.ch/pipermail/r-sig-geo/2013-January/017273.html – Remi Daigle Jul 10 '15 at 21:26
  • I just found out R and GRASS can interact! grasswiki.osgeo.org/wiki/R_statistics If you figure out how to fix your polygons in GRASS, you should be able to automate it in R – Remi Daigle Jul 21 '15 at 5:36

Applying a buffer with width=0 can be risky for some polygons, e.g. cases of bowtie polygons, altering a lot source geometries. You can try to use the cleangeo package which aims to fix spatial objects. You can install it from Github or from CRAN. You can use this simple code to correct your spatial object:

library(cleangeo)
polysData.clean <- clgeo_Clean(polysData)

#check geometry validity of output
gIsValid(polysData.clean)

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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