QGIS – Referencing Features Geometry in an Expression for Advanced Analysis

expressionfunctionpyqgisqgisqgis-custom-function

I have a layerA with roads and a layerB with one shortest path.

If I put an aggregate function into a virtual field in layerB to calculate the sum of fieldAA where the shortest path contains features of layerA, it returns the sum of every feature in layerA.

I tried different expression and they returned all the same result (sum of every feature):

aggregate( 'layerA', 'sum', "fieldAA", contains($geometry, $geometry))

aggregate( 'layerA', 'sum', "fieldAA", contains(geometry($currentfeature), $geometry))

aggregate( 'layerA', 'sum', "fieldAA", within($geometry, geometry($currentfeature)))

In the documentation you find

enter image description here

So I was wondering to which geometry they both refer to and how I can perform an expression like contains(geometry currentfeature layerA, geometry features layerB)

Suggestions on how to get the same result with a custom pythoncode in the function editor are also welcome.

Best Answer

As far as I know, the filter expression (fourth parameter of aggregate) can only access fields and values from the referenced layer (layerA in your example) and not from the layer you're going to store results into (i.e., layerB). Two possible solutions could be:

  1. To use the Join attributes by location from Processing (select contains as Geometric predicate and the value sum for Statistics for summary).
  2. To use the Function Editor and use Python code.

I'll elaborate more on the second solution.

A function with the following parameters would be handy in your case:

  • Name of the layer with all roads (layerA for you).
  • Name of the field from layerA to aggregate (fieldAA for you).
  • Current feature geometry (from layer with the shortest path).

So, go to Field Calculator, open the Function Editor, create a new file, and paste this code:

from qgis.core import ( qgsfunction, QgsSpatialIndex, QgsFeatureRequest, 
    QgsMapLayerRegistry )

index_cache = {}

def get_layer_index( layer ):
    if layer.name() not in index_cache:
        index_cache[layer.name()] = QgsSpatialIndex( layer.getFeatures() )
    return index_cache[layer.name()]

@qgsfunction(args='auto', group='Custom')
def sum_contains(road_layer_name, field_name, geometry, feature, parent):
    layers = QgsMapLayerRegistry.instance().mapLayersByName( road_layer_name )
    if not layers:
        raise Exception( "The layer could not be found." )
    layer = layers[0]
    intersectIds = get_layer_index( layer ).intersects( geometry.boundingBox() )
    features = layer.getFeatures( QgsFeatureRequest().setFilterFids( intersectIds ) )
    return sum( f[field_name] for f in features if geometry.contains(f.geometry()) )

Click on Load to register the functions above, go to the Expression tab and write this expression:

sum_contains( 'layerA', 'fieldAA',  $geometry )

Still in the Field Calculator, create/choose your target field (from layerB) and click OK.

You should get the sum of fieldAA from only those features in layerA contained by each feature in layerB.

NOTE 1: The sum_contains() function uses a cache to avoid building the whole spatial index for each feature (source: How to Use Function Editor in QGIS Field Calculator).

NOTE 2: Once you use the aforementioned expression the cache is built for the road layer, so, if you edit features from such layer and you want to use the sum_contains() function again, you would need to go to Function Editor, select the file you created, and click on Load so that the index is initialized again. Otherwise, you would be using an outdated index.

EDIT: Adding in-line comments to explain the code.

from qgis.core import ( qgsfunction, QgsSpatialIndex, QgsFeatureRequest, 
    QgsMapLayerRegistry )

index_cache = {} # Global variable to store spatial index for layers

def get_layer_index( layer ):
    # Python function to check whether we already built an index for
    # a given layer. If so, just return the built index, otherwise
    # build it now and store it in the global variable for further calls.
    if layer.name() not in index_cache:
        index_cache[layer.name()] = QgsSpatialIndex( layer.getFeatures() )
    return index_cache[layer.name()]

@qgsfunction(args='auto', group='Custom')
def sum_contains(road_layer_name, field_name, geometry, feature, parent):

    # First, check if we can access the road layer 
    layers = QgsMapLayerRegistry.instance().mapLayersByName( road_layer_name )
    if not layers:
        raise Exception( "The layer could not be found." )
    layer = layers[0]

    # Second, obtain ids of road features that intersect the
    # current geometry BBox. This uses a spatial index 
    # to optimize the spatial search. Features that intersect the BBox
    # likely intersect the geometry itself, but we need to check that later.
    intersectIds = get_layer_index( layer ).intersects( geometry.boundingBox() )

    # From the step above we get a list of matching feature ids,
    # Now we use them to request the feature objects.
    features = layer.getFeatures( QgsFeatureRequest().setFilterFids( intersectIds ) )

    # Now we perform a spatial search using the geometry predicate
    # 'contains' with current geometry and all features that 
    # intersect its BBox, because they are great candidates to match.
    # We do this in one line (see Python list comprehensions): 
    # iterate (road) features, validate the 'contains' predicate, 
    # and list field values of matching features. Finally, we sum 
    # the list of values and return the total.
    return sum( f[field_name] for f in features if geometry.contains(f.geometry()) )
Related Question