I have a folder with a gdal2tiles output: 6 levels of zoom (and folders), with subfolders and .png tiles.
Is there a gdal tool to merge in a single geotiff file these .png tiles?
My solution requires a fair amount of work to get there, but this should work.
Set up TileCache, GeoWebCache, or MapProxy to consume the tile cache and output a WMS service.
You could then consume this WMS service with QGis or another desktop GIS that supports WMS.
Create your printable map using the desktop GIS.
gdal_merge. It should be included with most distributions of GDAL.
You'll have to generate a world file for each image first. You would use just the images from the largest scale directory (6) each subdirectory corresponds to a column of the image and each the highest number image in the subdirectory is the top row. There should be an XML file from the gdal2tiles output that has the extents. You'll have to calculate the upper left pixel coordinate of each image and pixel size (probably 1) from that.
So you might have something like:
1.000000000000 0 0 -1.000000000000 130000 2180000
where the lines correspond to these values from the XML file:
units-per-pixel for the lowest level tileset 0 (no rotation) 0 (no rotation) negative units-per-pixel for the lowest level tileset Origin x BoundingBox maxY
This will take some scripting. The world files should end in .pgw, .pngw or .wld.