Generate cover layer for atlas

Question

Let's say I have data for the area in the picture.

I created a model that generates a grid that I can use for my atlas. In my script, I enter:

Area of interest (couche d'emprise)
Scale of my print layout (echelle du plan)
Horizontal size of the map canevas in the layout (dimension horizontal de la carte (mm))
Vertical size of the map canevas in the layout (dimension verticale de la carte (mm))
Overlay I want between the atlas (distance de recouvrement entre les cartes (m))

As you can see, it works perfectly.

But I'm trying to optimize it. You can see that every rows and columns are aligned, but it is not useful. What I really want is to create as little rectangle as possible, to cover the entire area. For exemple, the column on the left can be shifted a bit down, so that I have only 3 rectancle, and not 4 anymore. The first row can be shifted to the right so that I would have only 2 grid at the end.

This optimization will be very useful for exemple if I need to do an atlas along a river for exemple. I don't want the grid to be aligned, only as few rectancle as possible to cover the entire river.

[EDIT] I used the expression of @Babel (with the adaptation), which works perfectly. But 2 things can be better. First, as you can see, if I use these vert_dist and hor_dist, the grid cells in the last row are smaller (vert_dist). How to modify the expression, so that the cells would always be the correct size?

Secondly, at which line should I enter a "distance overlapping" between the cells. In most cases, when we do an atlas, a bit of overlapping is desired.

Not sure to understand... You want to see the result I want? — katagena, May 02 '23 at 14:43
No, I mean: you're looking for a solution that shifts the grid cells automatically or are you also interested how to move the polygons manually? — Babel, May 02 '23 at 14:50
No I don't want to do it manually, that's why I created a script... But I'm trying to optimize it! By the way, I not searching to shift the grid cells. I used "grid cell" until now, cause I didn't find something else. What I'm looking for, is to cover my area with as few rectancles (with the dimensions I choosed) as possible. — katagena, May 02 '23 at 14:54
So what you want is to optimize the size (length,width) of the rectangles (grid cells)? Otherwise, if you have a pre-defined size, I don't see what you want to achieve without moving some of the cells...? Like moving C1 and D1 to the right so that no E1 is necessary any more? Than I mean be "shift grids" — Babel, May 02 '23 at 14:55
Not at all... I want to optimize their position, to minimazie the number of rectancles. What I'm saying is I used "grid cell" cause I was the only thing I found... If I need to use something else to achieve what I want, I will! — katagena, May 02 '23 at 15:00
How can you optimize if you don't want to change the size and not want to move them? Can you share a screenshot of how to result should look like? You want to create rectancles of a fixed size or not? You write "the column on the left can be shifted a bit down" - so you want to shift or not? It is quite confusing to understand what you intend to do. — Babel, May 02 '23 at 15:01
Of course I want to move them... It is what I wrote in the question with "shift"... You said you want to move your "grid cells"... I'm just saying there is perhaps a solution without create a grid cells and moving them, but generate "directlly" the correct rectancles. — katagena, May 02 '23 at 15:04

Babel · Accepted Answer · 2023-05-03T13:18:22.523

Algorithmic approach

Define a fixed width x and hight y for the rectangles (grid cells). Split your polygon in horizontal slices (bands, belts) at the interval of y (red outlined in screenshot).
Get the intersection of the bands of step 1 with the initial polygon (highlighted in yellow).
Get the bounding box of the intersection (black dotted line).
Now split the black dotted polygon vertically with a horizontal distance of x to get the polygon in dark blue line pattern.
Repeat these steps to the right until you reach the end of the bounding box of step 3, then go to the next line (band) and continue until you have split up all bands like this.

Implementation using QGIS expressions

Use this expression with Geometry Generator or Geometry by expression to create the grid cells for the coverage layer. Define the width and length of the cells in lines 2 and 3 (see last screenshot at the bottom):

collect_geometries (
with_variable ('ver_dist',22800,  -- change polygon width here
with_variable ('hor_dist',31000,  -- change polygon length here
    array_foreach (
        generate_series (0,(y_max($geometry)-ymin($geometry))/@ver_dist),
        with_variable (
            'bb',
            bounds (
            intersection (
                $geometry,
                bounds(
                    make_line (
                        make_point (x_min($geometry), y_max($geometry)-@ver_dist*@element),
                        make_point (x_max($geometry), y_max($geometry)-@ver_dist*@element),
                        make_point (x_max($geometry), y_max($geometry)-@ver_dist*(@element+1))
        )))),
        collect_geometries (
            array_foreach (
                generate_series (0,(x_max(@bb)-x_min(@bb))/@hor_dist),
                    bounds(
                        make_line (
                            make_point (x_min(@bb)+@hor_dist*@element,y_min(@bb)),
                            make_point (x_min(@bb)+@hor_dist*@element,y_max(@bb)),
                            project (make_point (x_min(@bb)+@hor_dist*@element,y_max(@bb)), @hor_dist, radians (90))
)))))))))

Changing the values, you can easily create differnt cell sizes. Using geometry genrator, you can see the changes in realtime - good for testing an ideal cellsize:

Variants

If you choose smaller sizes for the rectangles, depending on the shape of the polygon to cover, there will be empty cell. To create only cells that intersect the polygon, modify the expression: replace the last 6 lines (starting with bounds(... with this expression:

                 with_variable(
                     'bs',
                     bounds(
                         make_line (
                             make_point (x_min(@bb)+@hor_dist*@element,y_min(@bb)),
                             make_point (x_min(@bb)+@hor_dist*@element,y_max(@bb)),
                             project (make_point (x_min(@bb)+@hor_dist*@element,y_max(@bb)), @hor_dist, radians (90)))),
                         case 
                         when intersects (@bs, $geometry)
                         then @bs
                         else centroid ($geometry)
                         end
))))))))

