Counting polygons in one layer that are overlapping/touching/contained in each polygon in separate layer of contiguous polygons using QGIS?

Question

I have a layer of contiguous (pink) polygons (UK Parliamentary constitutencies) and a separate layer containing (blue) polygons representing protected nature areas (Natura 2000 areas):

Each constituency consists of a single polygon. Each Natura 2000 area can consist of one or more polygons.

For each constituency I need a combined total for the number of Natura 2000 areas that are:

1. contained in
2. touched by
3. overlapping

that constituency.

I've added both as vector layers in QGIS and tried using Spatial Query and Topology Checker plugins and Geoprocessing tools (Clip, Intersect). But I don't really know what I'm doing. Although I was able to select constituencies where there is an intersect I couldn't see the number of intersects per constituency and this didn't solve my need to include touching and contained polygons in the totals.

And just to be clear about Natura 2000 areas consisting of one or more polygons. In this screenshot (contained in a box) is an example of a single Natura 2000 site consisting of multiple polygons. Because all of these polygons are contained in a single constituency we'd only want to count them once as a single Natura 2000 area. But if they were spread across the boundaries of constituencies we'd count them for each constituency any part of them touched/overlapped/was contained by. In the same way we'd count the single polygon Natura 2000 site that overlaps two constituencies (circled) twice - once for each constituency it covers.

(For background, I'm campaigning to protect Natura 2000 areas and my aim is to develop a target list of UK Members of Parliament prioritised by those who have the highest number of these areas in their constituency.

Answer 1

EDIT: there is actually a simpler answer. Compute the intersection between your two layers (Vector > Geoprocessing tools > Intersect), then use the Group Stats plugin to compute statistics on the constituencies. You can even compute the total area of Natura 2000 sites per constituency.

Minor issue: it will only take into account Natura 2000 sites that intersect (overlap + contained) the constituencies.

Original answer:

As stated in my last comment I wrote a small Python script to do the job. It is probably the most simple way to answer your question.

Following the comment of mgc, I assume your Natura 2000 layer is made of MultiPolygons. If it is not the case (i.e. if the layer is made of simple polygons), just use Vector > Geometry Tools > Singleparts to Multipart... to get a new MultiPolygon layer from your Natura 2000 layer. You will have to specify the field that identifies unique areas (cf. my last comment).

The simplest way to execute the script is to install the Script Runner plugin in QGIS. Then, copy the following text in a file that you will call e.g. script.py:

from qgis.core import *

def run_script(iface):

    lyr_n2000 = QgsVectorLayer("PATH/TO/YOUR/NATURA2000/SHAPEFILE", "", "ogr")
    lyr_constit = QgsVectorLayer("PATH/TO/YOUR/CONSTITUENCIES/SHAPEFILE", "", "ogr")

    list_n2000 = list(lyr_n2000.getFeatures())
    list_constit = list(lyr_constit.getFeatures())

    print "Constituency | contained | touched | overlapped | TOTAL"

    for constit in list_constit:

        n2000_contained = 0
        n2000_touched = 0
        n2000_overlapping = 0
        constit_geom = constit.geometry()

        for n2000 in list_n2000:

            try:
                n2000_geom = n2000.geometry()

                if constit_geom.contains( n2000_geom ):
                    n2000_contained = n2000_contained + 1

                if constit_geom.touches( n2000_geom ):
                    n2000_touched = n2000_touched + 1

                if constit_geom.intersects( n2000_geom ):
                    n2000_overlapping = n2000_overlapping + 1

            except Exception:
                print "Error with n2000 #{0} and constit #{1}".format(n2000.id(), constit.id())
                pass

        print "{0} | {1} | {2} | {3} | {4}".format(constit.attribute("id"), n2000_contained, n2000_touched, n2000_overlapping-n2000_contained-n2000_touched, n2000_overlapping)

Change the paths to your shapefiles, and replace the "id" in constit.attribute("id") (last line) by the name of the field of the constituency layer you want to see appear.

Open the script with the Script Runner plugin (launch the plugin and click on the green "+"). To run the script, click on the green arrow. The results will appear in the bottom panel, they look like that (I used homemade test layers) :

Constituency | contained | touched | overlapped | TOTAL
1 | 0 | 0 | 0 | 0
2 | 0 | 0 | 0 | 0
3 | 0 | 0 | 1 | 1
4 | 0 | 0 | 0 | 0
5 | 0 | 0 | 1 | 1
6 | 1 | 1 | 1 | 3
7 | 0 | 0 | 1 | 1

The first column represents the constituency. The second, third and fourth columns are respectively the number of N2000 areas completely contained, adjacent or simply overlapping (not intersecting!) the constituency. As you can guess, the last column is the sum of the previous three.

This script may need some adjustments depending on the precise nature of your data. Don't hesitate to post a comment if anything goes wrong (bad results, script not running, etc).

Test layer (N2000 in blue, constit. in orange):