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.

  • 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 Dec 6 '16 at 15:50
  • Are you looking for the intersection between two layers or a self-intersection between all the features from the same layer? – mgri Dec 6 '16 at 19:28
  • Oh sorry. Yes, a self intersection I need to find overlapping polygons. – Pimpel Dec 6 '16 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 Dec 7 '16 at 3:23
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

  • How long did the processing take with your script? – Pimpel Dec 7 '16 at 13:58
  • It took about 18 minutes. – xunilk Dec 7 '16 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 Dec 7 '16 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 Dec 7 '16 at 22:13

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.