I primarily use QGIS.

I have two Layers which both contain polygons.So one layer(A) contains only one polygon, the other contains a number of polygons(layer (B)).

I can run the Intersect-Tool and QGIS returns the polygons of layer (B) which completly or partly lie within the polygon of layer (A).

But is it possible that QGIS returns a classification of those intersecting polygons of layer (B) by "lying completly" inside the the polygon of layer (A) and "lying partly" inside the polygon of layer (A).

Maybe this could be done by calculating and comparing the area of the polygons of layer(B) before and after intersecting, but I would like to avoid that for reasons.

And how can I automate this with a script?

  • You can use 'Vector > Spatial Query' to select polygons that fall within other polygons. But I don't know how to do that in python – Knightshound Sep 16 '16 at 8:39

It can be done using spatial queries pretty easily, first you can select the polygons from B that lie completely within the polygon layer A by using the st_contains function. Classify these as being completely within the polygon A.

Next up you can select all the polygons from B that intersect A using the st_intersects function, but this will also select the polygons that are completely within A, so you will have to remove the results of the first query/classification.

As far as QGIS is concerned, you can use the Spatial Query plugin to do this, enter image description here

  • Isn't is so that the big polygon Contains and small ones are Within? – user30184 Sep 16 '16 at 9:21
  • Yup, that is how it is. – Hasan Mustafa Sep 16 '16 at 9:27

You can otherwise do it with the field calculator of layer B with the expression:

within( $geometry,geometry( get_feature('LayerA, 'id', 1))) 

Providing that your layer A is called 'LayerA' and that it has an 'id' field with a value of '1'. You can of course modify all those settings.

The structure of get_feature is the following: get_feature(layer, attribute, value).

The expression will return 1 for the polygons entirely inside and 0 for the others.

example of successful implementation


For automating with a script you can use next code:

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

feats_A = [ feat for feat in layers[0].getFeatures() ]

feats_B = layers[1].getFeatures().next()

for i, feat in enumerate(feats_A):
    if feat.geometry().intersects(feats_B.geometry()):
        if feats_B.geometry().contains(feat.geometry()):
            print "total overlapping, ", "feature: ", i
            areaA = feat.geometry().intersection(feats_B.geometry()).area()
            areaB = feats.geometry().area()
            print "partial overlapping, feature: %d, perc: %.2f" % (i, (areaA/areaB)*100)

I try it out with two shapefiles for a similar condition in your question:

enter image description here

After running the code at Python Console of QGIS, it can be observed at above image, that there were printed only the three possible cases and its respective ids and percentages (except total overlapping: 100 %)

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