Answering my own question, here's the solution I went with: roll-your-own in python.
The map in question uses scale-dependent rendering, so, I wanted to do a different render for each 'z' (integer TMS zoom level, 10 thru 17 in this case).
Overall workflow, once for each zoom level:
1) render from qgis to tif with a worldfile (i.e. the tif will have coordinates but no crs/transform info; the worldfile will have the standard worldfile info)
2) convert the rendered tif plus tfw into a real geotiff (i.e. crs/transform embedded in the tif file)
3) gdal2tiles to build the TMS tile set
some more detail:
- render from qgis
Since the various layers of the map have different CRSs, make sure to render via the composer instead of bypassing the composer. The composer capitalizes on the on-the-fly projection, so everything in the render looks like your map area in qgis when on-the-fly projection is enabled.
Since big render jobs will cause a qgis error, just do a bunch of smaller render jobs (hardcode some arbitrary 'render tile' size in paper units (millimeters or inches)), then plan on combining all the 'render tile' geotiffs for the current zoom level into a single VRT per zoom level in the next step.
- convert tif+tfw to GeoTIFF
for each 'render tile':
gdalwarp -s_srs <map_project_srs> <render_tile> <render_tile_out>
after all 'render tiles' are converted to real GeoTIFF:
gdalbuildvrt <output_vrt_name> <list/glob/listfile of render_tile_out files>
- gdal2tiles to build the TMS tile set
gdal2tiles -z <target_z_level> <output_vrt_name> <output_tile_dir>
which will put all the tms tiles into
<output_tile_dir>/<target_z_level> so you can specify the same output_tile_dir for each iteration of this step
I think I got it to the point where there's a limited set of configuration options that you have to get to the script by whatever method:
- x0,y0,x1,y1 - bounding box in map units of the entire area to build a TMS tile set of
- zmin, zmax - what TMS z zoom levels to build tiles for
- 'render tile' size in mm per side
- dpi - dots per inch, as used by the renderer
- a list of map scales that you want to correspond to the integer TMS 'z' zoom levels; this is for purposes of scale-dependent rendering in the source map layers
- base_dir - base directory where all the work will be saved
Would folks here recommend turning this into a plugin? It would be my first, but honestly it takes so friggin long and so friggin much disk space that it's been much better so far to do it in pieces. Probably the only difference between this and qtiles is the concept of chopping up the render task into small renders that won't error out, then building a vrt out of them - as well as the flexibility to determine what scale per z level you want, to account for scale dependent rendering.
Any ideas or suggestions for optimizing this flow?
UPDATE: sorry I forgot to include this credit from the script I did the first time around, referencing the link from zoneblue's solution:
# large portions of this code taken from
# and the qgis docs at