Sunday, 16 June 2019

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


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.



Answer



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!


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...