Finding nearest polygon in PyQGIS using QgsSpatialIndex

Question

I am trying to write a custom processing tool for QGIS which joins the nearest attributes with some conditions. Therefore I'd like to use a QgsSpatialIndex to increase performance. This works perfectly for single points or centroids. But since QgsSpatialIndex.nearestNeighbor() expects a QgsPointXY as input this wont work properly to find the nearest multipoint, polygon or line. To focus this question I am only asking for polygons here. A polygon inside another polygon would then of course be the closest polygon, otherwise, lets say the shortest distance from its outer ring to outer ring.

To provide some shortest code:

from datetime import datetime
starttime = datetime.now()
source_layer = QgsProject.instance().mapLayersByName('sourcelayer_polygon')[0]
join_layer = QgsProject.instance().mapLayersByName('joinlayer_polygon')[0]

join_layer_idx = QgsSpatialIndex(join_layer.getFeatures())

for current, source_feat in enumerate(source_layer.getFeatures()):
    source_feat_centroid = source_feat.geometry().centroid()
    nearest_neighbors = join_layer_idx.nearestNeighbor(source_feat_centroid, neighbors = -1, maxDistance = 1234)
    for join_feat_id in nearest_neighbors:
        join_feat = join_layer.getFeature(join_feat_id)
        if source_feat["name"] == join_feat["name"]:
            print('Source feat: ' + str(source_feat["id"]) + ' - ' + str(source_feat["name"]))
            print('Join feat: ' + str(join_feat["id"]) + ' - ' + str(join_feat["name"]))
            print('Nearest join feat dist: ' + str(join_feat.geometry().distance(source_feat_centroid)))
            print('\n')
            
print('Runtime: ' + str((datetime.now() - starttime).total_seconds()))

I have tried to learn from the processing tool "join attributes by nearest" here: https://github.com/qgis/QGIS/commit/95af4d4a45f48889358446a0c9d02eb9d54860d7#diff-86d605b8b46166ffcb11eee1c453712517f25712fe1166b3b6659963ccfab2fb But so far I was not able to figure how it uses the index exactly.

Also, of course I could build a dictionary with distances from boundary to boundary or similar, but I am afraid this would be extremely performance intensive.

How can I efficiently find the nearest polygons using PyQGIS?

For reference here to full code of my processing script:

import operator, processing
from PyQt5.QtCore import QCoreApplication, QVariant
from qgis.core import (QgsField, QgsFeature, QgsProcessing, QgsExpression, QgsSpatialIndex,
                       QgsFeatureSink, QgsFeatureRequest, QgsProcessingAlgorithm,
                       QgsProcessingParameterFeatureSink, QgsProcessingParameterField, QgsProcessingParameterFeatureSource, QgsProcessingParameterEnum, QgsProcessingParameterExpression, QgsProcessingParameterNumber, QgsProcessingParameterString)

