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

enter image description here

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.

enter image description here

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

  • related: gis.stackexchange.com/questions/60045/…
    – Mapperz
    Jan 15 '16 at 18:19
  • Are you open to python scripting? I'm not sure you can achieve what you want with a simple combination of QGIS functions.
    – ArMoraer
    Jan 15 '16 at 18:37
  • I'm just starting to write a simple script, but I need more info. The most important being : how do you know when multiple polygons are part of the same Natura 2000 area? Does the Natura 2000 layer have an attribute such as "area id" shared by polygons of a same area, for example?
    – ArMoraer
    Jan 15 '16 at 19:41

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:

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

        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): Test layer with natura 2000 areas in blue

  • Thanks so much for this! I tried running the script as suggested but got the following error: Traceback (most recent call last): File "C:/Users/sion elis.williams/.qgis2/python/plugins\scriptrunner\scriptrunner.py", line 553, in run_script user_script.run_script(self.iface) AttributeError: 'module' object has no attribute 'run_script' Jan 18 '16 at 13:58
  • Well, apparently you can't name your script script.py, for some reason... Can you try naming it myscript.py?
    – ArMoraer
    Jan 18 '16 at 14:29
  • Hello, thanks again. Renamed but this time got this: Traceback (most recent call last): File "C:/Users/sion elis.williams/.qgis2/python/plugins\scriptrunner\scriptrunner.py", line 553, in run_script user_script.run_script(self.iface) File "C:/Users/sion elis.williams/Desktop\myscript.py", line 24, in run_script if constit_geom.contains( n2000_geom ): Exception: unknown Jan 18 '16 at 18:31
  • Hi, sorry for being late. It's difficult to identify the issue from the message you got. It's possible that one of your polygons has a "wrong" geometry. I updated the script in order to ignore this kind of message, can you give it a try?
    – ArMoraer
    Jan 20 '16 at 15:17
  • Hey, no worries! Tried again but it was "Crash dumped" to C:\Users\SIONEL~1.WIL\AppData\Local\Temp\qgis-20160125-113527-5632-992-e2a51df.dmp which I'm struggling to find a way of opening. Will reply again once I've found a way. Thanks for your patience and persistence. Jan 25 '16 at 12:55

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.