Generate raster tiles from QGIS project, avoid creating redundant tiles (not overlapping a vector layer), containing only OSM basemap pixels

Question

Context: raster tiles webmap on private server

Using QGIS, I created raster tiles from a project with Generate XYZ Tiles (Directory). I basically have a vector layer with buildings and an OSM basemap in the background. I uploaded the generated tiles to a private webserver to create a webmap out from my QGIS project. This works fine.

Problem: a lot of redundant tiles, containing only OSM basemap

However, for the extent and zoom levels I use, over 2 million tiles were generated. Many tiles (for larger zoom levels probably the vast majority) do not contain any features (buildings) from my vector layer, but just the OSM background map (forests, uninhabited areas etc). See screenshots at the bottom.

So my idea was that it would be better just to create my own tiles for the areas where I have vector features (cities, villages etc). I could then serve these tiles from my own server and for the other areas (missing in my tiles), my webmap would have a regular OSM basemap in the background. I tested it, in principle this works: OSM basemap covering the whole area in the background, above it my tiles for inhabited areas.

As I use QGIS specific stylings (like blending modes etc.), vector tiles are not a soulution in this case.

Question: avoid redundant tiles to save processing time, number of files and storage

So now the question arises: how can I generate only those tiles that overlap with my vector layer - and not for the whole extent? In principle I could think of two possible approaches: avoid creating "empty" (containing only OSM basemap) tiles at all (ideal solution) so saving a lot of workload when updating the map. Or: creating all tiles and then find a way to delete those that are not needed.

Screenshot that illustrates the problem: I only have buildings (vector features) at the bottom; the upper half of QGIS map canvas contains only OSM basemap - still, for this part, tiles are generated as well:

Output directory, containing the generated tiles. In the upper part tiles containing only the OSM basemap (those I want to avoid), in the lower part those I need with vector features:

Answer 1

Have a go taking this code and modifying it to suit

This takes a full sized grid extent, and produces a subset of vector grid tiles with a level, row and column id for the full grid where the tiles intersect an area of relevance

I'm going to use this vector in FME for cutting out tiles from a larger raster layer for Spectrum Spatial, but if you were to do this in QGIS perhaps you could pair this with some sort of looping that runs through every tile and uses Clip raster by extent or Clip raster by mask layer

You could use a definition query (Extract by expression) for the looping perhaps, as well as to only render out a specific level

import time
startTime = time.time()
#Variable setting
minimumXFullExtent = -40867.905
maximumXFullExtent = 3329926.909
maximumYFullExtent = 8483438.288
minimumYFullExtent = 5112643.474
widthAndHeightOfSquareExtent = 3370794.814
theSpacing = widthAndHeightOfSquareExtent
coordSys = 'EPSG:7855'
extent = str(minimumXFullExtent) + ',' + str(maximumXFullExtent) + ',' + str(minimumYFullExtent) + ',' + str(maximumYFullExtent) + ' [' + coordSys + ']'
directory = 'D:/BasemapTiling/Bounds/'
name = 'BasemapTilingV1'
smallerAreaForTiling = 'D:/Whatever/YourRelevantArea.gpkg'
listOfLevels = []
#Start at the highest tiling level and work your way down
for x in range(1,17):
    print(x)
#Halve the tile dimensions each time
theSpacing = theSpacing/2

#Make a grid based on the original extent or on the extent of the latest subset
processing.run("native:creategrid", {'TYPE':2,'EXTENT':extent,
'HSPACING':str(theSpacing),'VSPACING':str(theSpacing),'HOVERLAY':0,'VOVERLAY':0,'CRS':QgsCoordinateReferenceSystem(coordSys),'OUTPUT':directory + name + 'Level' + str(x) + '.gpkg'})

#Remove fields to save space
processing.run("qgis:deletecolumn", {'INPUT':directory + name + 'Level' + str(x) + '.gpkg','COLUMN':['id','left','top','right','bottom'],
'OUTPUT':directory + name + 'Level' + str(x) + 'Reduce.gpkg'})

#Extract only the tiles relevant to the subset
processing.run("native:extractbylocation", {'INPUT':directory + name + 'Level' + str(x) + 'Reduce.gpkg','PREDICATE':[0,1,4,5,6],
'INTERSECT':smallerAreaForTiling,'OUTPUT':directory + name + 'Level' + str(x) + 'ReduceExtract.gpkg'})

#Calculate the column number 
processing.run("native:fieldcalculator", {'INPUT':directory + name + 'Level' + str(x) + 'ReduceExtract.gpkg','FIELD_NAME':'Column','FIELD_TYPE':1,'FIELD_LENGTH':0,'FIELD_PRECISION':0,
'FORMULA':'round (\r\n ( x ( centroid ( $geometry ) ) - (' + str(minimumXFullExtent) + '+ 1))\r\n / (\r\n range ( x ( centroid ( $geometry ) ) ) \r\n / (\r\n count_distinct ( x ( centroid ( $geometry ) ) )\r\n -1\r\n )\r\n )\r\n)\r\n',
'OUTPUT':directory + name + 'Level' + str(x) + 'ReduceExtractWColumn.gpkg'})

#Calculate the row number 
processing.run("native:fieldcalculator", {'INPUT':directory + name + 'Level' + str(x) + 'ReduceExtractWColumn.gpkg','FIELD_NAME':'Row','FIELD_TYPE':1,'FIELD_LENGTH':0,'FIELD_PRECISION':0,
'FORMULA':'round (\r\n ( ' + str(maximumYFullExtent) + '- 1 \r\n - \r\n y ( centroid ( $geometry ) ) \r\n  ) / \r\n ( range ( y ( centroid ( $geometry ) ) ) \r\n  / \r\n ( count_distinct ( \r\n ( y ( centroid ( $geometry ) ) ) , x ( centroid ( $geometry ) ) \r\n ) -1 )\r\n )\r\n)\r\n',
'OUTPUT':directory + name + 'Level' + str(x) + 'ReduceExtractWColumnWRow.gpkg'})

#Add the level
processing.run("native:fieldcalculator", {'INPUT':directory + name + 'Level' + str(x) + 'ReduceExtractWColumnWRow.gpkg','FIELD_NAME':'Level','FIELD_TYPE':1,'FIELD_LENGTH':0,'FIELD_PRECISION':0,
'FORMULA':str(x),'OUTPUT':directory + name + 'Level' + str(x) + 'ReduceExtractWColumnWRowWLevel.gpkg'})

#Add to the list of what has been made so far
listOfLevels.append(directory + name + 'Level' + str(x) + 'ReduceExtractWColumnWRowWLevel.gpkg')

#Get the extent of the latest grid extract to make the next grid
extentVector = QgsVectorLayer(directory + name + 'Level' + str(x) + 'ReduceExtract.gpkg')
extentRectangle = extentVector.extent()
xminExtentRectangle = extentRectangle.xMinimum()
xmaxExtentRectangle = extentRectangle.xMaximum()
yminExtentRectangle = extentRectangle.yMinimum()
ymaxExtentRectangle = extentRectangle.yMaximum()
coordsExtentRectangle = "%f,%f,%f,%f" %(xminExtentRectangle, xmaxExtentRectangle, yminExtentRectangle, ymaxExtentRectangle)
extent = coordsExtentRectangle + ' [' + coordSys + ']'

#Prevent this from running for ages if you're just testing
currentTime = time.time()
timeSoFar = currentTime - startTime
if timeSoFar > 20000:
    break


#Bring it all together
processing.run("native:mergevectorlayers", {'LAYERS':listOfLevels,
'CRS':QgsCoordinateReferenceSystem(coordSys),'OUTPUT':directory + 'MergedLevels' + name + '.gpkg'})

For this I used some modified code from Vector grid id generation

This is a visualisation of the levels of the grid and their attributes