class JoinAttributesByNearestCentroidWithCondition(QgsProcessingAlgorithm):
    SOURCE_LYR = 'SOURCE_LYR'
    SOURCE_FIELD = 'SOURCE_FIELD'
    SOURCE_EXPRESSION = 'SOURCE_EXPRESSION'
    JOIN_LYR = 'JOIN_LYR'
    JOIN_FIELDS = 'JOIN_FIELDS'
    JOIN_FIELD = 'JOIN_FIELD'
    JOIN_EXPRESSION = 'JOIN_EXPRESSION'
    OPERATION = 'OPERATION'
    JOIN_N = 'JOIN_N'
    JOIN_DIST = 'JOIN_DIST'
    JOIN_PREFIX = 'JOIN_PREFIX'
    OUTPUT = 'OUTPUT'

    def initAlgorithm(self, config=None):
        
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.SOURCE_LYR, self.tr('Source Layer')))
        self.addParameter(
            QgsProcessingParameterExpression(
                self.SOURCE_EXPRESSION, self.tr('Filter-Expression for Source-Layer'), parentLayerParameterName = 'SOURCE_LYR', optional = True))
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.JOIN_LYR, self.tr('Join Layer')))
        self.addParameter(
            QgsProcessingParameterExpression(
                self.JOIN_EXPRESSION, self.tr('Filter-Expression for Join-Layer'), parentLayerParameterName = 'JOIN_LYR', optional = True))
        self.addParameter(
            QgsProcessingParameterField(
                self.JOIN_FIELDS, self.tr('Add the following fields of join layer to result (if none are chosen, all fields will be added)'),parentLayerParameterName='JOIN_LYR', allowMultiple = True, optional = True))
        self.addParameter(
            QgsProcessingParameterNumber(
                self.JOIN_N, self.tr('Join x nearest neighbors'), type = 0, defaultValue = 3, minValue = 1, maxValue = 2147483647))
        self.addParameter(
            QgsProcessingParameterNumber(
                self.JOIN_DIST, self.tr('Maximum join distance (0 means unlimited)'), type = 1, defaultValue = 0, minValue = 0, maxValue = 2147483647))
        self.addParameter(
            QgsProcessingParameterString(
                self.JOIN_PREFIX, self.tr('Join Prefix'), defaultValue = 'join_'))
        self.addParameter(
            QgsProcessingParameterField(
                self.SOURCE_FIELD, self.tr('Source Layer compare Field'),parentLayerParameterName='SOURCE_LYR', optional = True))
        self.addParameter(
            QgsProcessingParameterEnum(
                self.OPERATION, self.tr('Comparison operator (if no operator is set, the comparison fields remain unused)'), [None,'!=','=','<','>','<=','>=','is','not','is not'], defaultValue = 0, allowMultiple = False))
        self.addParameter(
            QgsProcessingParameterField(
                self.JOIN_FIELD, self.tr('Join Layer compare Field'),parentLayerParameterName='JOIN_LYR', optional = True))
        self.addParameter(
            QgsProcessingParameterFeatureSink(
                self.OUTPUT, self.tr('Joined Layer')))

    def processAlgorithm(self, parameters, context, feedback):
        # Get Parameters
        source_layer = self.parameterAsSource(parameters, self.SOURCE_LYR, context)
        source_field = self.parameterAsFields(parameters, self.SOURCE_FIELD, context)
        if source_field:
            source_field = source_field[0]
        join_layer = self.parameterAsLayer(parameters, self.JOIN_LYR, context)
        join_field = self.parameterAsFields(parameters, self.SOURCE_FIELD, context)
        if join_field:
            join_field = join_field[0]
        join_fields = self.parameterAsFields(parameters, self.JOIN_FIELDS, context)
        operation = self.parameterAsInt(parameters, self.OPERATION, context)
        ops = { # get the operator by this index
            0: None,
            1: operator.ne,
            2: operator.eq,
            3: operator.lt,
            4: operator.gt,
            5: operator.le,
            6: operator.ge,
            7: operator.is_,
            8: operator.not_,
            9: operator.is_not
            }
        op = ops[operation]
        source_expression = self.parameterAsExpression(parameters, self.SOURCE_EXPRESSION, context)
        source_expression = QgsExpression(source_expression)
        join_expression = self.parameterAsExpression(parameters, self.JOIN_EXPRESSION, context)
        join_expression = QgsExpression(join_expression)
        join_n = self.parameterAsInt(parameters, self.JOIN_N, context)
        join_dist = self.parameterAsDouble(parameters, self.JOIN_DIST, context)
        join_prefix = self.parameterAsString(parameters, self.JOIN_PREFIX, context)


        source_layer_fields = source_layer.fields()
        if join_fields:
            join_layer = join_layer.materialize(QgsFeatureRequest().setSubsetOfAttributes(join_fields, join_layer.fields()))
        join_layer_fields = join_layer.fields()
        output_layer_fields = source_layer_fields
        for join_layer_field in join_layer.fields():
            output_layer_fields.append(QgsField(join_prefix + join_layer_field.name(), join_layer_field.type()))
        output_layer_fields.append(QgsField(join_prefix + 'dist', QVariant.Double, len=20, prec=5))
        
        (sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
                                               output_layer_fields, source_layer.wkbType(),
                                               source_layer.sourceCrs())
        
        if source_expression not in (QgsExpression(''),QgsExpression(None)):
            source_layer = source_layer.materialize(QgsFeatureRequest(source_expression))
        if join_expression not in (QgsExpression(''),QgsExpression(None)):
            join_layer = join_layer.materialize(QgsFeatureRequest(join_expression))
        
        total = 100.0 / source_layer.featureCount() if source_layer.featureCount() else 0
        
        if source_layer.sourceCrs() != join_layer.sourceCrs():
            reproject_params = {'INPUT': join_layer, 'TARGET_CRS': source_layer.sourceCrs(), 'OUTPUT': 'memory:Reprojected'}
            reproject_result = processing.run('native:reprojectlayer', reproject_params)
            join_layer = reproject_result['OUTPUT']
            
        join_layer_idx = QgsSpatialIndex(join_layer.getFeatures())
        
        for current, source_feat in enumerate(source_layer.getFeatures()):
            if feedback.isCanceled():
                break
            matches_found_counter = 0
            source_feat_centroid = source_feat.geometry().centroid()
            
            if op is None:
                nearest_neighbors = join_layer_idx.nearestNeighbor(source_feat_centroid, neighbors = join_n, maxDistance = join_dist)
            else:
                nearest_neighbors = join_layer_idx.nearestNeighbor(source_feat_centroid, neighbors = -1, maxDistance = join_dist)
            
            for join_feat_id in nearest_neighbors:
                if matches_found_counter >= join_n:
                    break
                join_feat = join_layer.getFeature(join_feat_id)
                if op is None:
                    matches_found_counter += 1
                    new_feat = QgsFeature(output_layer_fields)
                    new_feat.setGeometry(source_feat.geometry())
                    attridx = 0
                    for attr in source_feat.attributes():
                        new_feat[attridx] = attr
                        attridx += 1
                    for attr in join_feat.attributes():
                        new_feat[attridx] = attr
                        attridx += 1
                    new_feat[join_prefix + 'dist'] = source_feat_centroid.distance(join_feat.geometry())
                    sink.addFeature(new_feat, QgsFeatureSink.FastInsert)
                elif op(source_feat[source_field], join_feat[join_field]):
                    matches_found_counter += 1
                    new_feat = QgsFeature(output_layer_fields)
                    new_feat.setGeometry(source_feat.geometry())
                    attridx = 0
                    for attr in source_feat.attributes():
                        new_feat[attridx] = attr
                        attridx += 1
                    for attr in join_feat.attributes():
                        new_feat[attridx] = attr
                        attridx += 1
                    new_feat[join_prefix + 'dist'] = source_feat_centroid.distance(join_feat.geometry())
                    sink.addFeature(new_feat, QgsFeatureSink.FastInsert)
                
            if matches_found_counter == 0:
                new_feat = QgsFeature(output_layer_fields)
                new_feat.setGeometry(source_feat.geometry())
                attridx = 0
                for attr in source_feat.attributes():
                    new_feat[attridx] = attr
                    attridx += 1
                sink.addFeature(new_feat, QgsFeatureSink.FastInsert)
                    
            feedback.setProgress(int(current * total))
            

        return {self.OUTPUT: dest_id}


    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return JoinAttributesByNearestCentroidWithCondition()

    def name(self):
        return 'JoinAttributesByNearestCentroidWithCondition'

    def displayName(self):
        return self.tr('Join attributes by nearest centroid with condition')

    def group(self):
        return self.tr('FROM GISSE')

    def groupId(self):
        return 'from_gisse'

    def shortHelpString(self):
        return self.tr('This Algorithm finds the x nearest neighbors (centroids) by a given condition and joins them.')

