# Determining neighbouring tile ids in QGIS

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.

• what if there is no neighbor on one side. Do you want the value of the closeset cell in one direction or would you leave a void ? Oct 17, 2016 at 13:24
• I'm happy for a null in that case, I can easily set the label to only display when not null or empty. Oct 17, 2016 at 13:26
• How is it possible that your Grid is named from left to right instead of from top to down? I would like to do this too. even better, first line A, second line B,... and numbering from left to right. Is this possible? Apr 14, 2021 at 8:02
• @Carina - see gis.stackexchange.com/questions/330760/… Apr 14, 2021 at 8:35

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:

• Although useful, this does not answer his question. Oct 19, 2016 at 15:13
• @Victor Thanks for your comment, I updated my answer. Oct 20, 2016 at 3:18
• this works in your example (and his), since the sides of the keymap/grid are regular. If they were not straight it would not work since the number to add or subtract (6 in your example) would vary depending on the atlas page you are on. Oct 20, 2016 at 8:56
• I agree with you. If the grid is not regular the process will be more complicated. But since he wants to apply it on a regular grid the method applied in my suggested solution will work. Oct 20, 2016 at 9:02
• just noting the fact, in case you have another good idea! Especially since my grid is not regular! Oct 20, 2016 at 9:50

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!

• I really should have thought of just using the labels of the other polygons! :-) After a quick experiment with the Geometry Generator I can draw a bounding box, but it's harder to build a grid Oct 18, 2016 at 16:56
• I was thinking along the lines of generating label points offset into the adjacent polygons, rather than grids. Another option would be to expand the index feature's MBR into the adjacent features to allow labels to be drawn. Oct 19, 2016 at 13:47
• Just had a play, and it appears that the geometry generated by the geometry generator styling doesn't get labelled, so isn't the more elegant solution that I had hoped for. Oct 20, 2016 at 15:12

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.

• Germán Carrillo, your link to the script does not work. I would need it. Could you post the link? Thank you very much in advance! Apr 14, 2021 at 8:46
• You're right. The script is no longer there. I deleted it because it was written for QGIS v2.x So, please confirm that you really really :) need it, so I can make some time (soon) for migrating it to QGIS 3.x. Apr 14, 2021 at 13:56
• Germán Carrillo, thank you for your offer! It is not about life or death ; ) I will keep on trying. But if you happen to write a script for QGIS 3 in the next months, think of me ; ) Apr 15, 2021 at 6:43
• Great! If you don't mind, open an issue here, so that I can keep the script under the radar. Apr 15, 2021 at 19:57