Determining neighbouring tile ids in QGIS

Question

I was asked in a recent training course if QGIS could automatically calculate the next/previous and above/below page numbers for a map book created using the atlas generator. I managed to work out a fairly reasonable label expression for a regular grid if you know the width and height of the grid.

But we then started to think of realistic examples where we don't want to draw pages that don't contain our district of interest, such as this one of my home county:

So this afternoon I had a play at a Python script to work out the 4 neighbours I was interested in for each grid cell and added those values to my grid (this is heavily based on Ujaval Gandhi's tutorial):

for f in feature_dict.values():
    print 'Working on %s' % f[_NAME_FIELD]
    geom = f.geometry()
    # Find all features that intersect the bounding box of the current feature.
    # We use spatial index to find the features intersecting the bounding box
    # of the current feature. This will narrow down the features that we need
    # to check neighboring features.
    intersecting_ids = index.intersects(geom.boundingBox())
    # Initalize neighbors list and sum
    neighbors = []
    neighbors_sum = 0
    for intersecting_id in intersecting_ids:
        # Look up the feature from the dictionary
        intersecting_f = feature_dict[intersecting_id]
        int_geom = intersecting_f.geometry()
        centroid = geom.centroid()
        height = geom.boundingBox().height()
        width = geom.boundingBox().width()
        # For our purpose we consider a feature as 'neighbor' if it touches or
        # intersects a feature. We use the 'disjoint' predicate to satisfy
        # these conditions. So if a feature is not disjoint, it is a neighbor.
        if (f != intersecting_f and
            not int_geom.disjoint(geom)):
            above_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x(),
               centroid.asPoint().y()+height))
            below_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x(),
               centroid.asPoint().y()-height))
            left_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x()-width,
               centroid.asPoint().y()))
            right_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x()+width,
               centroid.asPoint().y()))
            above = int_geom.contains(above_point)   
            below = int_geom.contains(below_point)   
            left = int_geom.contains(left_point)
            right = int_geom.contains(right_point)
            if above:
                print "setting %d as above %d"%(intersecting_f['id'],f['id'])
                f['above']=intersecting_f['id']
                
            if below:
                print "setting %d as below %d"%(intersecting_f['id'],f['id'])
                f['below']=intersecting_f['id']
    
            if left:
                print "setting %d as left of %d"%(intersecting_f['id'],f['id'])
                f['left']=intersecting_f['id']
                
            if right:
                print "setting %d as right of %d"%(intersecting_f['id'],f['id'])
                f['right']=intersecting_f['id']
    
    # Update the layer with new attribute values.
    layer.updateFeature(f)

layer.commitChanges()

This works just fine.

But to be honest the whole creating a test point to the North and then testing all the possible neighbours seems wrong. However after an afternoon of wracking my brain I can't think of a better way to determine what is a particular grid cell's northern neighbour?

Ideally I'd like something simple enough to put in a print composer text box, but I suspect that's too much to ask for.

Answer 1

You already did most of the work required to determine the tiles that you want to print using an atlas. But the point is how to adjust everything together to show only the tile IDs that you need. To demonstrate my idea, I will use in this example a DEM image and a grid vector file, as you can see below:

First, we need to show the label of each grid.

In the layout view, I used the grid as the coverage layer in the atlas, I created two maps: the main view window map, and an index map that shows only the grid, as you can see below:

Then I did the following:

1. I adjusted the scale of the index map to show the entire grid extent then I fixed the scale
2. I fixed the view's extent to prevent the map from panning when using Preview atlas, and
3. I enabled the Overview to see the extent and location of the main view map, as you can see below:

For the main view window map, I fixed the scale to the extent of each grid block, to make sure that the scale will not be changed if anything happens, as you can see below;

Using an index map, you can easily see the ID and location of each tile concerning other tiles, even when you turn off the grid from the main view map window. For example, the following map has a tile ID=14, and you can see the surrounding tile IDs.

Update:

I will update my answer because I realized that you wanted to show the page number index of the surrounding layouts not the IDs of surrounding layouts.

To ease the understanding of the process, I will update the ID numbers in the index map to show the layout page number, as shown below:

Since the IDs that I have start from 0 (Zero), the ID of the first grid shown on the index map will start from 3. Therefore, I want to change the page number to start from 1 by subtracting 2 from the ID number in Atlas: Page number: ID -2, then I will use the current page number as a reference in the expression to create labels for the current page, previous page, next page, up page and below page, as follows:

• Current Page has this expression in the label text box: Current Page Number: [%@atlas_pagename%]

• Previous Page expression: [%if((@atlas_pagename = 1), Null, '↑ Page Number: ' || (@atlas_pagename - 1))%] since there are no pages before 1

• Next Page expression: [%if( (@atlas_pagename = 25), Null, '↓ Page Number: ' || (@atlas_pagename + 1))%] since there are no pages after 25

• Up Page expression: [%if((@atlas_pagename <= 6),NULL,'↑ Page Number: ' || (@atlas_pagename -6))%] since there are no pages before 6 in the upper direction

• Below Page expression: [%if((@atlas_pagename >= 20), Null, '↓ Page Number: ' || (@atlas_pagename + 6))%] since there are no pages after 20 in the lower direction

Some output results:

Answer 2

If you are not fitting each page extent (from the index layer) exactly into the composer, but instead having overlapping borders with adjacent pages (as shown in your second screenshot), then you could use labels from the index layer, with the downside that they would be inside the map border.

If there isn't any overlap, then you could replicate a technique that I used successfully in the past (coincidently across E & W Sussex!) in MapInfo, where I wrote a small script that generated a set of four points for each index feature, offset into the adjacent features, with attributes of both the sheet number, and the direction of offset. The point layer was then used to generate labels again, with the direction of offset allowing the orientation of the labels to be adjusted for a nicer effect.

I haven't tried this, but you may be able to avoid generating a separate data layer in QGIS through the use of the new geometry generator styling functionality, this would make for a more elegant and dynamic solution that wasn't achievable in MapInfo!

Answer 3

This solution will work for rectangular grids and it's automatic (should work for any scenario without adjusting anything manually).

Let's assume you have a grid with page numbers. You can run my Processing script selecting the grid layer and its page number field as parameters. The script creates four fields (right, left, above, below) in the grid layer and calculates the corresponding neighbor page id for each grid cell. Then you can use your expressions (e.g., [% if( "left" is not NULL, 'to page' || "left", "" ) %]) to show neighbor page labels.

Just add my repository (https://github.com/gacarrillor/QGIS-Resources.git) from QGIS Resource Sharing plugin and install the script:

How it works

The script determines the relation (right, left, above, or below) by comparing bounding boxes coordinates from both the current grid cell and each intersecting cell. It turns out that for each relation, one of the coordinates is missing.

If the relation is above, the missing coordinate is yMin, i.e., all other 3 coordinates from current grid cell's bounding box will be present in the above cell's bounding box. Remember that QGIS bounding boxes are defined in this order: [xMin, yMin, xMax, yMax].

For a numeric example let's take rectangles with sides of length 1. Say current cell's bounding box is defined as bbox1=[0,0,1,1]. The above cell's bounding box would be defined as bbox2=[0,1,1,2]. X coordinates from bbox1 are present in bbox2, whereas bbox1's yMin is missing in Y coordinates of bbox2.

We can define our 4 relations in this way (o: present, #: missing):

right: [#,o,o,o]
above: [o,#,o,o]
left:  [o,o,#,o]
below: [o,o,o,#]

As you can see, the missing index gives us all information we need.