# How to determine 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 ? – radouxju Oct 17 '16 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. – Ian Turton Oct 17 '16 at 13:26

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 – Ian Turton Oct 18 '16 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. – Andy Harfoot Oct 19 '16 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. – Andy Harfoot Oct 20 '16 at 15:12

Actually, you already did most of the work required to determine the tiles that you want to print using 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 grid as the coverage layer in 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 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 be sure that the scale will not be changed if anything happen, as you can see below; Using an index map, you can easily see the ID and location of each tile with reference to other tile, 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 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 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 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. – Victor Oct 19 '16 at 15:13
• @Victor Thanks for your comment, I updated my answer. – ahmadhanb Oct 20 '16 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. – Victor Oct 20 '16 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. – ahmadhanb Oct 20 '16 at 9:02
• just noting the fact, in case you have another good idea! Especially since my grid is not regular! – Victor Oct 20 '16 at 9:50

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.