Optimizing Intersection Area Loops on big layers using PyQGIS

Question

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 = 
QgsProject().instance().mapLayersByName('PHN_boundaries_AUS_May2017_V7')[0]

d = QgsDistanceArea()
d.setEllipsoid('WGS84')

#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 
POA_layer.getFeatures()}
parent_region_features = {feature.id(): feature for (feature) in 
parent_region_layer.getFeatures()}

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

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

Answer 1

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()
d.setEllipsoid('WGS84')

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():
    index.insertFeature(v)

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())
        engine.prepareGeometry()
        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'])
            else:
                print('There is no intersection')
                writer.writerow([str(POA_name)] + [str(parent_region_name)] + ['No intersection- 0'])
output_file.close()