referencia a la geometría de las características en una expresión

Question

Tengo una capaA con carreteras y una capaB con un camino más corto.

Si pongo una función agregada en un campo virtual en la capaB para calcular la suma del campoAA donde el camino más corto contiene características de la capaA, devuelve la suma de cada característica en la capaA.

He probado diferentes expresiones y todas devuelven el mismo resultado (suma de todas las características):

---

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

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

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

---

En la documentación encontrará

Así que me preguntaba a qué geometría se refieren ambas y cómo puedo realizar una expresión como contains(geometry currentfeature layerA, geometry features layerB)

También son bienvenidas las sugerencias sobre cómo obtener el mismo resultado con un pythoncode personalizado en el editor de funciones.

Answer 1

Hasta donde yo sé, la expresión del filtro (cuarto parámetro de aggregate ) sólo puede acceder a los campos y valores de la capa referenciada ( layerA en tu ejemplo) y no de la capa en la que vas a almacenar los resultados (es decir, layerB ). Dos posibles soluciones podrían ser:

1. Para utilizar el Join attributes by location de Procesamiento (seleccione contains como Geometric predicate y el valor sum para Statistics for summary ).
2. Para utilizar el Function Editor y utilizar código Python.

Me explayaré más sobre la segunda solución.

---

Una función con los siguientes parámetros sería útil en su caso:

• Nombre de la capa con todas las carreteras ( layerA para ti).
• Nombre del campo de layerA para agregar ( fieldAA para ti).
• Geometría actual del elemento (de la capa con el camino más corto).

Por lo tanto, vaya a Field Calculator abra el Function Editor crea un nuevo archivo y pega este código:

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()) )

Haga clic en Cargar para registrar las funciones anteriores, vaya a la sección Expression y escribe esta expresión:

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

Todavía en el Field Calculator , cree/elija su campo de destino (de layerB ) y haga clic en OK .

Debería obtener la suma de fieldAA sólo de las características de layerA que contiene cada característica en layerB .

---

_NOTA 1: El sum_contains() utiliza una caché para evitar la construcción de todo el índice espacial para cada característica (fuente: Cómo utilizar el editor de funciones en la calculadora de campos de QGIS )._

_NOTA 2: Una vez que utilice la expresión antes mencionada la caché se construye para la capa de carreteras, por lo que, si edita características de dicha capa y desea utilizar la función sum_contains() de nuevo, tendrá que ir a Function Editor seleccione el archivo creado y haga clic en Load para que el índice se inicialice de nuevo. De lo contrario, estarías utilizando un índice obsoleto._

---

EDIT: Añadir comentarios en línea para explicar el código.

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()) )