Saturday, 14 July 2018

qgis - Finding the common borders between polygons in the same shapefile


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.)



Answer



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)


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)


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}})

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...