[GIS] Determining neighbouring tile ids in QGIS

atlasprint-composerpyqgis

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:

enter image description here

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.

enter image description here

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.

Best Answer

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:

enter image description here

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:

enter image description here

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:

enter image description here

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;

enter image description here

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.

enter image description here

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:

enter image description here

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:

enter image description here

  • 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:

enter image description here

enter image description here

enter image description here

enter image description here

Related Question