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)?– mgriDec 6, 2016 at 15:50
-
dl.dropboxusercontent.com/u/17606602/example.zip– PimpelDec 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?– mgriDec 6, 2016 at 19:28
-
Oh sorry. Yes, a self intersection I need to find overlapping polygons.– PimpelDec 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 GDec 7, 2016 at 3:23
1 Answer
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:
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.
-
-
-
Then it is not really an improvement. Actually I wanted to improve the speed of intersection processing by using a max distance.– PimpelDec 7, 2016 at 14:24
-
2Your 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.– xunilkDec 7, 2016 at 22:13