Answer 1

According to the docs, you can also pass a QgsGeometry object to the QgsSpatialIndex.nearestNeighbor() method, which means it can be used on line and polygon geometries etc. You should also construct the spatial index using the FlagStoreFeatureGeometries flag.

If this QgsSpatialIndex object was not constructed with the FlagStoreFeatureGeometries flag, then the nearest neighbor test is performed based on the feature bounding boxes ONLY, so for non-point geometry features this method is not guaranteed to return the actual closest neighbors.

I tested with a polygon layer and a line layer and made a small script which calculates nearest neighbor/s in the join layer for every feature in the source layer.

Finally, it selects a random source feature and its nearest neighbor/neighbors in the join layer, then zooms to the selected features so that one may visually inspect the results.

In terms of timing, the runtime was consistently around 0.8 seconds with around 3300 features in the source layer and 27000 features in the join layer.

from datetime import datetime
from random import randrange

starttime = datetime.now()
source_layer = QgsProject.instance().mapLayersByName('Random_buildings')[0]
join_layer = QgsProject.instance().mapLayersByName('AUS_roads')[0]

join_layer_idx = QgsSpatialIndex(join_layer.getFeatures(), flags=QgsSpatialIndex.FlagStoreFeatureGeometries)

rand_id = randrange(source_layer.featureCount())

