2

I am trying to get the intersection of the following buffers in R:

Circles

I know that for two shapefiles with gIntersects or st_intersection I can get the intersected areas between two shapefiles. However, these belong to the same shapefile, and hence are stored within the same layer. How could I get back ALL of the intersections? For instance, in the image we can see that there are 3 intersecting circles. I would like to get all 4 intersections: 1U2,2U3,1U3 and 1U2U3. Is this possible?

  • Try gIntersection() with byid=TRUE – Loïc Dutrieux Jul 25 '17 at 20:32
  • With PyQGIS code, this is possible by using 'unary_union' and 'polygonize' methods from shapely python module (please, see my answer). – xunilk Jul 25 '17 at 21:53
  • @LoïcDutrieux how? with gIntersection you need two layers, and I have only got one! Repeating the layer does not produce the desired result – Miguel M. Jul 26 '17 at 7:47
1

With PyQGIS code, this is possible by using 'unary_union' and 'polygonize' methods from shapely python module. All intersections are saved as features of a polygon memory layer. Complete code is:

from shapely.wkt import loads
from shapely.geometry import LineString
from shapely.ops import unary_union, polygonize

layer = iface.activeLayer()

polygons = [ feat.geometry().exportToWkt() for feat in layer.getFeatures() ]

shapely_polygons = [ loads(pol) for pol in polygons ]

rings = [ LineString(pol.exterior.coords) for pol in shapely_polygons ]

union = unary_union(rings)

new_intersections = [ geom.wkt for geom in polygonize(union) ]

epsg = layer.crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'new_polygons',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(new_intersections)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromWkt(new_intersections[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

I tried it out with layer of next image:

enter image description here

After running the code at Python Console of QGIS, I got all intersections as expected. They are selected at next image. In attributes table, are pointed out all of them.

enter image description here

  • Worked great! but we still have the issue of matching the interesections to know which of the original polygons interact. We can do this via a spatial match in R or Python – Miguel M. Jul 26 '17 at 13:12
1

I've an unsatisfactory solution using R. Here goes.

First make some sample data - three intersecting circles and a detached circle:

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

pts = SpatialPoints(cbind(c(0,.5,1,2.5),c(0,1,0,2.5)))

circles = gBuffer(pts, width=1.2, byid=TRUE)

plot(circles, border="black", lwd=4)
text(coordinates(circles),as.character(1:4),col="blue",cex=3)

To get all the intersections, convert to lines, merge, then node them, then make new polygons:

merged = gLineMerge(as(circles,"SpatialLines"))
regions = gPolygonize(gNode(merged))

Show the 8 regions thus:

par(mar=c(0,0,0,0))
par(mfrow=c(2,4))
for(i in 1:8){plot(regions);plot(regions[i,],add=TRUE,col="red")}

enter image description here

Getting the relationship between the original circles and the 8 regions is the unsatisfactory bit. I would have thought this would work:

gContains(circles, regions, byid=TRUE)
      1     2     3     4
1 FALSE FALSE FALSE FALSE
2 FALSE FALSE FALSE FALSE
3 FALSE FALSE FALSE FALSE
4 FALSE FALSE FALSE FALSE
5 FALSE  TRUE FALSE FALSE
6 FALSE FALSE FALSE FALSE
7 FALSE FALSE FALSE FALSE
8 FALSE FALSE FALSE  TRUE

But that only shows region 5 overlaps circle 2. The circles and regions with numbers (blue for circles) look like this:

enter image description here

Which shows that region 5 does overlap circle 2, but we also want region 2 to overlap circle 3, and so on. Something isn't right.

The unsatisfactory fix is to check the relation with slightly shrunk regions using a negative buffer:

 rbuffs = gBuffer(regions, width=-0.001, byid=TRUE)
 gContains(circles, rbuffs, byid=TRUE)
      1     2     3     4
1  TRUE FALSE  TRUE FALSE
2 FALSE FALSE  TRUE FALSE
3  TRUE FALSE FALSE FALSE
4  TRUE  TRUE FALSE FALSE
5 FALSE  TRUE FALSE FALSE
6  TRUE  TRUE  TRUE FALSE
7 FALSE  TRUE  TRUE FALSE
8 FALSE FALSE FALSE  TRUE

Now see that circle 1 overlaps regions 1, 3, 4, and 6. Circle 2 overlaps 4,5,6,7, and circle 3 overlaps 1,2,6,7 - you an verify these from the diagram. From that matrix you can then tell which regions come from which circles.

How small should the buffer size be? I think its all to do with the precision of the geometry, so too small and you won't get the effect, too big and you'll might miss some overlaps for small areas. Anyway, until I find a better solution, this might do.

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.