5

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.

  • 1
    As far as I know, the filter expression (forth 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). One solution could be to use the Function Editor and write Python code. – Germán Carrillo Sep 29 '16 at 14:17
  • how would such python code look like? – PieterB Sep 29 '16 at 14:20
  • It would look as the solution given by Detlev here: gis.stackexchange.com/questions/178522/…, although the Python code there can be optimized with a spatial index. You could also have a look at the "Join attributes by location" tool, as explained by wiltomap in the same thread. – Germán Carrillo Oct 1 '16 at 12:48
  • I know the 'join attributes by location', but I need contain, not intersect. I tried Detlev's approach, but I have multiple items within 1 original. So I need to sum them up. – PieterB Oct 1 '16 at 19:24
  • Have a look at the Processing tool 'Join Attributes by Location', there you can select 'Contains'. I'll elaborate more on the Python approach, which I find very interesting. – Germán Carrillo Oct 2 '16 at 1:32
2
+50

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()) )
  • Wonderfull. It does what I was looking for. I would like to use your answer as a base for learning pyqgis. Can you break up your code and give some information about what the different codeblocks doeing? – PieterB Oct 2 '16 at 19:43
  • I added a section with in-line comments to explain what the code does. – Germán Carrillo Oct 3 '16 at 17:59
  • OK, I can start studying... – PieterB Oct 3 '16 at 18:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.