3

I'm using the intersection tool (processing toolbox) to intersect polygons. I have a huge number of small polygons to intersect and it needs so much time to do do that. Maybe it is caused by the great extent of the polygon shape. Is it possible to set a maximum distance which defines the max distance between polygons they will be intersected. I would like to put this option into a python script if it is possible.

7
  • Can you provide a graphical sample of your input layers? Therefore, do you need to necessarily run the intersection algorithm within your script or you only need to intersect features (without considering using the intersection tool)?
    – mgri
    Commented Dec 6, 2016 at 15:50
  • dl.dropboxusercontent.com/u/17606602/example.zip
    – Pimpel
    Commented Dec 6, 2016 at 16:19
  • Are you looking for the intersection between two layers or a self-intersection between all the features from the same layer?
    – mgri
    Commented Dec 6, 2016 at 19:28
  • Oh sorry. Yes, a self intersection I need to find overlapping polygons.
    – Pimpel
    Commented Dec 6, 2016 at 19:29
  • Have you considered PostGIS/Spatialite for this operation? It will be so much quicker - I tried it on your example data and I got 20s vs 3 minutes
    – Liam G
    Commented Dec 7, 2016 at 3:23

1 Answer 1

3

I tried out your shapefile with next script, at the Python Console of QGIS, to produce a memory layer with all intersections (intersection of intersections; 57,401 features) of original shapefile (20,730 features).

import fiona
from shapely.geometry import shape, LineString
from shapely.ops import unary_union, polygonize
import time

c = fiona.open('/home/zeito/pyqgis_data/example/example_1.shp')

collection = []

start = time.time()

print "wait..."

for i,item in enumerate(c):
    geom = shape(item['geometry'])
    collection.append(geom)

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

union = unary_union(rings)

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

epsg = int(c.crs.values()[0].split(':')[1])

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

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

prov = mem_layer.dataProvider()

for i,feat in enumerate(new_intersections):
    feat = QgsFeature()
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromWkt(new_intersections[i]))
    prov.addFeatures([feat])

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

end = time.time()

time_tot = end - start

print "total time = %.4f s" % time_tot 

At the first image, I selected one arbitrary feature in original shapefile:

enter image description here

At second image, where memory layer has a 50 % of transparency, it can be observed that this original feature has seven features in the new layer.

enter image description here

4
  • How long did the processing take with your script?
    – Pimpel
    Commented Dec 7, 2016 at 13:58
  • It took about 18 minutes.
    – xunilk
    Commented Dec 7, 2016 at 14:00
  • Then it is not really an improvement. Actually I wanted to improve the speed of intersection processing by using a max distance.
    – Pimpel
    Commented Dec 7, 2016 at 14:24
  • 2
    Your shapefile has 20,730 features. The number of combinations, taken 2 at a time, is 20,730*20,729/2 = 214,856,085 for evaluation of first intersections. For this reason, the intersection tool of processing toolbox spent a lot of time. In editing mode of your produced layer with this method, you could observe a lot of superposed features in some areas. So, only 18 minutes for producing all 57,401 intersections (between millions of possibilities with intersections of intersections) is really a very good performance and a great improvement.
    – xunilk
    Commented Dec 7, 2016 at 22:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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