Skip to content

Printing Map Part II: Generate Map Tiles with QGIS

August 29, 2011

After map design has been complete, I needed to print tiles of it in various zoom levels. As far as I know, currently there are no tools to print map as tiles  in Qgis, so I have written my own script that prints the whole map in desired zoom level.

Please note, that current script is written for meters, not degrees, but it can easily be changed.

I had to print very specific tif files. It had to be:

  • Tiled
  • Tiles size 512 x 512
  • Colors in PCT
  • Compressed in LZW
  • It had to be GTiff (not tiff with world file)
  • Pixels count in height and width had to be divisible by 512
  • Meters per pixel had to be an integer

To do first 5 points, I’ve used GDAL, I will not post the code here, it is quite simple and you can even do such task with gdal_translate.

For 6th point it is not much to say, after having a script that prints map, you can easily modify it for such or similar things, and giving this information would just make source code more unreadable and complicated for reason that is not needed for anyone else but me (I suppose).

And for 7th point it is simply the initial data I provide for the script.

To print a TIF I have used info from Qgis tutorial.

I have written my own World File output function, because I couldn’t find native one. And it is very simple anyway:

def createWorldFile(fileName, mainScale):
	#>>Create World File
	mapRect = iface.mapCanvas().mapRenderer().extent()
	f = open(fileName + ".tifw", 'w')
	f.write(str(mainScale) + '\n')
	f.write(str(0) + '\n')
	f.write(str(0) + '\n')
	f.write('-' + str(mainScale) + '\n')
	f.write(str(mapRect.xMinimum() + mainScale/2) + '\n')
	f.write(str(mapRect.yMaximum() - mainScale/2))
	f.close()
	#<

For high zoom levels, I always ended up having blank files, so I don’t even print them, by using this code:

First take all the features you want to print. For me it was country polygon:

#take border stored in c:\\border.shp, so later I would not print blank files
#I have country polygon in border.shp
border = QgsVectorLayer('c:\\border.shp', 'border', 'ogr')
border.invertSelection()
borderFeatures = border.selectedFeatures()

Then when you have a rectangle you want to print, just check if it is inside your features:

for borderFeature in borderFeatures:
    if borderFeature.geometry().intersects(QgsRectangle(newRect.xMinimum(), newRect.yMinimum(), newRect.xMaximum(), newRect.yMaximum())):
        intersects = 1

First rectangle is set to be in left bottom corner and finishing rectangle is top right corner:

#start loop at bottom left corner.
iface.mapCanvas().zoomToFullExtent()
mapRenderer = iface.mapCanvas().mapRenderer()
mapRect = mapRenderer.extent()
xStart = mapRect.xMinimum()
xEnd = mapRect.xMaximum()

yStart = mapRect.yMinimum()
yEnd = mapRect.yMaximum()

Calculate meters to be seen in  rectangle:

width = mapRenderer.width()
height = mapRenderer.height()

#printed raster will have (width * metersPerPixel) meters in X and
#(height * metersPerPixel) meters in Y
xDiff = width * metersPerPixel
yDiff = height * metersPerPixel

And set extent for it:

#>>Generate first rectangle
newRect = QgsRectangle(xStart, yStart, xStart + xDiff, yStart + yDiff)
mapRenderer.setExtent(newRect)

After having this, there is just a loop required to finish this task:

#loop it!
while (newRect.yMinimum() < yEnd):
    while (newRect.xMinimum() < xEnd):

and regenerate rectangle:

    #move along x
    newRect = QgsRectangle(newRect.xMaximum(), newRect.yMinimum(), newRect.xMaximum() + xDiff, newRect.yMaximum())
    mapRenderer.setExtent(newRect)
#move along y
newRect = QgsRectangle(xStart, newRect.yMaximum(), xStart + xDiff, newRect.yMaximum() + yDiff)
mapRenderer.setExtent(newRect)

Inside loop there is a simple print block as explained in previously mentioned Qgis tutorial:

#set up composition
composition = QgsComposition(mapRenderer)
composition.setPrintResolution(256)
composition.setPaperSize(mapRenderer.width()/(mapRenderer.outputDpi()/25.4), mapRenderer.height()/(mapRenderer.outputDpi()/25.4))
composition.setPlotStyle(QgsComposition.Print)
dpi = composition.printResolution()
dpmm = dpi / 25.4 #get dots per mm
#add a map to the composition
composerMap = QgsComposerMap(composition,0,0,composition.paperWidth(),composition.paperHeight())
composition.addItem(composerMap)
#create output image and initialize it
image = QImage(QSize(width, height), QImage.Format_ARGB32) #output image size
image.setDotsPerMeterX(dpmm * 1000) #mm to meters
image.setDotsPerMeterY(dpmm * 1000) #mm to meters
image.fill(0)
# render the composition
imagePainter = QPainter(image)
sourceArea = QRectF(0, 0, composition.paperWidth(), composition.paperHeight() )
targetArea = QRectF(0, 0, width, height)
composition.render(imagePainter, targetArea, sourceArea)
imagePainter.end()
print "c:\\a" + str(i) + ".tif"
image.save("c:\\a" + str(i) + ".tif" , "tif")

