La búsqueda de polígonos que están> 90% superpuestos por otro polígono vector capa utilizando QGIS?

Question

Estoy tratando de averiguar cómo utilizar python para extraer los polígonos en un vector, que se superponen por >90% por otro vector. Entonces me gustaría tener un vector/mapa del que sólo se muestran los polígonos. El ejemplo muestra la imagen de mis capas. Quiero que todo el gris de los polígonos que se encuentran > 90% de rojo.

Tengo que hacer todo esto a través de python (o similar métodos automatizados). He ~1000 mapas de proceso de la misma manera.

Answer 1

Aquí una solución sin python necesidades.

Añadir un nuevo virtual de la capa con una consulta como :

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

Con :

• Basins800 es la capa que desea filtrar con gris polígonos

• Gravedad: el rojo de la capa de superposición.

El resultado será una nueva capa con sólo todo el gris plolygons >90% con la que se superponen rojo polígonos, con un nuevo campo que contiene la superposición por ciento.

Espero que esto funcione. Puedo agregar más detalles sobre la consulta si es necesario.

Nota : Sus datos contiene muy pequeños polígonos (que viene de su trama de procesamiento y correspondiente a una trama de píxeles (en la imagen, podemos ver 4 polígonos, pero hay otros 25 pequeños polígonos). Esto hace que la consulta muy largo para ejecutar (Intersección de la función de generar una característica para cada par de características de las dos capas).

Answer 2

Siguiente código funciona en mi Consola de Python de QGIS. Produce una capa de memoria con los polígonos que son > 90% con la que se superponen las áreas rojas.

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)

He probado el código con estas dos capas vectoriales:

Después de ejecutar el código en la Consola de Python de QGIS, para la corroboración de los resultados, no fueron impresos los índices i, j de características, la intersección de las áreas, atributo de campo en polygons_intersects (1 áreas rojas y 2 para las áreas grises) y la superposición de criterio.

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

Al crear una capa de memoria (verde) se pueden observar en la siguiente imagen. Se fue como se esperaba.

Answer 3

Después de ver el enlace a la Gravedad y Basins800 shapefiles, yo podía entender el necesario geoprocess. He modificado el código en:

Mediante programación encontrar polígonos que son >90% superpuesta por otro vector capa de polígonos utilizando QGIS?

para llegar a este:

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)

Después de ejecutar el código, con estos archivos en la Consola de Python de QGIS, en un par de minutos me dieron un resultado similar al de Pierma; donde la capa de memoria tenía 31 (de las características de los diferentes, de 29 de polígonos tengo por él).

No voy a depurar los resultados, porque hay 1901*3528 = 6706728 de las interacciones de las características. Sin embargo, el código se ve prometedor.