Using 3.10.1-A Coruña on Windows 10

I am looping over some (10) national shapefiles and finding intersection and areas but the process is very slow.

I have two loops - one loop for each postcode and one loop for each feature in the Parent polygon layer.

It takes an hour to test 31 Parent polygons against 2,670 postcodes.

It takes 2.25 hours to test 43 Parent polygons against the same number of postcodes.

I have some shapefiles with hundreds and thousands of polygons and I fear this will literally take weeks of computation.

I'm using Ben W's code suggestion from my post PyQGIS layer intersection area seems really inaccurate. Ben W has been a great help to me on this site. I do not want to give the impression that I disrespect Ben W skills and helpful experience at all.

Are there optimization strategies for these kind big loop QGIS tasks?

I'm trying Ben W's suggestion of a spatial index. I've put together this code, but am stuck on the next step:

POA_layer = QgsProject().instance().mapLayersByName('POA_2016_AUST')[0]
parent_region_layer = 

d = QgsDistanceArea()

#POA_features = [f for f in POA_layer.getFeatures()]
#parent_region_features = [f for f in parent_region_layer.getFeatures()]

POA_features = {feature.id(): feature for (feature) in 
parent_region_features = {feature.id(): feature for (feature) in 

POA_index = QgsSpatialIndex()
map(POA_index.insertFeature, POA_features.values())

Parent_index = QgsSpatialIndex()
map(Parent_index.insertFeature, parent_region_features.values())

Try this modified code from your previous question below. I have reworked it a bit- creating a spatial index on the postcode layer and running the outer for loop on the parent region layer and the inner loop on the postcode layer.

One effect of the changes is you will get many fewer of the 'No intersection' rows in the resulting csv file. This is because the intersection test is not being run for every feature in one layer against every feature in the other layer. The inner loop is only testing the post code features which intersect the bounding box of the parent region features on each iteration.

Hopefully this will be a bit faster for you.

Edit after discussion in comments:

I have tweaked the code a bit more to try and achieve some further optimization by making the following changes:

-Using a geometry engine on the parent region layer.

-Storing both the feature ids and features of the postcode layer in a dictionary. This means we can remove a getFeatures() call from the outer loop (since these are expensive and should be minimized where possible) and should make retrieving the features a bit quicker.

import csv
output_file = open('D:\\Folder\\Intersection.csv', 'w', newline='')
writer = csv.writer(output_file)
writer.writerow(['POA Feature'] + ['Parent Region'] + ['Percentage'])

POA_layer = QgsProject().instance().mapLayersByName('POA_2016_AUST')[0]
parent_region_layer = QgsProject().instance().mapLayersByName('PHN_boundaries_AUS_May2017_V7')[0]

d = QgsDistanceArea()

POA_features = {f.id(): f for f in POA_layer.getFeatures()}
parent_region_features = [f for f in parent_region_layer.getFeatures()]

index = QgsSpatialIndex()
for k, v in POA_features.items():

for parent_region_feature in parent_region_features:
        parent_region_name = parent_region_feature['FIRST_PHN_']
        region_geom = parent_region_feature.geometry()
        engine = QgsGeometry.createGeometryEngine(region_geom.constGet())
        candidate_ids = index.intersects(region_geom.boundingBox())
        candidate_features = [v for k, v in POA_features.items() if k in candidate_ids]
        for POA_feature in candidate_features:
            POA_name = POA_feature['POA_NAME16']
            POA_geom = POA_feature.geometry()
            total_area = d.convertAreaMeasurement(d.measureArea(POA_geom), QgsUnitTypes.AreaSquareKilometers)
            if engine.intersects(POA_geom.constGet()):
                intersection = POA_geom.intersection(region_geom)
                intersection_km2 = d.convertAreaMeasurement(d.measureArea(intersection), QgsUnitTypes.AreaSquareKilometers)
                pcnt = (intersection_km2/total_area)*100
                print('Percentage of {} in parent region {}: {}%'.format(POA_name, parent_region_name, pcnt))
                writer.writerow([str(POA_name)] + [str(parent_region_name)] + [str(pcnt)])
            elif engine.contains(POA_geom.constGet()):
                print('{} is fully enclosed by {}'.format(POA_name, parent_region_name))
                writer.writerow([str(POA_name)] + [str(parent_region_name)] + ['Fully enclosed- 100'])
                print('There is no intersection')
                writer.writerow([str(POA_name)] + [str(parent_region_name)] + ['No intersection- 0'])
  • Seriously man, you've been like a guardian angel looking over my shoulder and whispering in my ear for the past few days. This code is 10X faster. Smart idea to flip the index bounding box/ outer loop. Best wishes to you @BenW ! Jan 2 '20 at 7:24
  • Hi @Corpuscular. No worries, you're more than welcome- glad that helped too. You could probably get a further gain in performance by creating a geometry engine on the parent region layer. I think the spatial index makes the biggest difference though. If I get time I will try and update my answer. I would be interested to see how much further difference it would make. Cheers.
    – Ben W
    Jan 2 '20 at 8:22
  • I've just kicked off my biggest problem layer with 2300 features. I think the source of the shapefile has some influence on the speed. "Offical" shapefiles seem to be processing much faster than "less than official" shapefiles - I'm guessing this might have something to do with the rounding of the polygon coordinates - I'm seeing some "0 percent" and "99.9999x percent" results. <edit> :) In the time it took me to write this the problem layer of 2310 polygons just finished in 6 minutes and produced 17,308 results. Thanks man. Jan 2 '20 at 8:56
  • Cool- that's pretty reasonable time. We could probably cut it down a bit more with the geometry engine if needed. For spatial analysis of very large layers I don't think you can beat SQL queries in Postgis for speed, but at this point I don't have much expertise in that area!
    – Ben W
    Jan 2 '20 at 9:34
  • 1
    Hi @BenW thank you for your repeated and extensive help. I can see a gain using the GeomEngine approach. Using SpatialIndex execution UTC time was begin 22:41:26, end 22:48:17 = 7 min 15 sec. Using GeomEngine execution UTC time was begin 22:32:12, end 22:38:39 = 6 min 45 sec. So a 10% gain. :) I'm using the largest datasets in my project with this test but I guess they are not so big as some shapefiles where a 10% gain could save heaps of time. I too hope to one day master Postgis as I have no knowledge of it at the moment. Thanks Ben W. Cheers! Jan 2 '20 at 23:00

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.