I will post complete source code of the script.  It has comments inside it for those who still need them:

from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
import qgis.utils
from qgis.utils import iface
from qgis.gui import *

import sys
import os
import os.path
import shutil

def createWorldFile(fileName, mainScale):
    #>>Create World File
    mapRect = mapRenderer.extent()
    f = open(fileName + ".tifw", 'w')
    f.write(str(mainScale) + '\n')
    f.write(str(0) + '\n')
    f.write(str(0) + '\n')
    f.write('-' + str(mainScale) + '\n')
    f.write(str(mapRect.xMinimum() + mainScale/2) + '\n')
    f.write(str(mapRect.yMaximum() - mainScale/2))
    f.close()

metersPerPixel = 16
#take border stored in c:\\border.shp, so later I would not print blank files
#I have country polygon in border.shp
border = QgsVectorLayer('c:\\border.shp', 'border', 'ogr')
border.invertSelection()
borderFeatures = border.selectedFeatures()

#start loop at bottom left corner.
iface.mapCanvas().zoomToFullExtent()
mapRenderer = iface.mapCanvas().mapRenderer()
mapRect = mapRenderer.extent()
xStart = mapRect.xMinimum()
xEnd = mapRect.xMaximum()

yStart = mapRect.yMinimum()
yEnd = mapRect.yMaximum()

width = mapRenderer.width()
height = mapRenderer.height()

#printed raster will have (width * metersPerPixel) meters in X and
#(height * metersPerPixel) meters in Y
xDiff = width * metersPerPixel
yDiff = height * metersPerPixel

#Generate first rectangle
newRect = QgsRectangle(xStart, yStart, xStart + xDiff, yStart + yDiff)
mapRenderer.setExtent(newRect)

i = 0
#loop it!
while (newRect.yMinimum() < yEnd):
    while (newRect.xMinimum() < xEnd):
        i += 1
        #We don't want to print empty files
        intersects = 0
        for borderFeature in borderFeatures:
            if borderFeature.geometry().intersects(QgsRectangle(newRect.xMinimum(), newRect.yMinimum(), newRect.xMaximum(), newRect.yMaximum())):
                intersects = 1
        #if current rect intersects with country border, lets print it
        if intersects:
            #set up composition
            composition = QgsComposition(mapRenderer)
            composition.setPrintResolution(256)
            composition.setPaperSize(mapRenderer.width()/(mapRenderer.outputDpi()/25.4), mapRenderer.height()/(mapRenderer.outputDpi()/25.4))
            composition.setPlotStyle(QgsComposition.Print)
            dpi = composition.printResolution()
            dpmm = dpi / 25.4 #get dots per mm
            # add a map to the composition
            composerMap = QgsComposerMap(composition,0,0,composition.paperWidth(),composition.paperHeight())
            composition.addItem(composerMap)
            # create output image and initialize it
            image = QImage(QSize(width, height), QImage.Format_ARGB32) #output image size
            image.setDotsPerMeterX(dpmm * 1000) #mm to meters
            image.setDotsPerMeterY(dpmm * 1000) #mm to meters
            image.fill(0)
            # render the composition
            imagePainter = QPainter(image)
            sourceArea = QRectF(0, 0, composition.paperWidth(), composition.paperHeight() )
            targetArea = QRectF(0, 0, width, height)
            composition.render(imagePainter, targetArea, sourceArea)
            imagePainter.end()
            print "c:\\a" + str(i) + ".tif"
            image.save("c:\\a" + str(i) + ".tif" , "tif")
            createWorldFile("c:\\a" + str(i), metersPerPixel)
        #move along x
        newRect = QgsRectangle(newRect.xMaximum(), newRect.yMinimum(), newRect.xMaximum() + xDiff, newRect.yMaximum())
        mapRenderer.setExtent(newRect)
    #move along y
    newRect = QgsRectangle(xStart, newRect.yMaximum(), xStart + xDiff, newRect.yMaximum() + yDiff)
    mapRenderer.setExtent(newRect)
    print i

Note that this will print 16 meters per pixel.
It did not work as good as I wanted, so I had to return to mapnik in the end. But that will be in part III.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

%d bloggers like this: