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.

  • 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
  • 1
    Pierma - I just need an automated method to find the polygons. – dnormous May 30 '16 at 2:25
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

5

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

  • 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
  • Here's a link to download the shapefiles: link – dnormous May 30 '16 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. – 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
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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.