5

Background:

At 3.20.3, there are two intersect tools available: Vector > Geoprocessing tools > Intersection and Processing Toolbox > SAGA > Features - Polygons > Intersect.

In both cases, only two layers can be intersected at a time. However, I have numerous layers that need to be intersected, and it would be much more efficient to conduct the task once, rather than two-at-a-time.

I notice that the Advanced ArcGIS Desktop 10.8.1 Intersection tool allows more than two input layers. I've searched the internet, plugins, and the documentation without success for a similar QGIS solution.

Question:

How do I intersect more than two layers at a time using QGIS?

3

Short answer: the principles

You can use QGIS expressions to generate custom overlay like intersection of as many layer as you want. Use Geometry generator or Geoemtry by expression (see here for differences).

How to do it

Let's say you have 5 polygon layers, name p1 to p5. Say you look for the areas where all five intersect. Split up this task in several steps. This is because the expression intersect() only takes two arguments (input layers), so you have to create a "nested" (combined) statement, say like (pseudocode):

intersect the (interesction of A and B) with the (interesction of C and D) 

Or more detailed:

  1. Intersection of p2 with p3 (lines 4-7 in the expression below)
  2. Intersection of p4 with p5 (lines 8-11)
  3. Intersect the result of steps 1 and 2 (you get: intersection of layers p2, p3, p4 and p5)
  4. Intersect p1 with the result of step 3.

The expression to use

The expression looks like this (see below at the bottom for detailed explanation of the functions used):

intersection (
    $geometry, 
    intersection (
        intersection (
            buffer (collect_geometries(overlay_intersects('p2', $geometry)),0), 
            buffer (collect_geometries (overlay_intersects( 'p3', $geometry)),0)
        ),
        intersection (
            buffer (collect_geometries(overlay_intersects('p4', $geometry)),0), 
            buffer (collect_geometries (overlay_intersects( 'p5', $geometry)),0)
        )
    )
)

Screenshot: the expression with Geometry generator: the diagonal X hached area is created with the expression above - the area where all 5 layers intersect:

enter image description here

Variant: how to adapt the expression to changing number of layers

For an even number of intersecting polygons, adapt the expression as follows (here: 6 polygons p1 to p6). As you see, I added a buffer([geometry],0) statement. Sometimes, without this, the expression returns NULL.

intersection (
    intersection (
        buffer (collect_geometries (overlay_intersects ('p5',$geometry)),0), 
        buffer (collect_geometries (overlay_intersects ('p6',$geometry)),0)
    ),
    intersection (
        intersection (
            buffer (collect_geometries (overlay_intersects ('p4',$geometry)),0), 
            buffer (collect_geometries (overlay_intersects ('p3',$geometry)),0)
        ),
        intersection(
            $geometry, 
            collect_geometries (overlay_intersects ('p2',$geometry))
        )
    )
)

Screenshot: the diagonal X hached area, representing the common space of all 6 polygons, is bordered by all 6 polgons as you can see from the color of the outer boundaries: red, blue, green, yellow, black and pink: enter image description here

Explanation: How this expression works

  • The intersection is done using the intersection() function, this only takes two arguments: intersection(geometry1, geometry2).
  • So we have to combine intermediate steps with nested (combined) intersection functions. Geometry generator as well as Geoemtry by expression are always based on a specific layer, so to get the geometries of the other layers to combine (intersect) them, we use the overlay_intersect() function.
  • This returns an array of all geometries, so to get geometries, we must convert from array to geometries using collect_geometries().
  • To work properly, we must buffer these geometries by 0, using buffer() function. For the very layer we work on, it's enough to get the geometries with $geometry.
8
  • This is very interesting! I just successfully tried it with five layers, using Geometry Generator. Now I want to use Geometry by Expression to create an output layer using six input layers. Since I'm not much of a programmer, how would your code be modified for an even number of input layers?
    – Stu Smith
    Oct 9 '21 at 21:59
  • The expression works the same also with Geoemtry by expression - just copy/paste, no need to program something.
    – Babel
    Oct 9 '21 at 22:01
  • Here's where I'm confused: The Geometry by Expression tool starts with one input layer. Then your code works with the remaining layers to be intersected in pairs. That means that there must be an odd number of total layers to be intersected (the required input layer, plus the pairs). However, in my case, there is an even number of layers to be intersected. How do I modify your code to accommodate an even number of layers to be intersected?
    – Stu Smith
    Oct 9 '21 at 22:08
  • Added a version using the intersection of 6 polygons: does that help? Also added a buffer() statement (buffering with 0) as without this, it doesn't always work properly.
    – Babel
    Oct 10 '21 at 10:02
  • Yes, this helps! Some comments: 1) I added a buffer() statement to the last line (the P2 layer), for consistency, and it works fine. 2) You can simply change the term intersection with union and Voila! 3) How to include overlay attribute fields? The only fields in the output layer are those from P1. I'll need to include fields from the overlay layers. Thoughts?
    – Stu Smith
    Oct 10 '21 at 19:06
1

As your initial question focused on intersection only, without attributes, my first answer did just that. As you specified later, you're interested in the attributes as well. So this is best done using a virtual layer, based on the same principles: as st_intersection() takes just two arguments, you have to combine several st_intersection() statements.

To create a new layer with the common area of four layers p1 to p4 plus all attributes, use this query in the dialog creating the "Virtual Layer" through Layer > Add Layer > Add/Edit Virtual Layer...:

select
    st_intersection(st_intersection(a.geometry, b.geometry), st_intersection(c.geometry, d.geometry)),
    *
from
    p1 as a, p2 as b, p3 as c, p4 as d

Screenshot: the hatched area in pink is the area where all four layers intersect. It's created as a virtual layer. As you can see from the attribute tables, it does not only create the geometry, but also adds all attributes, thanks to the , * statement at the end of the first line: enter image description here

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.