Найти объекты других слоев, которые покрывают объект (обнаружение наложений)

Question

На следующем рисунке вы видите 2 слоя. Слой 1 разделен на более мелкие части, чем слой 2, но все объекты в слое 2 состоят из одного или нескольких объектов из слоя 1.

Что я хотел бы сделать, так это с помощью python повторить функции слоя 2 и обнаружить все функции уровня 1, которые эта функция "покрывает" (пересечение полигонов?). Таким образом, я смог бы определить числа, напечатанные красным цветом для каждой особенности слоя 2.

Есть ли простой фрагмент кода для выполнения этой задачи?

Я пробовал использовать feature.geometry().intersection(base_feature.geometry()) и проверка на != None, но, похоже, это не помогает. Он проверяет только точки, насколько я могу судить, и мне нужно проверить, перекрывается ли область / пересекается.

Answer 1

Краткий ответ: вместоintersection()использоватьintersects().

---

Вот высокоэффективный код для PyQGIS 3 с использованием QgsSpatialIndex() , а также QgsGeometryEngine . Объяснения в виде комментариев в коде.

Обратите внимание, что в зависимости от желаемого результата вы можете использовать другой геометрический предикат, чемintersects, такой какoverlaps, в строке 25 (if layer2_feat_geometryengine.intersects(layer1_feat_geom.constGet()):). Подробное описание различий смотрите здесь.

      ### Settings:
layer1 = QgsProject.instance().mapLayersByName('intersection_polygon1')[0]
layer2 = QgsProject.instance().mapLayersByName('intersection_polygon2')[0]

### Run:
# Create the spatial index and store the geometries in it, so we dont need a slow .getFeature() request
layer1_idx = QgsSpatialIndex(layer1.getFeatures(), flags=QgsSpatialIndex.FlagStoreFeatureGeometries)
# Iterate over all features of layer2
for layer2_feat in layer2.getFeatures():
    # Get the features geometry
    layer2_feat_geom = layer2_feat.geometry()
    # Create a geometryengine, which is much faster for mass-intersection tests than normal geometry, see https://api.qgis.org/api/classQgsGeometryEngine.html
    layer2_feat_geometryengine = QgsGeometry.createGeometryEngine(layer2_feat_geom.constGet())
    layer2_feat_geometryengine.prepareGeometry()
    # Get all features of layer1 that intersect with the current layer2 features bounding box
    # It does return a list of the feature id of the bbox-intersecting features, nothing more
    overlapping_features = layer1_idx.intersects(layer2_feat_geom.boundingBox())
    overlapping_ids = []
    # Iterate over the list of bbox-intersecting features
    for layer1_feat_id in overlapping_features:
        # Get the geometry of the bbox-intersecting feature from the spatial index, which is much faster than a feature request from the layer
        layer1_feat_geom = layer1_idx.geometry(layer1_feat_id)
        # Check if the two geometries really intersect, not just their bounding boxes
        # You can also use another predicate than intersects depending on your desired result, such as overlaps.
        if layer2_feat_geometryengine.intersects(layer1_feat_geom.constGet()):
            overlapping_ids.append(str(layer1_feat_id))
    print('Overlapping features for ' + str(layer2_feat.id()) + ':')
    print(','.join(overlapping_ids))