Czy programowo znajdujesz wielokąty, które w ponad 90% pokrywają się z inną warstwą wielokątów wektorowych za pomocą QGIS?

Próbuję wymyślić, jak użyć Pythona do wyodrębnienia wielokątów w jednym wektorze, które nakładają się o> 90% przez inny wektor. Chciałbym wtedy mieć wektor / mapę, która pokaże tylko te wielokąty. Przykładowe zdjęcie pokazuje moje warstwy. Chcę wszystkie szare wielokąty, które są> 90% czerwone.

Muszę to wszystko zrobić za pomocą Pythona (lub podobnie zautomatyzowanych metod). Mam ~ 1000 map do przetworzenia w ten sam sposób.

Następny kod działa w mojej konsoli Python QGIS. Tworzy warstwę pamięci z wielokątami, które w ponad 90% pokrywają się z czerwonymi obszarami.

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)

Wypróbowałem kod z tymi dwiema warstwami wektorowymi:

Po uruchomieniu kodu w konsoli Python w QGIS, w celu potwierdzenia wyników, wydrukowano indeksy i, j zaangażowanych cech, obszary przecięcia, atrybut pola w polygons_intersects (1 dla obszarów czerwonych i 2 dla obszarów szarych) oraz kryterium nakładania się .

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

Utworzoną warstwę pamięci (elementy zielone) można zaobserwować na następnym obrazie. Tak było zgodnie z oczekiwaniami.

Tutaj rozwiązanie, które nie wymaga Pythona.

Dodaj nową warstwę wirtualną za pomocą zapytania takiego jak:

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

Z :

• Basins800 jako warstwę, którą chcesz filtrować za pomocą szarych wielokątów

• Istotność: nakładanie się czerwonej warstwy.

Rezultatem będzie nowa warstwa z tylko wszystkimi szarymi ploligonami> 90% nakładającymi się na czerwone wielokąty, z nowym polem zawierającym procent nakładania się.

Mam nadzieję, że to zadziała. W razie potrzeby mogę dodać więcej szczegółów do zapytania.

Uwaga: Twoje dane zawierają bardzo małe wielokąty (pochodzące z przetwarzania rastrowego i odpowiadające pikselowi rastrowemu (na zdjęciu widzimy 4 wielokąty, ale jest 25 innych małych wielokątów). To sprawia, że ​​zapytanie jest bardzo powolne do wykonania (funkcja Przecięcie generuje jeden element dla każdej pary elementów z dwóch warstw).

Po obejrzeniu linku do plików kształtu Severity i Basins800 mogłem zrozumieć niezbędny geoproces. Zmodyfikowałem kod w:

Czy programowo znajdujesz wielokąty, które w ponad 90% pokrywają się z inną warstwą wielokątów wektorowych za pomocą QGIS?

za uzyskanie tego:

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)

Po uruchomieniu kodu z tymi plikami kształtów w Konsoli Pythona QGIS, w ciągu kilku minut uzyskałem podobny wynik jak Pierma ; gdzie warstwa pamięci miała 31 cech (różnych z 29 uzyskanych przez niego wielokątów).

Nie zamierzam debugować wyników, ponieważ istnieją interakcje 1901 * 3528 = 6706728 dla funkcji. Jednak kod wygląda obiecująco.