Efficient way to compare fields or expressions of intersecting features: bulk-loading vs. get features one-by-one

Question

I want to find an efficient (especially in terms of memoryusage and processing speed) method to compare intersecting features by attribute/expression.

For this script I am iterating over a source_layer_vl and for each source-feature I am iterating over the intersecting overlay-features from an overlay_layer_vl which I take from a QgsSpatialIndexKDBush(). To compare the expression results I then need to get the feature or expression result somehow. This is where I am unsure of the advantages and disadvantages of the methods I have tried:

Bulk loading overlay-features into a dictionary before iterating over the source layer, and getting the intersecting features from this dictionary instead by feature request
Bulk loading overlay-expression-results into a dictionary before iterating over the source layer, and getting the expression-results of the intersecting features from this dictionary instead by feature request
Loading the matching overlay-features via getFeature(id) one-by-one as they intersect and evaluate the expression

Lets say about 50% of the points do actually intersect with the polygons.

The first method seems to be by far the fastest, almost 50% faster than method 3 and 25% faster than method 2. But: Isnt this very memory expensive, especially on larger layers? And why is bulk-loading so much faster than getting features one-by-one?

I am aware this question sounds opinion-based, but I am pretty sure it can be answered with pure facts.

Here is a reproducible code, containing all three methods, just change layer and fieldnames:

from datetime import datetime
import operator

# Vars
source_layer_vl = QgsProject.instance().mapLayersByName('earthquakepolygons')[0] # a polygonlayer
overlay_layer_vl = QgsProject.instance().mapLayersByName('earthquake2')[0] # a 2D-Single-Point-Layer
source_compare_expression = QgsExpression('"fid"') # just some expression or field of the polygonlayer
overlay_compare_expression = QgsExpression('"gap"') # just some expression or field of the pointlayer
op = operator.lt # the comparisonoperator (operator.lt means where source < overlay)

# All methods
overlay_layer_idx = QgsSpatialIndexKDBush(overlay_layer_vl.getFeatures())
if overlay_layer_idx.size() == 0:
    print('Overlay layer is empty or not of type 2D-Single-Point!')

### Method 1 ###########
method_01_starttime = datetime.now()
method_01_result = {}

overlay_layer_dict = {feat.id():feat for feat in overlay_layer_vl.getFeatures()}

for source_feat in source_layer_vl.getFeatures():
    method_01_result[source_feat.id()] = 0
    
    source_feat_geom = source_feat.geometry()
    source_feat_geometryengine = QgsGeometry.createGeometryEngine(source_feat_geom.constGet())
    source_feat_geometryengine.prepareGeometry()
    
    overlay_features = overlay_layer_idx.intersects(source_feat_geom.boundingBox())
    
    source_compare_expression_context = QgsExpressionContext()
    source_compare_expression_context.setFeature(source_feat)
    source_compare_expression_context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(source_layer_vl))
    source_compare_expression_result = source_compare_expression.evaluate(source_compare_expression_context)
    
    for overlay_feat in overlay_features:
        overlay_feat_geom = QgsGeometry.fromPointXY(overlay_feat.point()).constGet()
        if source_feat_geometryengine.intersects(overlay_feat_geom):
            overlay_real_feat = overlay_layer_dict[overlay_feat.id]
            overlay_compare_expression_context = QgsExpressionContext()
            overlay_compare_expression_context.setFeature(overlay_real_feat)
            overlay_compare_expression_context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(overlay_layer_vl))
            overlay_compare_expression_result = overlay_compare_expression.evaluate(overlay_compare_expression_context)
            if op(source_compare_expression_result, overlay_compare_expression_result):
                method_01_result[source_feat.id()] += 1
            
method_01_endtime = datetime.now()
method_01_runtime = method_01_endtime - method_01_starttime


### Method 2 ###########
method_02_starttime = datetime.now()
method_02_result = {}

overlay_layer_dict = {}
for overlay_feat in overlay_layer_vl.getFeatures():
    overlay_compare_expression_context = QgsExpressionContext()
    overlay_compare_expression_context.setFeature(overlay_feat)
    overlay_compare_expression_context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(overlay_layer_vl))
    overlay_compare_expression_result = overlay_compare_expression.evaluate(overlay_compare_expression_context)
    overlay_layer_dict[overlay_feat.id()] = overlay_compare_expression_result

