9

Example layers

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.

Improve this question
6
  • 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/… Commented May 30, 2016 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? Commented May 30, 2016 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 ? Commented May 30, 2016 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. Commented May 30, 2016 at 2:14
  • 1
    Pierma - I just need an automated method to find the polygons. Commented May 30, 2016 at 2:25

3 Answers 3

Reset to default
6

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.

enter image description here

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

Improve this answer
9
  • 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! Commented May 30, 2016 at 3:48
  • Here's a link to download the shapefiles: link Commented May 30, 2016 at 4:01
  • 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. Commented May 30, 2016 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. Commented May 30, 2016 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" Commented May 30, 2016 at 4:38
3

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:

enter image description here

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.

enter image description here

Improve this answer
2
  • If I need to run the same in a single file, where my features are having >90% overlapping? Commented Jun 7, 2022 at 10:47
  • It is another question but you can use combinations method from itertools python module for same file (with only one features list). Commented Jun 7, 2022 at 12:57
2

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

enter image description here

I am not going to debug the results because there are 1901*3528 = 6706728 interactions for features. However, the code looks promising.

Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.