Programmatically finding polygons which are >90% overlapped by another vector polygon layer using QGIS?
Here a solution that does not require python.
Add a new virtual layer with a query like :
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
With :
Basins800 as your layer you want filter with grey polygons
Severity: your red layer overlapping.
The result will be a new layer with only all the grey plolygons >90% overlapped by red polygons, with a new field containing the overlap percent.
Hope this works. I can add more details on the query if needed.
Note : Your data contains very small polygons (coming from your raster processing and corresponding to a raster pixel (on the picture, we can see 4 polygons but there are 25 other small polygons). This make the query very slow to execute (Intersection function generates one feature for each couple of features from the two layers).
Next code works in my Python Console of QGIS. It produces a memory layer with polygons which are > 90% overlapped by red areas.
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)
I tried out the code with these two vector layers:
After running the code at the Python Console of QGIS, for corroborating results, there were printed the indexes i, j of involved features, intersection areas, attribute of field in polygons_intersects (1 for red areas and 2 for gray areas) and the overlapping criterion.
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
The created memory layer (green features) can be observed at the next image. It was as expected.
After seeing the link to Severity and Basins800 shapefiles, I could understand the necessary geoprocess. I modified the code in:
Programmatically finding polygons which are >90% overlapped by another vector polygon layer using QGIS?
for getting this one:
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)
After running the code, with these shapefiles at the Python Console of QGIS, in a few minutes I got a similar result as Pierma; where the memory layer had 31 features (different of 29 polygons got by him).
I am not going to debug the results because there are 1901*3528 = 6706728 interactions for features. However, the code looks promising.