for source_feat in source_layer_vl.getFeatures():
    method_02_result[source_feat.id()] = 0
    
    source_feat_geom = source_feat.geometry()
    source_feat_geometryengine = QgsGeometry.createGeometryEngine(source_feat_geom.constGet())
    source_feat_geometryengine.prepareGeometry()
    
    overlay_features = overlay_layer_idx.intersects(source_feat_geom.boundingBox())
    
    source_compare_expression_context = QgsExpressionContext()
    source_compare_expression_context.setFeature(source_feat)
    source_compare_expression_context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(source_layer_vl))
    source_compare_expression_result = source_compare_expression.evaluate(source_compare_expression_context)
    
    for overlay_feat in overlay_features:
        overlay_feat_geom = QgsGeometry.fromPointXY(overlay_feat.point()).constGet()
        if source_feat_geometryengine.intersects(overlay_feat_geom):
            overlay_compare_expression_result = overlay_layer_dict[overlay_feat.id]
            if op(source_compare_expression_result, overlay_compare_expression_result):
                method_02_result[source_feat.id()] += 1
            
method_02_endtime = datetime.now()
method_02_runtime = method_02_endtime - method_02_starttime


### Method 3 ###########
method_03_starttime = datetime.now()
method_03_result = {}

for source_feat in source_layer_vl.getFeatures():
    method_03_result[source_feat.id()] = 0
    
    source_feat_geom = source_feat.geometry()
    source_feat_geometryengine = QgsGeometry.createGeometryEngine(source_feat_geom.constGet())
    source_feat_geometryengine.prepareGeometry()
    
    overlay_features = overlay_layer_idx.intersects(source_feat_geom.boundingBox())
    
    source_compare_expression_context = QgsExpressionContext()
    source_compare_expression_context.setFeature(source_feat)
    source_compare_expression_context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(source_layer_vl))
    source_compare_expression_result = source_compare_expression.evaluate(source_compare_expression_context)
    
    for overlay_feat in overlay_features:
        overlay_feat_geom = QgsGeometry.fromPointXY(overlay_feat.point()).constGet()
        if source_feat_geometryengine.intersects(overlay_feat_geom):
            overlay_real_feat = overlay_layer_vl.getFeature(overlay_feat.id)
            overlay_compare_expression_context = QgsExpressionContext()
            overlay_compare_expression_context.setFeature(overlay_real_feat)
            overlay_compare_expression_context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(overlay_layer_vl))
            overlay_compare_expression_result = overlay_compare_expression.evaluate(overlay_compare_expression_context)
            if op(source_compare_expression_result, overlay_compare_expression_result):
                method_03_result[source_feat.id()] += 1
            
method_03_endtime = datetime.now()
method_03_runtime = method_03_endtime - method_03_starttime


### Results ###########
print('Method 01 Results:\n' + str(method_01_result))
print('Method 02 Results:\n' + str(method_02_result))
print('Method 03 Results:\n' + str(method_03_result))
print('\n')
print('Runtime Method 01 in Microseconds: ' + str(method_01_runtime.microseconds))
print('Runtime Method 02 in Microseconds: ' + str(method_02_runtime.microseconds))
print('Runtime Method 03 in Microseconds: ' + str(method_03_runtime.microseconds))

Kalak · Accepted Answer · 2023-01-13 13:44:54Z

Bulk loading will, in most cases, be faster than lazy loading unless your process can break/return early in the iteration (this is true for any iterative process not just python/qgis/getfeature/getfeatures) at the expense of some amount of memory usage.

If you look at the code of QgsVectorLayer.getFeature() here you can see that it actually calls QgsVectorLayer.getFeatures() filtered by the id.

Which "basically" iterates over all the features of the layer.

This means that your process can be approximated to doing ->

for poly in layer_poly.getFeatures():
    ids = idx_points.intersects(poly)
    for id in ids:
        for f in layer.getFeatures():
            if f.id() == id:
                # your comparaison code

As you can see that is fairly exponential. while with bulk loading it looks more like this:

for poly in layer_poly.getFeatures():
    ids = idx_points.intersects(poly)
    for id in ids:
        f = dictionary[id]
        # your comparaison code

