Trouver par programme des polygones qui sont> 90% chevauchés par une autre couche de polygones vectoriels en utilisant QGIS?

J'essaie de comprendre comment utiliser python pour extraire les polygones d'un vecteur qui se chevauchent de> 90% par un autre vecteur. J'aimerais alors avoir un vecteur / carte qui n'affichera que ces polygones. L'image d'exemple montre mes calques. Je veux tous les polygones gris qui sont> 90% rouges.

J'ai besoin de tout faire via python (ou des méthodes automatisées similaires). J'ai ~ 1000 cartes à traiter de la même manière.

Le code suivant fonctionne dans ma console Python de QGIS. Il produit une couche mémoire avec des polygones qui sont> 90% chevauchés par des zones rouges.

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

#for polygon_intersects
feats_lyr1 = [ feat for feat in layers[0].getFeatures() ]

#for xwRcl
feats_lyr2 = [ feat for feat in layers[1].getFeatures() ]

selected_feats = []

for i, feat1 in enumerate(feats_lyr1):
    area1 = 0
    area2 = 0
    for j, feat2 in enumerate(feats_lyr2):
        if feat1.geometry().intersects(feat2.geometry()):
            area = feat1.geometry().intersection(feat2.geometry()).area()
            print i, j, area, feat2.attribute('class')
            if feat2.attribute('class') == 1:
                area1 += area
            else:
                area2 += area
    crit = area1/(area1 + area2)
    print crit
    if crit > 0.9:
        selected_feats.append(feat1)

epsg = layers[0].crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           "mem_layer",
                           "memory")

prov = mem_layer.dataProvider()

for i, feat in enumerate(selected_feats):
    feat.setAttributes([i])

prov.addFeatures(selected_feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

J'ai essayé le code avec ces deux couches vectorielles:

Après avoir exécuté le code sur la console Python de QGIS, pour corroborer les résultats, les index i, j des entités impliquées, les zones d'intersection, l'attribut de champ dans polygons_intersects (1 pour les zones rouges et 2 pour les zones grises) et le critère de chevauchement ont été imprimés. .

0 0 9454207.56892 1
0 1 17429206.7906 2
0 2 10326705.2376 2
0 4 40775341.6814 1
0 5 26342803.0964 2
0 7 11875753.3216 2
0.432253120382
1 6 1198411.02558 2
1 7 1545489.96614 2
1 10 27511427.9909 1
0.90930850584
2 7 750262.940888 2
2 8 12012343.5859 1
0.941213972294
3 6 23321277.5158 2
0.0

La couche de mémoire créée (traits verts) peut être observée à l'image suivante. C'était comme prévu.

Voici une solution qui ne nécessite pas de python.

Ajoutez une nouvelle couche virtuelle avec une requête comme:

WITH r AS (
SELECT 
    Basins800.rowid AS idGray, 
    area(Basins800.geometry) AS areaGray, 
    area(Intersection(Basins800.geometry, Severity.geometry)) AS aeraInter, 
    Basins800.geometry AS geomGray 
  FROM Basins800, Severity
)

SELECT *, areaInterSum/areaGray  AS overlap , geomGray 
    FROM (
        SELECT 
           idGray, 
           areaGray, 
           sum(areaInter) AS areaInterSum, 
           geomGray 
        FROM r 
        GROUP BY idGray) 
     WHERE areaInterSum/areaGray > 0.9

Avec :

• Basins800 comme couche que vous souhaitez filtrer avec des polygones gris

• Gravité: votre couche rouge se chevauche.

Le résultat sera une nouvelle couche avec seulement tous les plolygones gris> 90% chevauchés par des polygones rouges, avec un nouveau champ contenant le pourcentage de chevauchement.

J'espère que cela fonctionne. Je peux ajouter plus de détails sur la requête si nécessaire.

Remarque: Vos données contiennent de très petits polygones (provenant de votre traitement raster et correspondant à un pixel raster (sur l'image, nous pouvons voir 4 polygones mais il y a 25 autres petits polygones). Cela rend la requête très lente à exécuter (fonction Intersection génère une entité pour chaque couple d'entités à partir des deux couches).

Après avoir vu le lien vers les fichiers de formes Severity et Basins800 , j'ai pu comprendre le géoprocessus nécessaire. J'ai modifié le code dans:

Trouver par programme des polygones qui sont> 90% chevauchés par une autre couche de polygones vectoriels en utilisant QGIS?

pour obtenir celui-ci:

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

#for Severity
feats_lyr1 = [ feat for feat in layers[0].getFeatures() ]

#for Basins800
feats_lyr2 = [ feat for feat in layers[1].getFeatures() ]

selected_feats = []

print "processing..."

for i, feat1 in enumerate(feats_lyr1):
    for j, feat2 in enumerate(feats_lyr2):
        if feat1.geometry().intersects(feat2.geometry()):
            area1 = feat1.geometry().intersection(feat2.geometry()).area()
            area2 = feat1.geometry().area()
            print i, j, area1, area2
    crit = area1/area2
    print crit
    if crit > 0.9:
        selected_feats.append(feat1)

epsg = layers[0].crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           "mem_layer",
                           "memory")

prov = mem_layer.dataProvider()

for i, feat in enumerate(selected_feats):
    feat.setAttributes([i])

prov.addFeatures(selected_feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

Après avoir exécuté le code, avec ces fichiers de formes sur la console Python de QGIS, en quelques minutes, j'ai obtenu un résultat similaire à Pierma ; où la couche mémoire avait 31 entités (différentes des 29 polygones obtenus par lui).

Je ne vais pas déboguer les résultats car il existe des interactions 1901 * 3528 = 6706728 pour les fonctionnalités. Cependant, le code semble prometteur.