4

The question is essentially the same as here, except that I want to do it in qgis. Thus, is there a qgis function equivalent to the ArcMap intersect function that will take a polygon layer and output the line segments corresponding to the intersections (common boundaries) of two polygons. (I realize that since polygons are independent objects in shapefiles, satisfying the intent of the query is more challenging than in a topolically modeled division of the surface.)

  • I came here to file that exact same question, only to find it on top of the newest questions page. Crossing fingers someone will answer. – nirvn May 11 '14 at 6:37
9

The problem is easily solved in Python with the module Shapely, but, as far as I know, a solution does not exist natively in QGIS (only between 2 layers). I will try to explain a solution in the Python console with PyQGIS (sorry if you don't know Python, but I have no other solution, generally I use the module Fiona to do it, simpler and faster and without GIS).

The problem (you want the red line):

enter image description here

is the same as the intersection of the two exterior rings = a LineString:

enter image description here

In the Python console:

# select the active layer
layer = qgis.utils.iface.activeLayer()
# first feature of the layer
elem1 = layer.getFeatures().next()
poly1 = elem1.geometry()
# second feature of the layer
elem2 = layer.getFeatures().next()
poly2 = elem2.geometry()

Now, you can apply the predicates and operations (PyQGIS: Geometry Handling- the correct predicate is geometry1.touches(geometry2)but, I do not know why, I am not getting the same result as Shapely, so I use geometry1.intersects(geometry2). If the the layer is topologically incorrect, you can add a conditionif i[0].intersects(i1).wkbType() == QGis.WKBLineString:`).

With the polygons:

poly1.intersects(poly2)
True
# so the intersection line
poly1.intersection(poly2)
<qgis.core.QgsGeometry object at 0x12c610680>

With the exterior ring:

# compute the rings
ring1= QgsGeometry.fromPolyline(poly1.asPolygon()[0])
ring2= QgsGeometry.fromPolyline(poly2.asPolygon()[0])
# predicate
ring1.intersects(ring2)
True
ring1.intersection(ring2)
<qgis.core.QgsGeometry object at 0x12c610680>

Therefore, with a layer with many polygons,

enter image description here

the solution is to iterate through pairs of geometries in the layer. For that, i will use the itertools module:

import itertools
# list of exterior rings
rings = [QgsGeometry.fromPolyline(elem.geometry().asPolygon()[0]) for elem in layer.getFeatures()]
# compute the intersection lines
for i in  itertools.combinations(rings, 2):
    if i[0].intersects(i[1]):
        i[0].intersection(i[1]).exportToWkt()
u'LINESTRING(106.79212560184775782 -208.05159683595101683, 111.57578217667511922 -224.87687168534384341)'
u'LINESTRING(161.91452099965721345 -249.48826040815123406, 184.32247129585158518 -238.84448401745891033)'
u'LINESTRING(186.5298528721075968 -183.15279895785877784, 163.86119984940140171 -173.64006063583030937)'

and it is easy to create a memory layer with the results. `

With Fiona and Shapely:

import fiona
from shapely.geometry import shape, mapping
Multi = MultiPolygon([shape(poly['geometry']) for poly in fiona.open("polygons.shp")])
import itertools
schema = {'geometry': 'LineString','properties': {'test': 'int'}}
# creation of the new shapefile
with fiona.open('intersection.shp','w','ESRI Shapefile', schema) as e:
    for i in  itertools.combinations(Multi, 2):
       if i[0].touches(i[1]):
           e.write({'geometry':mapping(i[0].intersection(i[1])), 'properties':{'test':1}})
  • Having trouble following - are poly1 and poly2 later written as geom1 and geom2? – Simbamangu May 11 '14 at 16:58

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.