Mit QGIS programmgesteuert Polygone finden, die zu> 90% von einer anderen Vektorpolygonschicht überlappt werden?

Ich versuche herauszufinden, wie man mit Python die Polygone in einem Vektor extrahiert, die von einem anderen Vektor um> 90% überlappt werden. Ich hätte dann gerne einen Vektor / eine Karte, die nur diese Polygone zeigt. Das Beispielbild zeigt meine Ebenen. Ich möchte alle grauen Polygone, die> 90% rot sind.

Ich muss dies alles über Python (oder ähnlich automatisierte Methoden) tun. Ich habe ~ 1000 Karten, um auf die gleiche Weise zu verarbeiten.

Der nächste Code funktioniert in meiner Python-Konsole von QGIS. Es wird eine Speicherschicht mit Polygonen erzeugt, die zu> 90% von roten Bereichen überlappt sind.

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)

Ich habe den Code mit diesen beiden Vektorebenen ausprobiert:

Nachdem der Code in der Python-Konsole von QGIS ausgeführt wurde, wurden zur Bestätigung der Ergebnisse die Indizes i, j der beteiligten Features, Schnittbereiche, Feldattribute in polygons_intersects (1 für rote Bereiche und 2 für graue Bereiche) und das Überlappungskriterium gedruckt .

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

Die erstellte Speicherschicht (grüne Merkmale) kann beim nächsten Bild beobachtet werden. Es war wie erwartet.

Hier eine Lösung, die kein Python benötigt.

Fügen Sie eine neue virtuelle Ebene mit einer Abfrage wie der folgenden hinzu:

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

Mit:

• Basins800 als Ihre Ebene möchten Sie mit grauen Polygonen filtern

• Schweregrad: Ihre rote Schicht überlappt sich.

Das Ergebnis ist eine neue Schicht, in der nur alle grauen Plolygone> 90% von roten Polygonen überlappt werden. Ein neues Feld enthält den Überlappungsprozentsatz.

Hoffe das funktioniert. Ich kann bei Bedarf weitere Details zur Abfrage hinzufügen.

Hinweis: Ihre Daten enthalten sehr kleine Polygone (die aus Ihrer Rasterverarbeitung stammen und einem Rasterpixel entsprechen (auf dem Bild sehen Sie 4 Polygone, aber es gibt 25 andere kleine Polygone). Dadurch wird die Ausführung der Abfrage sehr langsam (Schnittfunktion) generiert ein Feature für jedes Feature-Paar aus den beiden Ebenen).

Nachdem ich den Link zu den Shapefiles Severity und Basins800 gesehen hatte , konnte ich den erforderlichen Geoprozess verstehen. Ich habe den Code geändert in:

Mit QGIS programmgesteuert Polygone finden, die zu> 90% von einer anderen Vektorpolygonschicht überlappt werden?

für dieses bekommen:

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)

Nachdem ich den Code mit diesen Shapefiles in der Python-Konsole von QGIS ausgeführt hatte, erhielt ich in wenigen Minuten ein ähnliches Ergebnis wie Pierma . wo die Speicherschicht 31 Merkmale hatte (verschiedene von 29 Polygonen, die er bekommen hatte).

Ich werde die Ergebnisse nicht debuggen, da es 1901 * 3528 = 6706728 Interaktionen für Features gibt. Der Code sieht jedoch vielversprechend aus.