Based on your comment: to avoid the grid cells in the last row to be smaller, instead of setting a fixed size for @ver_dist (vertical distance y), you can divide the whole north/south extent of the polygon in an even number so you have exactly matching hight of the grids to cover the whole polygon. You can than adopt the horizontal x distance to be the factor square root of 2 to have an aspect ratio corresponding to paper size A4, A3 etc. Modify lines 2 and 3 like this (multiply with sqrt(2) for landscape, divide by the same value for portrait orientation):
```
 with_variable ('ver_dist', (y_max($geometry)-y_min($geometry))/5,
 with_variable ('hor_dist',@ver_dist*sqrt(2), -- change multiply operator to divide operator for portrait (instead of landscape) orientation
```
Distance overlapping is normally done in the print composer when you set up the atlas, see next screenshot. If you want, you can do this already in the expression, to generate slightly overlapping polygons. Simply add a buffer at the last bounds function (6th last line in the above expression). I now set it to 10% of vertical distance, you can change the factor of 0.1 in the expression or set an absolute value for buffer size instead, as you like. See below for the modified expression below, including all variants discussed here.
I added an array_filter() function to avoid having an empty geometry, created to avoid empty polygons.

So the final expression, including the four variants, looks like this:

Screenshot, showing the overlapping buffers with 8 horizontal bands:

collect_geometries (
    with_variable ('ver_dist', (y_max($geometry)-y_min($geometry))/8,  -- change number of horizontal bands here
    with_variable ('hor_dist',@ver_dist*sqrt(2),  -- hange multiply operator to divide operator for portrait (instead of landscape) orientation
    array_foreach (
        generate_series (0,(y_max($geometry)-ymin($geometry))/@ver_dist),
        with_variable (
            'bb',
            bounds (
            intersection (
                $geometry,
                bounds(
                    make_line (
                        make_point (x_min($geometry), y_max($geometry)-@ver_dist*@element),
                        make_point (x_max($geometry), y_max($geometry)-@ver_dist*@element),
                        make_point (x_max($geometry), y_max($geometry)-@ver_dist*(@element+1))
        )))),
        collect_geometries (
            array_filter (
                array_foreach (
                    generate_series (0,(x_max(@bb)-x_min(@bb))/@hor_dist),
                    buffer (
                        with_variable(
                            'bs',
                            bounds(
                                make_line (
                                    make_point (x_min(@bb)+@hor_dist*@element,y_min(@bb)),
                                    make_point (x_min(@bb)+@hor_dist*@element,y_max(@bb)),
                                    project (make_point (x_min(@bb)+@hor_dist*@element,y_max(@bb)), @hor_dist, radians (90))
                            )),
                            case 
                            when intersects (@bs, $geometry)
                            then @bs
                            else NULL
                            end
                        ),
                        @ver_dist*0.1,  -- change size of overlap here, now set to 10% of vertical distance
                        join:='miter'
                )),
                @element is not NULL
)))))))

Variant with overlap without buffer