results = {}

for source_feat in source_layer.getFeatures():
    nearest_neighbors = join_layer_idx.nearestNeighbor(source_feat.geometry())
    results[source_feat.id()] = nearest_neighbors
            
print('Runtime: ' + str((datetime.now() - starttime).total_seconds()))

source_layer.selectByIds([rand_id])
join_layer.selectByIds(results[rand_id])
iface.mapCanvas().zoomToSelected([source_layer, join_layer])

Example result:

Also, as minor side point- I think that if you use print statements in a for loop (unless your print calls are connected to signals emitted from a background thread) it can bog down the execution and increase your runtime. It is not really noticeable on small datasets, but if you have thousands of features it really slows down.

---

In case of your shortest code it could be:

from datetime import datetime
starttime = datetime.now()
source_layer = QgsProject.instance().mapLayersByName('sourcelayer_polygon')[0]
join_layer = QgsProject.instance().mapLayersByName('joinlayer_polygon')[0]

join_layer_idx = QgsSpatialIndex(join_layer.getFeatures(), flags=QgsSpatialIndex.FlagStoreFeatureGeometries)

for source_feat in source_layer.getFeatures():
    source_feat_geom = source_feat.geometry()
    nearest_neighbors = join_layer_idx.nearestNeighbor(source_feat_geom, neighbors = -1, maxDistance = 1234)
    for join_feat_id in nearest_neighbors:
        join_feat = join_layer.getFeature(join_feat_id)
        if source_feat["name"] == join_feat["name"]:
            print(f'Source feat: {source_feat["id"]} - {source_feat["name"]}')
            print(f'Join feat: {join_feat["id"]} - {join_feat["name"]}')
            print(f'Nearest join feat dist: {join_feat.geometry().distance(source_feat_geom)}\n')

print(f'Runtime: {(datetime.now() - starttime).total_seconds()}')

Answer 2

For now I solved it by building my own dictionary, sorted by distances. As stated in the docs QgsGeometry.distance() returns the minimum distance between this geometry and another geometry, using GEOS. So that should do what I am after. Even though the attribute filter is done in the request and saving a lot of iterations for my usecase, performance is a lot slower than using the spatial index, so I am still interested in faster methods.

from datetime import datetime
starttime = datetime.now()
source_layer = QgsProject.instance().mapLayersByName('sourcelayer_polygon')[0]
join_layer = QgsProject.instance().mapLayersByName('joinlayer_polygon')[0]

for current, source_feat in enumerate(source_layer.getFeatures()):
    request = QgsFeatureRequest().setFilterExpression("\"name\" = '" + source_feat["name"] + "'")
    distances = {join_feat: source_feat.geometry().distance(join_feat.geometry())
                 for join_feat in join_layer.getFeatures(request)}
    distances_sorted = dict(sorted(distances.items(), key=lambda item: item[1]))
    for join_feat in distances_sorted:
        print('Source feat: ' + str(source_feat["id"]) + ' - ' + str(source_feat["name"]))
        print('Join feat: ' + str(join_feat["id"]) + ' - ' + str(join_feat["name"]))
        print('Nearest join feat dist: ' + str(join_feat.geometry().distance(source_feat_centroid)))
        print('\n')
            
print('Runtime: ' + str((datetime.now() - starttime).total_seconds()))