Finding nearest neighbor with smaller attribute value in QGIS for a huge dataset: make expression more efficient

Here is a more universal processing tool for this task with some advantages over my basic script:

• Its responsive while executing it
• Ability to cancel the execution
• Usable in graphical modeler
• It creates a copy of the layer instead of editing it
• Simple GUI with progress bar
• Adds ID and attribute of the near neighbor as well as the distance between feature and matching neighbor
• Free choice of fields and operators
• Setting for maximum nearest neighbors to compare (most likely increases performance a lot on large layers but most likely also leads to more non-matches!)
• Setting for maximum search distance (can increase performance a lot on large layers but can also lead to more non-matches!)
• Setting for a value where no comparison is done (can increase performance a lot on large layers, if chosen properly. Typically this is the minimum value of the available attributes when using < or <= operator and the maximum value of the available attributes when using the > or >= operator; only usable when comparing numerical values!)

# V1.3.2

from PyQt5.QtCore import QCoreApplication, QVariant
from qgis.core import (QgsSpatialIndex, QgsProcessingParameterFeatureSink, QgsFeatureSink, QgsField, QgsFields, QgsFeature, QgsGeometry, QgsPoint, QgsWkbTypes, 
                       QgsProcessingAlgorithm, QgsProcessingParameterField, QgsProcessingParameterBoolean, QgsProcessingParameterVectorLayer, QgsProcessingOutputVectorLayer, QgsProcessingParameterEnum, QgsProcessingParameterNumber)
import operator