For your last comment, use this expression, where you can set fixed width/length for the rectangle and a value for overlap:

collect_geometries (
    with_variable ('ver_dist', 20000,  -- change number of horizontal bands here
    with_variable ('hor_dist',25000,  -- change multiply operator to divide operator for portrait (instead of landscape) orientation
    with_variable ('overlap',1000,  -- change overlap distance
    array_foreach (
        generate_series (1,ceil((y_max($geometry)-ymin($geometry))/@ver_dist)+1),
        with_variable (
            'bo',
             bounds(   
                intersection (
                    $geometry,
                    bounds(
                        extrude( 
                            make_line (
                                make_point (x_min($geometry), y_max($geometry)+@overlap*(@element)-@ver_dist*@element),
                                make_point (x_max($geometry), y_max($geometry)+@overlap*(@element)-@ver_dist*@element)
                            ), 
                            0, 
                            @ver_dist
            )))),
            with_variable(
                'bx',
                bounds(
                    make_line (
                        project (make_point (x_min(@bo)-@overlap,y_max(@bo)+@overlap*(@element=1)),@ver_dist, radians(180)),
                        make_point (x_min(@bo)-@overlap,y_max(@bo)+@overlap*(@element=1)),
                        make_point (x_max(@bo)+@overlap,y_max(@bo)+@overlap*(@element=1))
                )),
                collect_geometries (
                    array_filter (
                        array_foreach (
                            generate_series (0,((x_max(@bx)-x_min(@bx))/@hor_dist)+1),
                            with_variable(
                                'bs',
                                bounds(
                                    make_line (
                                        make_point (x_min(@bx)-@overlap*(@element+1)+@hor_dist*@element,y_min(@bx)),
                                        make_point (x_min(@bx)-@overlap*(@element+1)+@hor_dist*@element,y_max(@bx)),
                                        project (make_point (x_min(@bx)-@overlap*(@element+1)+@hor_dist*@element,y_max(@bx)), @hor_dist, radians (90))
                                )),
                                case 
                                when intersects (@bs, $geometry)
                                then @bs
                                else NULL
                                end
                        )),
                        @element is not NULL
)))))))))

You can repeat this process starting with vertical bands. For some exotic source shapes it can reduce count of pages two times. — RainForest, May 02 '23 at 16:16
Not quite sure what you mean: creating vertical bands as alternative to horizontal ones, depending on the shape of the polygon? Or vertical bands in addition to horizontal ones? — Babel, May 02 '23 at 17:16
Just repeat this algorytm, but start it with vertical bands. Then compare count of pages after first and second processes. In most cases it will be almost the same. The difference can occurs with some extraordinar shapes like long slanted lines for example. — RainForest, May 03 '23 at 04:20
@Babel: It is quiet perfect, always very impressive your codes. Just 2 improvments. See my edit above. — katagena, May 03 '23 at 06:20
With all these modifications, it is not anymore what I need. When we do a layout in QGIS, the map canevas is almost never A4, A0 or else. This is the size of the page. But the map canevas can be exotic. This is why I use the expression the size of the map canevas and the scale of the map caneva. So given a number of horizontal bands means nothing in this case. — katagena, May 03 '23 at 09:18
OK, no problem. Feel free to simply add absoulte values for the size, as before, for variables @ver_dist and @hor_dist in lines 2 and 3 to get what you want. — Babel, May 03 '23 at 09:21
Using buffer for the overlapping seems incorrect, cause the size of the rectangle must be constant. I think it is more in the "make_point" to ajust the x and y position. And using the maring aroud feature in "control by atlas" would change the scale. — katagena, May 03 '23 at 09:21
You add the same buffer to each grid cell, so rectangle size is constant (same with control by atlas option). — Babel, May 03 '23 at 09:25
Yes, but with this, I can't control the scale to be for example 1: 100'000. Let says I have a landscape A4 for my map canevas with a scale of 1: 100'000, the vertical distance is (210/10)*(100000/100), so 21000. If I add a buffer to this, at the end, my scale won't be 1: 100'000 anymore in my atlas. So the overlapping must be done in the position of the cells, without changing the size of them. — katagena, May 03 '23 at 09:30
See edited answer for overlap by moving the cells accordingly — Babel, May 03 '23 at 13:19

Generate cover layer for atlas

1 Answers1

Algorithmic approach

Implementation using QGIS expressions

Variants

Variant with overlap without buffer

Linked