# Programmatically finding polygons which are >90% overlapped by another vector polygon layer using QGIS?

I'm trying to figure out how to use python to extract the polygons in one vector that are overlapped by >90% by another vector. I would then like to have a vector/map that will only show those polygons. The example picture shows my layers. I want all the grey polygons that are > 90% red.

I need to do this all via python (or similarly automated methods). I have ~1000 maps to process the same way.

• You want to do an overlay 'union' (see infogeoblog.wordpress.com/2013/01/08/geo-processing-in-qgis for some basics) then for each original polygon calculate the 'in' statistics and 'out' statistics gis.stackexchange.com/questions/43037/… to determine the overlay percent... hint: you need to have an area measurement gis.stackexchange.com/questions/23355/… – Michael Stimson May 30 '16 at 0:43
• Thanks for the tips. That's the same approach I was just attempting. I can do the union through the python console easy enough. Already added in the area attribute values. It's the next step I'm unsure of. How do I use python to calculate the 'in' and 'out' statistics so that I can identify/select/clip etc. the >90% polygons? – dnormous May 30 '16 at 0:50
• I think it's possible without python. Do you need absolutly python or a solution with virtual layers is good for you ? – Pierma May 30 '16 at 1:10
• The 'in' areas will have attributes from both polygons, the 'out' areas only have attributes from one set of polygons. Get both sets of area statistics and join back to the original polygons, add a field for 'in', 'out' and coverage, calculate the values for 'in' and 'out' from the sum of areas then divide 'in' by the original area (or 'in' + 'out') to calculate percent. – Michael Stimson May 30 '16 at 2:14
• Pierma - I just need an automated method to find the polygons. – dnormous May 30 '16 at 2:25

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])

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.

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).

• I'm getting an error when I run the query via the 'create a virtual layer' button. "Query execution error on CREATE TEMP VIEW _tview AS WITH r AS ("....rest of code here... followed by: "1 - near "WITH": syntax error" I'm fairly new to QGIS. Can I create this virtual layer programatically? Thanks for your help! – dnormous May 30 '16 at 3:48
• Sorry, a bad copy between grey and gray (sorry for my aproximative english). I edited the query. It should works now. Why do you want to create the layer progamatically ? The advantages of virtual layer is that it is non-destructive, and if you edit your data (the gray or the red polygons), the virtual layer will autoatically update. – Pierma May 30 '16 at 4:17
• This is only one small piece of the process. I have ~1000 of these maps to do, so automating the process will be extremely helpful. – dnormous May 30 '16 at 4:28
• I'm still getting the same error --> "1 - near "WITH": syntax error". I plugged the local names for each layer in for the grayLayer and redLayer. Is the local name what I should be using? ie: grey layer is labeled as "Basins_800", so I have code like "Basins_800.geometry" – dnormous May 30 '16 at 4:38

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])