class NearNeighborAttributeByAttributeComparison(QgsProcessingAlgorithm):
    SOURCE_LYR = 'SOURCE_LYR'
    SOURCE_FIELD = 'ID_FIELD'
    ATTRIBUTE_FIELD = 'ATTRIBUTE_FIELD'
    MAX_NEIGHBORS = 'MAX_NEIGHBORS'
    MAX_DISTANCE = 'MAX_DISTANCE'
    DONOT_COMPARE_VALUE = 'DONOT_COMPARE_VALUE'
    DONOT_COMPARE_BOOL = 'DONOT_COMPARE_BOOL'
    OPERATOR = 'OPERATOR'
    OUTPUT = 'OUTPUT'

    def initAlgorithm(self, config=None):
        
        self.addParameter(
            QgsProcessingParameterVectorLayer(
                self.SOURCE_LYR, self.tr('Source Layer'))) # Take any source layer
        self.addParameter(
            QgsProcessingParameterField(
                self.SOURCE_FIELD, self.tr('Attribute field containing unique IDs'),'id','SOURCE_LYR'))
        self.addParameter(
            QgsProcessingParameterField(
                self.ATTRIBUTE_FIELD, self.tr('Attribute field for comparison'),'year','SOURCE_LYR'))
        self.addParameter(
            QgsProcessingParameterNumber(
                self.MAX_NEIGHBORS, self.tr('Maximum number of nearest neighbors to compare (use -1 to compare all features of the layer)'),defaultValue=1000,minValue=-1))
        self.addParameter(
            QgsProcessingParameterNumber(
                self.MAX_DISTANCE, self.tr('Only take nearest neighbors within this maximum distance into account for comparison'),defaultValue=10000,minValue=0))
        self.addParameter(
            QgsProcessingParameterEnum(
                self.OPERATOR, self.tr('Operator to compare the attribute value (If attribute is of type string, only == and != do work)'),
                    ['<','<=','==','!=','>=','>'],defaultValue=0))
        self.addParameter(
            QgsProcessingParameterNumber(
                self.DONOT_COMPARE_VALUE, self.tr('Do not search for matches on features having a value (insert chosen operator here) x n Only works for numerical values and dependent on the chosen operator. n Typically this should be the max or min value available in the attributes and therefore there cant be a match.'),defaultValue=0))
        self.addParameter(
            QgsProcessingParameterBoolean(
                self.DONOT_COMPARE_BOOL,self.tr('Check this Box to actually use the previous option ("Do not search for matches on ....")'),defaultValue=0))
        self.addParameter(
            QgsProcessingParameterFeatureSink(
                self.OUTPUT, self.tr('Near Neighbor Attributes'))) # Output

    def processAlgorithm(self, parameters, context, feedback):
        # Get Parameters and assign to variable to work with
        layer = self.parameterAsLayer(parameters, self.SOURCE_LYR, context)
        idfield = self.parameterAsString(parameters, self.SOURCE_FIELD, context)
        idfield_index = layer.fields().indexFromName(idfield) # get the fieldindex of the id field
        idfield_type = layer.fields()[idfield_index].type() # get the fieldtype of this field
        attrfield = self.parameterAsString(parameters, self.ATTRIBUTE_FIELD, context)
        attrfield_index = layer.fields().indexFromName(attrfield) # get the fieldindex of the attribute field
        attrfield_type = layer.fields()[attrfield_index].type() # get the fieldtype of this field
        maxneighbors = self.parameterAsDouble(parameters, self.MAX_NEIGHBORS, context)
        maxdistance = self.parameterAsDouble(parameters, self.MAX_DISTANCE, context)
        donotcomparevalue = self.parameterAsDouble(parameters, self.DONOT_COMPARE_VALUE, context)
        donotcomparebool = self.parameterAsBool(parameters, self.DONOT_COMPARE_BOOL, context)
        op = self.parameterAsString(parameters, self.OPERATOR, context)
        op = int(op[0]) # get the index of the chosen operator
        #import operator
        ops = { # get the operator by this index
            0: operator.lt,
            1: operator.le,
            2: operator.eq,
            3: operator.ne,
            4: operator.ge,
            5: operator.gt
            }
        op_func = ops[op] # create the operator function
        
        total = 100.0 / layer.featureCount() if layer.featureCount() else 0 # Initialize progress for progressbar
        
        # if -1 has been chosen for maximum features to compare, use the amount of features of the layer, else use the given input
        if maxneighbors == -1:
            maxneighbors = layer.featureCount()
        
        fields = layer.fields() # get all fields of the inputlayer
        fields.append(QgsField("near_id", idfield_type)) # create new field with same type as the inputfield
        fields.append(QgsField("near_attr", attrfield_type)) # same here for the attribute field
        fields.append(QgsField("near_dist", QVariant.Double, len=20, prec=5)) # add a new field of type double
        
        idx = QgsSpatialIndex(layer.getFeatures()) # create a spatial index

        (sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
                                               fields, layer.wkbType(),
                                               layer.sourceCrs())
                                               
        for current, feat in enumerate(layer.getFeatures()): # iterate over source 
            new_feat = QgsFeature(fields) # copy source fields + appended
            attridx = 0 # reset attribute fieldindex
            for attr in feat.attributes(): # iterate over attributes of source layer for the current feature
                new_feat[attridx] = attr # copy attribute values over to the new layer
                attridx += 1 # go to the next field
            new_feat.setGeometry(feat.geometry()) # copy over the geometry of the source feature
            if ((not(op_func(feat[attrfield], donotcomparevalue))) or (not donotcomparebool)): # only search for matches if not beeing told to not do to so
                nearestneighbors = idx.nearestNeighbor(feat.geometry(), neighbors=maxneighbors, maxDistance=maxdistance) # get the featureids of the maximum specified number of near neighbors within a maximum distance
                try:
                    nearestneighbors.remove(feat.id()) # remove the current feature from this list (otherwise the nearest feature by == operator would always be itself...)
                except:
                    pass # ignore on error
                for near in nearestneighbors: # for each feature iterate over the nearest ones (the index is already sorted by distance, so the first match will be the nearest match)
                    if op_func(layer.getFeature(near)[attrfield], feat[attrfield]): # if the current nearest attribute is (chosen operator here) than the current feature ones, then
                        new_feat['near_id'] = layer.getFeature(near)[idfield] # get the near matchs's id value and fill the current feature with its value
                        new_feat['near_attr'] = layer.getFeature(near)[attrfield] # also get the attribute value of this near feature
                        new_feat['near_dist'] = feat.geometry().distance(layer.getFeature(near).geometry()) # and finally calculate the distance between the current feature and the nearest matching feature
                        break # break the for loop of near features and continue with the next feat
            else: # do not search for near neighbor matches if given value is (operator here) than x
                pass # do nothing and continue adding the feature
                    
            sink.addFeature(new_feat, QgsFeatureSink.FastInsert) # add feature to the output
            feedback.setProgress(int(current * total)) # Set Progress in Progressbar
            
            if feedback.isCanceled(): # Cancel algorithm if button is pressed
                break

        return {self.OUTPUT: dest_id} # Return result of algorithm

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

    def createInstance(self):
        return NearNeighborAttributeByAttributeComparison()

    def name(self):
        return 'NearNeighborAttributeByAttributeComparison'

    def displayName(self):
        return self.tr('Add near neighbor attribute by comparing attribute values')

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

    def groupId(self):
        return 'from_gisse'

    def shortHelpString(self):
        return self.tr(
        'This Algorithm searches the n'
        '- x nearest neighbors n'
        '- within a given maximum distance n'
        'of the current feature and compares a given attribute. n'
        'If this comparison returns true, it adds the id, and the attribute of this neighbor to the current feature as well as the distance to this neighbor. n n '
        'Further explanations available on https://gis.stackexchange.com/a/396856/107424'
        )

To use it, just copy paste the code without any changes and place the python file in C:UsersyourusernameAppDataRoamingQGISQGIS3profilesdefaultprocessingscripts. See the docs for more informations. Once saved, you will find it in your processing toolbox within "Scripts" -> "FROM GISSE"

Only tested in QGIS 3.18.2. I can already tell this wont work in QGIS 3.2 (as intended), as this version has no option to limit the maximum search distance. Don't know when this feature was introduced though. However, feel free to just test in your version, worst thing that can happen is an error or a crash. In this case, upgrade your QGIS or edit the script accordingly.

---

Some further explanations on this script and the used method:

First, a spatial index is built on the layer. Then we iterate over the whole layer. For each "actual feature " (naming it this way to not confuse with my used term "near feature") we are getting the nearest neighbors of the current actual feature by using the QgsSpatialIndex().nearestNeighbor(<current_point>,<maximum_neighbors>,<maximum_distance>) method. This method returns an, by distance ascendingly, ordered list of the nearest feature id's. Then, to prevent iterating over itself and therefore finding itself when using the == operator, we remove the current actual feature's id from this list.

Now we iterate of this ordered list of near feature id's (if there is no given limitation by using the "Do not search for matches on..." option). So first we get the very closest other feature by its id. Check if it fulfills our given attribute comparison criteria. If not, we go to the second nearest feature. Check again. Go to the third nearest... and so on. If we finally find the first one fulfilling the criteria, we are grabbing the two attributes by a focused layer.getFeature(<near_matching_feat_id>)[<field>] request and calculate the distance by using <current_actual_feat>.geometry().distance(layer.getFeature(<near_matching_feat_id>).geometry()). These two operations should not cost anything worth mentioning the calculation time. Also these things are always be done only once for an actual feature as we do not iterate here. As soon as we have done this, we stop (break command) the iteration over the ordered list of near features and go to the next actual feature.

---

Runtime and optimal settings:

So if I am not completely mistaken, the calculation time in the best case should be: t = time_needed_for((number_of_actual_features-1)*1) + time_needed_for(creating_index). The best case is, if for all actual features the very closest other (near) feature fulfills the given attribute comparison criteria. In the worst case it should be t = time_needed_for((number_of_actual_features-1)^2) + time_needed_for(creating_index). The worst case is, if for all actual features the very farthest other (near) feature fulfills the given attribute comparison criteria. So we have three options to prevent the worst case from happening: limit the maximum number of near features to compare (when limiting to max. 1000 comparisons the worst case is t = time_needed_for((number_of_actual_features-1)*1000) + time_needed_for(creating_index) if there are at least 1000 actual features) or less helpful, but still, limit the maximum distance of near features to compare. The second option is less helpful because in this distance, theoretically, still all other near features could be included. As third option we can limit the search by giving an attribute. This makes sense when using the available minimum attribute here when using < or <= operator. Or when using the available maximum attribute when using > or >= operator. This option prevents from searching for a match if the criterion is fulfilled. So e.g. when using < and the minimum year available, prevent the algorithm from searching where there cannot be a match anyway in that case.

So the optimal settings strongly depend on your individual layer and your individual desired result. On the individual layer, because the runtime is dependent on the spatial and attributional distribution of your features. The algorithm searches for the nearest feature, matching the attribute condition, ordered by their distance. If everytime the very nearest feature already is a match, the runtime is short. If not, the more features need to be compared, the greater the runtime is. Here comes your personal desired result in play: you can limit the number of features to compare and/or their maximum distance. Lets say you limit the number of features to compare to 100: If there is no attributional match within the 100 closest features, the loop will be skipped and the current actual feature will simply get no result. So no result, but increased runtime. But if there is a match at lets say the 51st closest feature, this setting does have absolutely no effect, as the loop is skipped anyway after the 51st comparison.

Here an example I have tested; given 1000 randomly distributed points within a given extent of 10km*10km and random years between 1850 and 2020. Not using attribute search limitations in this test.

• Max. Features to compare: 10; Max. searchdistance: 10000m; Runtime = 0.2 seconds; No match for 103 features
• Max. Features to compare: 100; Max. searchdistance: 10000m; Runtime = 0.34 seconds; No match for 11 features
• Max. Features to compare: 1000; Max. searchdistance: 10000m; Runtime = 1.13 seconds; No match for 7 features (all of them have the lowest year 1850, so there cant be a match, and the algorithm definitely searched the whole layer for these 7 features: 7*1000 iterations (minimum, dont know about the other features...))

The newly implemented option in V1.3 ("Do not search for matches on ...") can prevent a bad case (like in the last example setting): When there can't be a match, it can prevent from searching at all in that case. But this option can also be useful in other cases where one wants to limit the match searching.

---

Possible improvements:

So far, not beeing an expert in Python, I think the performance maybe could be improved by using an attribute index. This could allow to directly filter out (non) matches. However, my question if such a thing exists is still unanswered, and I am not sure if my "workaround" would increase the performance a lot. However, I am not sure if an attribute index improves performance here in combination with a spatial index. (But thinking about the whole thing again, the performance and result can be improved by adding another input, where you would manually add a year, where below no search at all would be done, e.g. the minimum year where there cant be a result anyways. --> just implemented that one in V1.3)

Surley there may be other things as well, a Python expert can optimize on this script. Processing tools run as background threads, thats why QGIS is still responsive while executing it.

---

Why not using expressions here?

Now lets come to the question on why expressions are a lot slower here: To my opinion/knowledge the reasons are: they are not using a spatial index and especially, the calculation is feature based. So basically overlay_nearest() is aggregating the whole layer as many times as the layer features has. The big difference to the Python script is, that we are only building an spatial index once and do not aggregate anything. We get the nearest neighbors by using this index on each iteration. As already said, I am not 100% sure here, so I welcome if my statement can be confirmed or disproved.