source: subversion/applications/rendering/toposm/toposm.py @ 34399

Last change on this file since 34399 was 25435, checked in by ahlzen, 9 years ago

Minor refactoring and cleanup of the TopOSM code.

File size: 12.2 KB
Line 
1#!/usr/bin/python
2
3"""toposm.py: Functions to control TopOSM rendering."""
4
5import sys, os, time, threading
6import numpy
7import multiprocessing
8from mapnik import Coord, Envelope
9from Queue import Queue
10from os import path
11from subprocess import call
12
13from env import *
14from coords import *
15from common import *
16import NED
17from JobManager import JobManager
18
19__author__      = "Lars Ahlzen and contributors"
20__copyright__   = "(c) Lars Ahlzen and contributors 2008-2011"
21__license__     = "GPLv2"
22
23
24# Print greeting/info on module load
25print "TopOSM. Using mapnik:", mapnik.mapnik_version()
26
27
28class RenderThread:
29    def __init__(self, q, maxz, threadNumber):
30        self.q = q
31        self.maxz = maxz
32        self.threadNumber = threadNumber
33        self.currentz = 0
34       
35    def init_zoomlevel(self, z):
36        self.currentz = z
37        self.tilesize = getTileSize(NTILES[z], True)
38        self.maps = {}
39        for mapName in MAPNIK_LAYERS:
40            self.maps[mapName] = mapnik.Map(self.tilesize, self.tilesize)
41            mapnik.load_map(self.maps[mapName], mapName + ".xml")
42
43    def runAndLog(self, message, function, args):
44        message = '[%02d] %s' % (self.threadNumber+1,  message)
45        console.printMessage(message)
46        try:
47            function(*args)       
48        except Exception as ex:
49            console.printMessage('Failed: ' + message)
50            errorLog.log('Failed: ' + message, ex)
51            raise
52
53    def renderTopoMetaTile(self, z, x, y, ntiles):
54        for mapname in ['hillshade', 'colormap']:   
55            msg = "Rendering %s %s %s %s" % (mapname, z, x, y)
56            self.runAndLog(msg, renderTopoMetaTile, (z, x, y, ntiles, mapname))
57           
58    def renderMapnikMetaTile(self, z, x, y, ntiles):
59        if (z != self.currentz):
60            self.init_zoomlevel(z)
61        for mapname in self.maps.keys():
62            msg = "Rendering %s %s %s %s" % (mapname, z, x, y)
63            self.runAndLog(msg, renderMapnikMetaTile, \
64                (z, x, y, ntiles, mapname, self.maps[mapname]))
65
66    def combineAndSliceTopoMetaTile(self, z, x, y, ntiles):
67        msg = "Combining topo tiles %s %s %s" % (z, x, y)
68        layers = (('color-relief', '.jpg'),)
69        self.runAndLog(msg, combineAndSlice, (z, x, y, ntiles, \
70            './combine-color-relief-tiles', layers))
71
72    def combineAndSliceMapnikMetaTile(self, z, x, y, ntiles):
73        msg = "Combining mapnik tiles %s %s %s" % (z, x, y)
74        layers = (('contours', '.png'), ('features', '.png'))
75        self.runAndLog(msg, combineAndSlice, (z, x, y, ntiles, \
76            './combine-mapnik-tiles', layers))
77   
78    def mergeTopoTiles(self, z, x, y, ntiles):
79        for dx in range(x*ntiles, (x+1)*ntiles):
80            for dy in range(y*ntiles, (y+1)*ntiles):
81                if not tileExists('color-relief', z, dx, xy, '.jpg'):
82                    msg = "Merging topo tiles %s %s %s" % (z, x, y)
83                    self.runAndLog(msg, mergeSubtiles, \
84                        (z, x, y, 'color-relief'))
85                   
86    def createJpegTiles(self, z, x, y, quality, ntiles):
87        for dx in range(x*ntiles, (x+1)*ntiles):
88            for dy in range(y*ntiles, (y+1)*ntiles):
89                if not tileExists('jpeg' + str(quality), z, dx, dy, '.jpg'):
90                    msg = "Creating JPEG%s %s %s %s" % (quality, z, dx, dy)
91                    self.runAndLog(msg, createJpegTile, (z, dx, dy, quality))
92           
93    def renderLoop(self):
94        self.currentz = 0
95        while True:
96            r = self.q.get()
97            if (r == None):
98                self.q.task_done()
99                break
100            self.render(*r)
101            self.q.task_done()
102
103    def render(self, action, z, x, y):
104        ntiles = NTILES[z]
105        if action == 'render':
106            self.renderMapnikMetaTile(z, x, y, ntiles)
107            self.combineAndSliceMapnikMetaTile(z, x, y, ntiles)
108            if z == self.maxz or \
109                not allConstituentTopoTilesExist(z, x, y, ntiles):
110                self.renderTopoMetaTile(z, x, y, ntiles)
111                self.combineAndSliceTopoMetaTile(z, x, y, ntiles)
112            else:
113                self.mergeTopoTiles(z, x, y, ntiles)
114            self.createJpegTiles(z, x, y, JPEG_QUALITY, ntiles)
115
116
117def getMetaTileDir(mapname, z):
118    return path.join(BASE_TILE_DIR, mapname, str(z))
119
120def getMetaTilePath(mapname, z, x, y, suffix = ".png"):
121    return path.join(getMetaTileDir(mapname, z), \
122        's' + str(x) + '_' + str(y) + suffix)
123
124def metaTileExists(mapname, z, x, y, suffix = ".png"):
125    return path.isfile(getMetaTilePath(mapname, z, x, y, suffix))
126
127def getTileDir(mapname, z, x):
128    return path.join(getMetaTileDir(mapname, z), str(x))
129
130def getTilePath(mapname, z, x, y, suffix = ".png"):
131    return path.join(getTileDir(mapname, z, x), str(y) + suffix)
132
133def tileExists(mapname, z, x, y, suffix = ".png"):
134    return path.isfile(getTilePath(mapname, z, x, y, suffix))
135
136def getTileSize(ntiles, includeBorder = True):
137    if includeBorder:
138        return TILE_SIZE * ntiles + 2 * BORDER_WIDTH
139    else:
140        return TILE_SIZE * ntiles
141       
142def allTilesExist(mapname, z, fromx, tox, fromy, toy, suffix = ".png"):
143    for x in range(fromx, tox+1):
144        for y in range(fromy, toy+1):
145            if not tileExists(mapname, z, x, y, suffix):
146                return False
147    return True
148   
149def allConstituentTopoTilesExist(z, x, y, ntiles):
150    subx = ntiles * x
151    suby = ntiles * y
152    return allTilesExist('color-relief', z+1, \
153        2*subx, 2*(subx+ntiles)-1, 2*suby, 2*(suby+ntiles)-1)
154
155def renderMapnikMetaTile(z, x, y, ntiles, mapname, map):
156    env = getMercTileEnv(z, x, y, ntiles, True)
157    tilesize = getTileSize(ntiles, True)
158    if not metaTileExists(mapname, z, x, y, '.png'):
159        map.zoom_to_box(env)
160        image = mapnik.Image(tilesize, tilesize)
161        mapnik.render(map, image)
162        view = image.view(BORDER_WIDTH, BORDER_WIDTH, tilesize, tilesize)
163        ensureDirExists(getMetaTileDir(mapname, z))
164        view.save(getMetaTilePath(mapname, z, x, y, '.png'))
165
166def renderTopoMetaTile(z, x, y, ntiles, mapname):
167    env = getMercTileEnv(z, x, y, ntiles, False)
168    envLL = getLLTileEnv(z, x, y, ntiles, False)
169    destdir = getMetaTileDir(mapname, z)
170    # NOTE: gdalwarp won't save as png directly, hence this "hack"
171    destTilePath = getMetaTilePath(mapname, z, x, y, '.tif')
172    finalTilePath = getMetaTilePath(mapname, z, x, y, '.png')
173    if os.path.isfile(finalTilePath):
174        pass
175    else:
176        ensureDirExists(destdir)
177        nedSlices = NED.getSlices(mapname, envLL, '.png')
178        tilesize = getTileSize(ntiles, False)
179        if len(nedSlices) > 0:
180            cmd = 'gdalwarp -q -t_srs "%s" -te %f %f %f %f -ts %d %d -r lanczos -dstnodata 255 ' % \
181                (MERCATOR_PROJECTION_DEF, env.minx, env.miny, env.maxx, env.maxy, tilesize, tilesize)
182            for slice in nedSlices:
183                cmd = cmd + '"' + slice + '" '
184            cmd = cmd + '"' + destTilePath + '"'
185            os.system(cmd)
186            # convert to png and remove tif (to conserve space)
187            cmd = 'convert "%s" "%s" && rm "%s"' % \
188                (destTilePath, finalTilePath, destTilePath)
189            os.system(cmd)
190        else:
191            return False
192    return True
193
194def sliceMetaTile(z, x, y, ntiles, mapname, destSuffix = '.png'):
195    srcTilePath = getMetaTilePath(mapname, z, x, y, '.png')
196    for dx in range(0, ntiles):
197        destDir = path.join(getMetaTileDir(mapname, z), str((x*ntiles)+dx))
198        ensureDirExists(destDir)
199        for dy in range(0, ntiles):
200            destfile = path.join(destDir, str((y*ntiles)+dy) + destSuffix)
201            if not path.isfile(destfile):
202                offsetx = dx * TILE_SIZE
203                offsety = dy * TILE_SIZE
204                cmd = 'convert "%s" -crop %dx%d+%d+%d +repage "%s"' % \
205                    (srcTilePath, TILE_SIZE, TILE_SIZE, offsetx, offsety, \
206                     path.join(destDir, str((y*ntiles)+dy) + destSuffix))
207                os.system(cmd)
208            else:
209                pass
210    tryRemove(srcTilePath)
211
212def combineAndSlice(z, x, y, ntiles, script, namesAndExts):
213    allExist = True
214    for n in namesAndExts:
215        if not allTilesExist(n[0], z, x*ntiles, (x+1)*ntiles-1, \
216            y*ntiles, (y+1)*ntiles-1, n[1]):
217            allExist = False
218            break
219    if not allExist:
220        cmd = "%s %s %d %d %d %d" % \
221            (script, BASE_TILE_DIR, z, x, y, getTileSize(ntiles, False))
222        os.system(cmd)
223        for n in namesAndExts:
224            sliceMetaTile(z, x, y, ntiles, n[0], n[1])
225
226def mergeSubtiles(z, x, y, mapname, suffix = '.jpg'):
227    """Merges (up to) four subtiles from the next higher
228    zoom level into one subtile at the specified location"""
229    cmd = 'convert -size 512x512 xc:white'
230    for dx in range(0,2):
231        for dy in range(0,2):
232            srcx = x*2 + dx
233            srcy = y*2 + dy
234            srcpath = getTilePath(mapname, z+1, srcx, srcy, suffix)
235            if os.path.isfile(srcpath):
236                cmd = cmd + ' "' + srcpath + '"'
237                cmd = cmd + ' -geometry +' + str(dx*256) + '+' + str(dy*256)
238                cmd = cmd + ' -composite'
239    cmd = cmd + ' -scale 256x256'
240    ensureDirExists(getTileDir(mapname, z, x))
241    destpath = getTilePath(mapname, z, x, y, suffix)
242    cmd = cmd + ' "' + destpath + '"'
243    os.system(cmd)
244
245def createJpegTile(z, x, y, quality):
246    colorreliefsrc = getTilePath('color-relief', z, x, y, '.jpg')
247    contourssrc = getTilePath('contours', z, x, y, '.png')
248    featuressrc = getTilePath('features', z, x, y, '.png')
249    desttile = getTilePath('jpeg' + str(quality), z, x, y, '.jpg')
250    if path.isfile(colorreliefsrc) and path.isfile(featuressrc):
251        ensureDirExists(path.dirname(desttile))       
252        # PIL generates internal errors opening the JPEG
253        # tiles so it's back to ImageMagick for now...
254        cmd = "convert " + colorreliefsrc;
255        if path.isfile(contourssrc):
256            cmd = cmd + " " + contourssrc + " -composite"
257        cmd = cmd + " " + featuressrc + " -composite"
258        cmd = cmd + " -quality " + str(quality) + " -strip " + desttile
259        os.system(cmd)
260 
261
262##### Public methods
263
264def prepareData(envLLs):
265    if not hasattr(envLLs, '__iter__'):
266        envLLs = (envLLs,)
267    manager = JobManager()
268    for envLL in envLLs:
269        tiles = NED.getTiles(envLL)       
270        for tile in tiles:
271            manager.addJob("Preparing %s" % (tile[0]), NED.prepDataFile, tile)
272    manager.finish()
273   
274    console.printMessage("Converting m to ft")
275    templ = 'echo "UPDATE %s SET height_ft = CAST(height * 3.28085 AS INT) ' \
276        + 'WHERE height_ft IS NULL;" | psql -q "%s"'
277    cmd = templ % (CONTOURS_TABLE, DATABASE)
278    os.system(cmd)
279   
280    # NOTE: This removes a LOT of artifacts around shorelines. Unfortunately,
281    # it also removes any legitimate 0ft contours (which may exist around
282    # areas below sea-level).
283    # TODO: Find a better workaround.
284    console.printMessage('Removing sea-level contours')
285    templ = 'echo "DELETE FROM %s WHERE height = 0;" | psql -q "%s"'
286    cmd = templ % (CONTOURS_TABLE, DATABASE)
287    os.system(cmd)
288
289def renderTiles(envLLs, minz, maxz):
290    if not hasattr(envLLs, '__iter__'):
291        envLLs = (envLLs,)
292    # Create mapnik rendering threads.
293    queue = Queue(32)
294    renderers = {}
295    for i in range(NUM_THREADS):
296        renderer = RenderThread(queue, maxz, i)
297        renderThread = threading.Thread(target=renderer.renderLoop)
298        renderThread.start()
299        renderers[i] = renderThread
300    # Queue up jobs. High-to-low zoom, so we can merge topo-tiles rather
301    # than render from scratch at lower zoom levels.
302    for envLL in envLLs:
303        for z in range(maxz, minz-1, -1):
304            ntiles = NTILES[z]
305            (fromx, tox, fromy, toy) = getTileRange(envLL, z, ntiles)
306            for x in range(fromx, tox+1):
307                for y in range(fromy, toy+1):
308                    queue.put(('render', z, x, y))
309            # Join threads after each completed level, since the
310            # color-relief layer on a lower zoom level depends on that of
311            # the next higher level being finished.
312            queue.join()
313    # Signal render threads to exit and join threads
314    for i in range(NUM_THREADS):
315        queue.put(None)
316    queue.join()
317    for i in range(NUM_THREADS):
318        renderers[i].join()       
319
Note: See TracBrowser for help on using the repository browser.