Which is obviously quicker. you can't beat a dictionary lookup (having a complexity of O(1) Average (complexity)).

Also, the biggest performance impact in the code you have provided is actually these fourth lines

...
source_compare_expression_context = QgsExpressionContext()
...
source_compare_expression_context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(source_layer_vl))
...
overlay_compare_expression_context = QgsExpressionContext()
...
overlay_compare_expression_context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(overlay_layer_vl))

These doesn't need calculated at every iteration, they can be taken out of the loops for a significant performance improvement.

### Method 4 ###########
method_04_starttime = datetime.now()
method_04_result = {}

overlay_layer_dict = {feat.id(): feat for feat in overlay_layer_vl.getFeatures()}

source_compare_expression_context = QgsExpressionContext()
source_compare_expression_context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(source_layer_vl))

overlay_compare_expression_context = QgsExpressionContext()
overlay_compare_expression_context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(overlay_layer_vl))

for source_feat in source_layer_vl.getFeatures():
    method_04_result[source_feat.id()] = 0

    source_feat_geom = source_feat.geometry()
    source_feat_geometryengine = QgsGeometry.createGeometryEngine(source_feat_geom.constGet())
    source_feat_geometryengine.prepareGeometry()

    overlay_features = overlay_layer_idx.intersects(source_feat_geom.boundingBox())

    source_compare_expression_context.setFeature(source_feat)
    source_compare_expression_result = source_compare_expression.evaluate(source_compare_expression_context)

    for overlay_feat in overlay_features:
        overlay_feat_geom = QgsPoint(overlay_feat.point())
        if source_feat_geometryengine.intersects(overlay_feat_geom):
            overlay_real_feat = overlay_layer_dict[overlay_feat.id]

            overlay_compare_expression_context.setFeature(overlay_real_feat)
            overlay_compare_expression_result = overlay_compare_expression.evaluate(overlay_compare_expression_context)
            if op(source_compare_expression_result, overlay_compare_expression_result):
                method_04_result[source_feat.id()] += 4

method_04_endtime = datetime.now()
method_04_runtime = method_04_endtime - method_04_starttime

Runtime Method 01 in Microseconds: 550555

Runtime Method 02 in Microseconds: 635280

Runtime Method 03 in Microseconds: 749011

Runtime Method 04 in Microseconds: 43001

This method is over 10x faster than Method01 and 20x faster than method03.

Thanks for this nice and extensive answer. Just two questions on this: In your fourth method you use the said method of bulkloading all features into a dictionary. I was wondering about the memory usage of doing so. Pretty sure its really fast and efficient for small data, but what about big datasets? Can this somehow cause an overload? I mean thats basically a copy of the whole layer. And secondly I read that placing QgsExpressionContext() outside of the loop and assigning it only once for the whole layer can cause crashes. Never tried it though. Will get back to this and some testings soon. — MrXsquared
– MrXsquared, Commented Jan 13, 2023 at 14:12
I am not an expert on QgsExpressionScopes I rarely use QgsExpression.evaluate() but I don't see the problem with setting a layer scope outside the loop, if it was a feature scope it would probably be necessary however. Storing all features in a dictionary does take a fair bit of space you are right a 1 000 000 point geopackage with 1 int field weighing 157MB makes a dictionary weighing 75MB — Kalak
– Kalak, Commented Jan 13, 2023 at 14:22
Sometimes you have to decides a right balance between speed / performance, few improvements can help with both at once — Kalak
– Kalak, Commented Jan 13, 2023 at 14:26
About big datasets the dictionary will become faster and faster compared to getFeature() because while getFeature() is iterative the dictionary lookup is always a lookup of average complexity O(1) since the you don't actually look for the id in the dictionary keys but compute a hash and access the value corresponding to that hash in the dictionary — Kalak
– Kalak, Commented Jan 13, 2023 at 14:42

Stack Exchange Network

Efficient way to compare fields or expressions of intersecting features: bulk-loading vs. get features one-by-one

1 Answer 1

Hot Network Questions

Efficient way to compare fields or expressions of intersecting features: bulk-loading vs. get features one-by-one

1 Answer 1

Related

Hot Network Questions