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

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

Render thread bugfix (currentz.patch from Phil Gold)

File size: 26.1 KB
Line 
1#!/usr/bin/python
2import sys, os, time, threading
3import numpy
4#import Image # PIL
5from Queue import Queue
6from os import path
7from math import pi,cos,sin,log,exp,atan,ceil,floor
8from subprocess import call
9from mapnik import *
10
11
12# Check that the environment is set and import configuration
13if not 'TOPOSM_ENV_SET' in os.environ:
14    print "Error: TopOSM environment not set."
15    sys.exit(1)
16BASE_TILE_DIR = os.environ['BASE_TILE_DIR']
17CONTOURS_TABLE = os.environ['CONTOURS_TABLE']
18DATABASE = os.environ['DB_NAME']
19TEMPDIR = os.environ['TEMP_DIR']
20NED13DIR = os.environ['NED13_DIR']
21HILLSHADE = os.environ['HILLSHADE']
22COLORRELIEF = os.environ['COLORRELIEF']
23COLORFILE = os.environ['COLORFILE']
24NUM_THREADS = int(os.environ['RENDER_THREADS'])
25SUBTILE_SIZE = int(os.environ['SUBTILE_SIZE'])
26BORDER_WIDTH = int(os.environ['BORDER_WIDTH'])
27ERRORLOG = os.environ['ERROR_LOG']
28
29# constants
30CONTOUR_INTERVAL = 15.24 # 50 ft in meters
31MAPNIK_LAYERS = ['watermask', 'area', 'areansh', 'contourlines', 'contourlabels',
32                 'features-main', 'features-fill', 'labels', 'labels-nohalo']
33OSM_SRS = '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 ' + \
34    '+y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over'
35
36# Optimal supertile size (N x N subtiles) by zoom level.
37# A too low number is inefficient. A too high number uses
38# large amounts of memory and sometimes breaks the gdal tools.
39NTILES = {2:1, 3:1, 4:1, 5:1, 6:1, 7:1, 8:1, 9:1, 10:1,
40    11:2, 12:4, 13:6, 14:8, 15:10, 16:12, 17:12, 18:12 }
41
42# side length (in degrees) for slicing up NED files into more
43# manageable pieces during preprocessing
44STEP=0.5
45
46DEG_TO_RAD = pi/180
47RAD_TO_DEG = 180/pi
48
49def minmax (a,b,c):
50    a = max(a,b)
51    a = min(a,c)
52    return a
53
54class GoogleProjection:
55    def __init__(self,levels=20):
56        self.Bc = []
57        self.Cc = []
58        self.zc = []
59        self.Ac = []
60        c = 256
61        for d in range(0,levels):
62            e = c/2;
63            self.Bc.append(c/360.0)
64            self.Cc.append(c/(2 * pi))
65            self.zc.append((e,e))
66            self.Ac.append(c)
67            c *= 2
68               
69    def fromLLtoPixel(self,ll,zoom):
70         d = self.zc[zoom]
71         e = round(d[0] + ll[0] * self.Bc[zoom])
72         f = minmax(sin(DEG_TO_RAD * ll[1]),-0.9999,0.9999)
73         g = round(d[1] + 0.5*log((1+f)/(1-f))*-self.Cc[zoom])
74         return (e,g)
75     
76    def fromPixelToLL(self,px,zoom):
77         e = self.zc[zoom]
78         f = (px[0] - e[0])/self.Bc[zoom]
79         g = (px[1] - e[1])/-self.Cc[zoom]
80         h = RAD_TO_DEG * ( 2 * atan(exp(g)) - 0.5 * pi)
81         return (f,h)
82
83GOOGLE_PROJECTION = GoogleProjection()
84
85LATLONG_PROJECTION_DEF = "+proj=latlong"
86LATLONG_PROJECTION = Projection(LATLONG_PROJECTION_DEF)
87
88MERCATOR_PROJECTION_DEF = "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over"
89MERCATOR_PROJECTION = Projection(MERCATOR_PROJECTION_DEF)
90
91# global locking object for file system ops (e.g. creating directories)
92fslock = threading.Lock()
93
94# thread-safe console printing access
95printlock = threading.Lock()
96def print_message(message):
97    printlock.acquire()
98    print message
99    printlock.release()
100
101# error log
102errorloglock = threading.Lock()
103def log_error(message, exception = None):
104    errorloglock.acquire()
105    try:
106        file = open(ERRORLOG, 'a')
107        file.write(message)
108        if exception != None:
109            file.write(" %s %s\n" % \
110                (str(exception), exception.message))
111        file.write('\n')
112        file.close()
113    except:
114        print "Failed to write to the error log:"
115        print "%s %s" % (sys.exc_info()[0], sys.exc_info()[1])
116    finally:
117        errorloglock.release()
118
119class RenderThread:
120    def __init__(self, q, maxz, threadNumber):
121        self.q = q
122        self.maxz = maxz
123        self.threadNumber = threadNumber
124        self.currentz = 0
125       
126    def init_zoomlevel(self, z):
127        self.currentz = z
128        self.ntiles = NTILES[z]
129        self.tilesize = getTileSize(self.ntiles, True)
130        self.maps = {}
131        for mapName in MAPNIK_LAYERS:
132            self.maps[mapName] = Map(self.tilesize, self.tilesize)
133            load_map(self.maps[mapName], mapName + ".xml")
134
135    def print_message(self, message):
136        print_message("[%02d] %s" % (self.threadNumber+1, message))
137
138    def prepare_data(self, basename, env):
139        self.print_message("Preparing %s at %s" % (basename, env))
140        try:
141            prep_ned13_data_file(basename, env)
142        except Exception as ex:
143            message = "Failed prepare_data %s %s" % (basename, env)
144            self.print_message(message)
145            log_error(message, ex) 
146
147    def render_topo_tiles(self, z, x, y):
148        ntiles = NTILES[z]
149        for mapname in ['hillshade', 'colormap']:   
150            self.print_message("Rendering %s %s %s %s" % (mapname, z, x, y))
151            try:
152                render_geotiff_tile(z, x, y, ntiles, mapname)
153            except Exception as ex:
154                message = "Failed render_topo_tiles %s %s %s %s" % (mapname, z, x, y)
155                self.print_message(message)
156                log_error(message, ex)
157
158    def render_mapnik_tiles(self, z, x, y):
159        if (z != self.currentz):
160            self.init_zoomlevel(z)
161        for mapname in self.maps.keys():
162            self.print_message("Rendering %s %s %s %s" % (mapname, z, x, y))
163            try:
164                render_mapnik_tile(z, x, y, self.ntiles, mapname, self.maps[mapname])
165            except Exception as ex:
166                message = "Failed render_mapnik_tiles %s %s %s %s" % (mapname, z, x, y)
167                self.print_message(message)
168                log_error(message, ex)
169
170    def combine_topo_tiles(self, z, x, y):
171        ntiles = NTILES[z]
172        if allSubtilesExist('color-relief', z, x, x+ntiles-1, y, y+ntiles-1, '.jpg'):
173            self.print_message("All topo subtiles exist. Skipping.")
174            return
175        try:
176            self.print_message("Combining topo tiles at %s %s %s" % (z, x, y))
177            command = "./combine-color-relief-tiles %s %d %d %d %d" % \
178                (BASE_TILE_DIR, z, x, y, getTileSize(ntiles, False));
179            call(command, shell=True)
180            slice_tile(z, x, y, ntiles, 'color-relief', '.jpg')
181        except Exception as ex:
182            message = "Failed combine_topo_tiles %s %s %s" % (z, x, y)
183            self.print_message(message)
184            log_error(message, ex)
185
186    def combine_mapnik_tiles(self, z, x, y):
187        ntiles = NTILES[z]
188        if allSubtilesExist('contours', z, x, x+ntiles-1, y, y+ntiles-1, '.png') and \
189            allSubtilesExist('features', z, x, x+ntiles-1, y, y+ntiles-1, '.png'):
190            self.print_message("All mapnik subtiles exist. Skipping.")
191            return
192        try:
193            self.print_message("Combining mapnik tiles at %s %s %s" % (z, x, y))
194            command = "./combine-mapnik-tiles %s %d %d %d %d" % \
195                (BASE_TILE_DIR, z, x, y, getTileSize(ntiles, False));
196            call(command, shell=True)
197            slice_tile(z, x, y, ntiles, 'contours', '.png')
198            slice_tile(z, x, y, ntiles, 'features', '.png')
199        except Exception as ex:
200            message = "Failed combine_mapnik_tiles %s %s %s" % (z, x, y)
201            self.print_message(message)
202            log_error(message, ex)
203
204    def merge_topo_subtiles(self, z, x, y):
205        ntiles = NTILES[z]
206        for dx in range(x, x+ntiles):
207            for dy in range(y, y+ntiles):
208                if subtileExists('color-relief', z, dx, dy, '.jpg'):
209                    self.print_message('Subtile exists. Skipping.')
210                    continue
211                try:
212                    self.print_message("Merging subtiles at %s %s %s" % (z, dx, dy))
213                    merge_subtiles(z, dx, dy, 'color-relief', '.jpg')
214                except Exception as ex:
215                    message = "Failed merge_subtiles %s %s %s" % (z, dx, dy)
216                    self.print_message(message)
217                    log_error(message, ex)
218           
219    def render_loop(self):
220        self.currentz = 0
221        while True:
222            r = self.q.get()
223            if (r == None):
224                self.q.task_done()
225                break
226            (action, z, x, y) = r
227            self.render(action, z, x, y)
228            self.q.task_done()
229
230    def render(self, action, z, x, y):
231        if action == 'prepare_data':
232            self.prepare_data(x, y)
233        elif action == 'render':
234            # NOTE: Mapnik tiles must be rendered before topo
235            if allFinalMapnikSubtilesExist(z, x, y):
236                self.print_message('Mapnik tiles exist. Skipping.')
237            else:
238                self.render_mapnik_tiles(z, x, y)
239                self.combine_mapnik_tiles(z, x, y)
240            if allFinalTopoSubtilesExist(z, x, y):
241                self.print_message('Topo exist at %d %d %d.' % (z, x, y))
242            else:
243                if z == self.maxz:
244                    self.render_topo_tiles(z, x, y)
245                    self.combine_topo_tiles(z, x, y)
246                else:
247                    self.merge_topo_subtiles(z, x, y)     
248
249def getTileDir(mapname, z):
250    return path.join(BASE_TILE_DIR, mapname, str(z))
251
252def getTilePath(mapname, z, x, y, suffix = ".png"):
253    return path.join(getTileDir(mapname, z), 's' + str(x) + '_' + str(y) + suffix)
254
255def getSubtileDir(mapname, z, x):
256    return path.join(getTileDir(mapname, z), str(x))
257
258def getSubtilePath(mapname, z, x, y, suffix = ".png"):
259    return path.join(getSubtileDir(mapname, z, x), str(y) + suffix)
260
261def getTileSize(ntiles, includeBorder = True):
262    if includeBorder:
263        return SUBTILE_SIZE * ntiles + 2 * BORDER_WIDTH;
264    else:
265        return SUBTILE_SIZE * ntiles;
266
267def tileExists(mapname, z, x, y, suffix = ".png"):
268    return path.isfile(getTilePath(mapname, z, x, y, suffix))
269
270def subtileExists(mapname, z, x, y, suffix = ".png"):
271    return path.isfile(getSubtilePath(mapname, z, x, y, suffix))
272   
273def allSubtilesExist(mapname, z, fromx, tox, fromy, toy, suffix = ".png"):
274    for x in range(fromx, tox+1):
275        for y in range(fromy, toy+1):
276            if not subtileExists(mapname, z, x, y, suffix):
277                return False
278    return True
279
280def allFinalTopoSubtilesExist(z, x, y):
281    ntiles = NTILES[z]
282    return allSubtilesExist('color-relief', z, x, x+ntiles-1, y, y+ntiles-1, '.jpg')
283
284def allFinalMapnikSubtilesExist(z, x, y):
285    ntiles = NTILES[z]
286    return \
287        allSubtilesExist('contours', z, x, x+ntiles-1, y, y+ntiles-1, '.png') and \
288        allSubtilesExist('features', z, x, x+ntiles-1, y, y+ntiles-1, '.png')
289
290def ensureDirExists(path):
291    fslock.acquire()
292    try:
293        if not os.path.isdir(path):
294            os.makedirs(path)
295    finally:
296        fslock.release()
297
298def tryRemove(filename):
299    fslock.acquire()
300    try:
301            if path.isfile(filename):
302                    os.remove(filename)
303    finally:
304        fslock.release()
305
306def writeEmpty(filename):
307        "Overwrites the specified filename with a new empty file."
308        fslock.acquire()
309        try:
310            open(filename, 'w').close();
311        finally:
312            fslock.release()
313       
314def render_mapnik_tile(z, x, y, ntiles, mapname, map):
315    bboxMerc = get_bbox_merc(z, x, y, ntiles, True)
316    to_x = x + ntiles
317    to_y = y + ntiles
318    if tileExists(mapname, z, x, y, '.png'):
319        pass
320    else:
321        map.zoom_to_box(bboxMerc)
322        tilesize = getTileSize(ntiles, True)
323        image = Image(tilesize, tilesize)
324        render(map, image)
325        view = image.view(BORDER_WIDTH, BORDER_WIDTH, tilesize, tilesize)
326        ensureDirExists(getTileDir(mapname, z))
327        view.save(getTilePath(mapname, z, x, y, '.png'))
328
329def get_bbox_ll(z, x, y, ntiles, includeBorder = True):
330    "Returns the lat/lon bbox for the specified tile"
331    border = 0
332    if includeBorder:
333        border = BORDER_WIDTH
334    p0 = GOOGLE_PROJECTION.fromPixelToLL(
335        (x * SUBTILE_SIZE - border,
336         (y + ntiles) * SUBTILE_SIZE + border), z)
337    p1 = GOOGLE_PROJECTION.fromPixelToLL(
338        ((x + ntiles) * SUBTILE_SIZE + border,
339         y * SUBTILE_SIZE - border), z)
340    return Envelope(p0[0], p0[1], p1[0], p1[1])
341                                       
342def get_bbox_merc(z, x, y, ntiles, includeBorder = True):
343    "Returns the mercator bbox for the specified tile"
344    border = 0
345    if includeBorder:
346        border = BORDER_WIDTH
347    p0 = GOOGLE_PROJECTION.fromPixelToLL(
348        (x * SUBTILE_SIZE - border,
349         (y + ntiles) * SUBTILE_SIZE + border), z)
350    p1 = GOOGLE_PROJECTION.fromPixelToLL(
351        ((x + ntiles) * SUBTILE_SIZE + border,
352         y * SUBTILE_SIZE - border), z)
353    c0 = MERCATOR_PROJECTION.forward(Coord(p0[0], p0[1]))
354    c1 = MERCATOR_PROJECTION.forward(Coord(p1[0], p1[1]))           
355    return Envelope(c0.x, c0.y, c1.x, c1.y)
356
357def get_tile_from_ll(lat, lon, z):
358    "Returns the subtile number (x, y) at the specified lat/lon/zoom"
359    (px,py) = GOOGLE_PROJECTION.fromLLtoPixel((lon,lat),z)
360    return ((int)(px/SUBTILE_SIZE), (int)(py/SUBTILE_SIZE))
361
362def get_ned13_tilepath(basename):
363    return os.path.join(NED13DIR, basename, 'grd' + basename + '_13', 'w001001.adf')
364
365def get_ned13_tiles(envLL):
366    "Gets the (basename, envelope) of all (existing) 1/3 NED tiles to cover the specified envelope"
367    fromx = int(floor(envLL.minx))
368    tox = int(floor(envLL.maxx))
369    fromy = int(ceil(envLL.miny))
370    toy = int(ceil(envLL.maxy))
371    tiles = []
372    for y in range(fromy, toy+1):
373        for x in range(fromx, tox+1):
374            basename = 'n%02dw%03d' % (y, -x)
375            tilepath = get_ned13_tilepath(basename)
376            if os.path.isfile(tilepath):
377                tiles.append((basename, Envelope(x, y-1, x+1, y)))
378    return tiles
379
380def get_ned13_slice(prefix, x, y, suffix = '.tif'):
381    filename = prefix + '_%.1f_%.1f%s' % (float(x), float(y), suffix)
382    return path.join(TEMPDIR, filename)
383
384def get_ned13_slices(prefix, envLL, suffix = '.tif'):
385    fromx = floor(envLL.minx/STEP)*STEP
386    fromy = floor(envLL.miny/STEP)*STEP
387    tox = ceil(envLL.maxx/STEP)*STEP
388    toy = ceil(envLL.maxy/STEP)*STEP
389    slices = []
390    for y in numpy.arange(fromy, toy, STEP):
391        for x in numpy.arange(fromx, tox, STEP):
392            slicefile = get_ned13_slice(prefix, x, y, suffix)
393            if os.path.isfile(slicefile):
394                slices.append(slicefile)
395    return slices
396
397def get_ned13_tilecoords(lat, lon, step = 1):
398    return (int(ceil(lat/float(step)))*float(step), int(floor(lon/float(step)))*float(step))
399
400def get_ned13_tilename(lat, lon, step = 1):
401    (y, x) = get_ned13_tilecoords(lat, lon, step)
402    if step == 1:
403        return 'n%02dw%03d' % (y, -x)
404    return 'n%02.2fw%03.2f' % (y, -x)
405                       
406def prep_ned13_data_file(basename, env):
407    print 'Preparing NED 1/3" data file:', basename
408    print '  Converting to GeoTIFF...'
409    sourcepath = get_ned13_tilepath(basename)
410    geotiff = path.join(TEMPDIR, basename + '.tif')
411    if not path.isfile(geotiff):
412        cmd = 'gdal_translate "%s" "%s"' % (sourcepath, geotiff)
413        call(cmd, shell = True)
414
415    print '  Generating contour lines...'
416    # split the GeoTIFF, since it's often too large otherwise
417    for y in numpy.arange(env.miny, env.maxy, STEP):
418        for x in numpy.arange(env.minx, env.maxx, STEP):
419            print '  Cutting geotiff slice...'
420            nedslice = get_ned13_slice('ned', x, y)
421            if not path.isfile(nedslice):
422                cmd = 'gdalwarp -q -te %f %f %f %f "%s" "%s"' % \
423                    (x, y, x+STEP, y+STEP, geotiff, nedslice)
424                call(cmd, shell=True)
425           
426            print '  Generating contour lines...'
427            contourbasefile = path.join(TEMPDIR, 'contours_' + str(x) + '_' + str(y))
428            contourfile = contourbasefile + '.shp'
429            if not path.isfile(contourfile):
430                cmd = 'gdal_contour -i %f -snodata 32767 -a height "%s" "%s"' % \
431                    (CONTOUR_INTERVAL, nedslice, contourfile)
432                call(cmd, shell=True)
433
434                print '  Importing contour lines...'
435                # NOTE: this assumes that the table is already set up
436                cmd = 'shp2pgsql -a -g way "%s" "%s" | psql -q "%s"' % \
437                    (contourfile, CONTOURS_TABLE, DATABASE)
438                call(cmd, shell=True)
439               
440                # Clear contents (but keep file to prevent us from importing
441                # these contours again).
442                writeEmpty(contourfile)
443               
444                # Remove shape index and attribute files.
445                tryRemove(contourbasefile + ".shx")
446                tryRemove(contourbasefile + ".dbf")
447
448            print '  Generating hillshade slice...'
449            hillshadeslice = get_ned13_slice('hillshade', x, y)
450            hillshadeslicePng = get_ned13_slice('hillshade', x, y, '.png')
451            if not path.isfile(hillshadeslicePng):
452                cmd = '"%s" "%s" "%s" -z 0.00001' % \
453                    (HILLSHADE, nedslice, hillshadeslice)
454                call(cmd, shell = True)
455                # convert to PNG + world file to save space
456                cmd = 'gdal_translate -quiet -of PNG -co WORLDFILE=YES "%s" "%s"' % \
457                    (hillshadeslice, hillshadeslicePng)
458                call(cmd, shell = True)
459                tryRemove(hillshadeslice)
460
461            print '  Generating colormap slice...'
462            colormapslice = get_ned13_slice('colormap', x, y)
463            colormapslicePng = get_ned13_slice('colormap', x, y, '.png')
464            if not path.isfile(colormapslicePng):
465                cmd = '"%s" "%s" "%s" "%s"' % \
466                    (COLORRELIEF, nedslice, COLORFILE, colormapslice)
467                call(cmd, shell = True)
468                # convert to PNG + world file to save space
469                cmd = 'gdal_translate -quiet -of PNG -co WORLDFILE=YES "%s" "%s"' % \
470                    (colormapslice, colormapslicePng)
471                call(cmd, shell = True)
472                tryRemove(colormapslice)
473
474            writeEmpty(nedslice) # don't need the raw slice anymore.
475               
476    writeEmpty(geotiff) # done with this GeoTIFF.
477
478def pad_envelope(envLL, z, ntiles):
479    """Extends the specified area with the width/height of a subtile toward
480    increasing x/y, plus one unit of padding on all sides. Used to ensure that
481    enough data is included for rendering a set of tiles."""
482    (xmin, ymax) = GOOGLE_PROJECTION.fromLLtoPixel((envLL.minx, envLL.miny), z)
483    (xmax, ymin) = GOOGLE_PROJECTION.fromLLtoPixel((envLL.maxx, envLL.maxy), z)
484    xmin -= BORDER_WIDTH
485    ymin -= BORDER_WIDTH
486    xmax += BORDER_WIDTH + SUBTILE_SIZE
487    ymax += BORDER_WIDTH + SUBTILE_SIZE
488    (lonmin, latmin) = GOOGLE_PROJECTION.fromPixelToLL((xmin, ymax), z)
489    (lonmax, latmax) = GOOGLE_PROJECTION.fromPixelToLL((xmax, ymin), z)
490    return Envelope(lonmin, latmin, lonmax, latmax)
491
492# TODO: rename method... e.g. render_topo_tile
493def render_geotiff_tile(z, x, y, ntiles, mapname):
494    bbox = get_bbox_merc(z, x, y, ntiles, False)
495    bboxLL = get_bbox_ll(z, x, y, ntiles, False)
496    destdir = getTileDir(mapname, z)
497    # NOTE: gdalwarp won't save as png directly, hence this "hack"
498    destTilePath = getTilePath(mapname, z, x, y, '.tif')
499    finalTilePath = getTilePath(mapname, z, x, y, '.png')
500    if os.path.isfile(destTilePath):
501        pass
502    else:
503        ensureDirExists(destdir)
504        nedSlices = get_ned13_slices(mapname, bboxLL, '.png')
505        tilesize = getTileSize(ntiles, False)
506        if len(nedSlices) > 0:
507            cmd = 'gdalwarp -q -t_srs "%s" -te %f %f %f %f -ts %d %d -r lanczos -dstnodata 255 ' % \
508                (OSM_SRS, bbox.minx, bbox.miny, bbox.maxx, bbox.maxy, tilesize, tilesize)
509            for slice in nedSlices:
510                cmd = cmd + '"' + slice + '" '
511            cmd = cmd + '"' + destTilePath + '"'
512            print cmd
513            call(cmd, shell=True)
514            # convert to png and remove tif (to conserve space)
515            cmd = 'convert "%s" "%s" && rm "%s"' % (destTilePath, finalTilePath, destTilePath)
516            call(cmd, shell=True)
517        else:
518            return False
519    return True
520
521def slice_tile(z, x, y, ntiles, mapname, destSuffix = '.png'):
522    srcTilePath = getTilePath(mapname, z, x, y, '.png')
523    for dx in range(0, ntiles):
524        destDir = path.join(getTileDir(mapname, z), str(x+dx))
525        ensureDirExists(destDir)
526        for dy in range(0, ntiles):
527            destfile = path.join(destDir, str(y+dy) + destSuffix)
528            if not path.isfile(destfile):
529                offsetx = dx * SUBTILE_SIZE
530                offsety = dy * SUBTILE_SIZE
531                cmd = 'convert "%s" -crop %dx%d+%d+%d +repage "%s"' % \
532                    (srcTilePath, SUBTILE_SIZE, SUBTILE_SIZE, offsetx, offsety,
533                     path.join(destDir, str(y+dy) + destSuffix))
534                call(cmd, shell=True)
535            else:
536                pass
537    tryRemove(srcTilePath)
538
539
540def merge_subtiles(z, x, y, mapname, suffix = '.jpg'):
541    """Merges (up to) four subtiles from the next higher
542    zoom level into one subtile at the specified location"""
543#    print 'Merging tile ', z, x, y, mapname
544    cmd = 'convert -size 512x512 xc:white'
545    for dx in range(0,2):
546        for dy in range(0,2):
547            srcx = x*2 + dx
548            srcy = y*2 + dy
549            srcpath = getSubtilePath(mapname, z+1, srcx, srcy, suffix)
550            if os.path.isfile(srcpath):
551                cmd = cmd + ' "' + srcpath + '"'
552                cmd = cmd + ' -geometry +' + str(dx*256) + '+' + str(dy*256)
553                cmd = cmd + ' -composite'
554    cmd = cmd + ' -scale 256x256'
555    ensureDirExists(getSubtileDir(mapname, z, x))
556    destpath = getSubtilePath(mapname, z, x, y, suffix)
557    cmd = cmd + ' "' + destpath + '"'
558    call(cmd, shell=True)
559
560def create_jpeg_tile(z, x, y, quality):
561    print 'Creating jpeg tile at', z, x, y, '...',
562    colorreliefsrc = getSubtilePath('color-relief', z, x, y, '.jpg')
563    contourssrc = getSubtilePath('contours', z, x, y, '.png')
564    featuressrc = getSubtilePath('features', z, x, y, '.png')
565    desttile = getSubtilePath('jpeg' + str(quality), z, x, y, '.jpg')
566    if path.isfile(colorreliefsrc) and path.isfile(featuressrc):
567        ensureDirExists(path.dirname(desttile))
568
569        #tile = Image.open(colorreliefsrc)
570        #if path.isfile(contourssrc):
571        #    contours = Image.open(contourssrc)
572        #    tile.paste(contours, None)
573        #features = Image.open(featuressrc)
574        #tile.paste(features, None)
575        #ensure_dir_exists(path.dirname(desttile))
576        #tile.save(desttile, 'JPEG', quality=quality, optimize=True)
577       
578        # PIL generates internal errors opening the JPEG
579        # tiles so it's back to ImageMagick for now...
580        cmd = "convert " + colorreliefsrc;
581        if path.isfile(contourssrc):
582            cmd = cmd + " " + contourssrc + " -composite"
583        cmd = cmd + " " + featuressrc + " -composite"
584        cmd = cmd + " -quality " + str(quality) + " " + desttile
585        call(cmd, shell=True)
586       
587        print 'done.'
588    else:
589        print 'source tiles not found. skipping.', colorreliefsrc, featuressrc
590
591
592# public methods
593
594def prepare_data_single(envLL, minz, maxz):
595    ntiles = NTILES[maxz]
596    envLL = pad_envelope(envLL, minz, ntiles)
597    tiles = get_ned13_tiles(envLL)
598    for tile in tiles:
599        print tile
600        prep_ned13_data_file(tile[0], tile[1])
601    print '  Converting meters to feet...'
602    cmd = 'echo "UPDATE %s SET height_ft = CAST(height * 3.28085 AS INT) WHERE height_ft IS NULL;" | psql -q "%s"' % (CONTOURS_TABLE, DATABASE)
603    call(cmd, shell=True)
604
605def prepare_data(envLL, minz, maxz):
606    ntiles = NTILES[maxz]
607    envLL = pad_envelope(envLL, minz, ntiles)
608    tiles = get_ned13_tiles(envLL)
609    queue = Queue(32)
610    renderers = {}
611    for i in range(NUM_THREADS):
612        renderer = RenderThread(queue, maxz, i)
613        render_thread = threading.Thread(target=renderer.render_loop)
614        render_thread.start()
615        renderers[i] = render_thread
616    for tile in tiles:
617        queue.put(('prepare_data', minz, tile[0], tile[1]))
618    for i in range(NUM_THREADS):
619        queue.put(None)
620    queue.join()
621    for i in range(NUM_THREADS):
622        renderers[i].join()
623    print 'Converting meters to feet...'
624    cmd = 'echo "UPDATE %s SET height_ft = CAST(height * 3.28085 AS INT) WHERE height_ft IS NULL;" | psql -q "%s"' % (CONTOURS_TABLE, DATABASE)
625    call(cmd, shell=True)
626    print 'Done.'
627
628def render_tiles(envLL, minz, maxz):
629    print 'Using mapnik version:', mapnik_version()
630    # Create mapnik rendering threads.
631    queue = Queue(32)
632    renderers = {}
633    for i in range(NUM_THREADS):
634        renderer = RenderThread(queue, maxz, i)
635        render_thread = threading.Thread(target=renderer.render_loop)
636        render_thread.start()
637        renderers[i] = render_thread
638    # Queue up jobs. High-to-low zoom, so we can merge topo-subtiles rather
639    # than render from scratch at lower zoom levels.
640    for z in range(maxz, minz-1, -1):
641        ntiles = NTILES[z]
642        (fromx, fromy) = get_tile_from_ll(envLL.maxy, envLL.minx, z)
643        (tox, toy) = get_tile_from_ll(envLL.miny, envLL.maxx, z)
644        for x in range(fromx, tox+1, ntiles):
645            for y in range(fromy, toy+1, ntiles):
646                queue.put(('render', z, x, y))
647        # Join threads after each completed level, because
648        # the color-relief layer of a lower zoom level depends
649        # on that of the next higher level being finished.
650        queue.join()
651    # Signal render threads to exit and join threads
652    for i in range(NUM_THREADS):
653        queue.put(None)
654    queue.join()
655    for i in range(NUM_THREADS):
656        renderers[i].join()
657
658def render_tiles_single(envLL, minz, maxz):
659    print 'Using mapnik version:', mapnik_version()
660    renderer = RenderThread(None, maxz, 1)
661    for z in range(maxz, minz-1, -1):
662        ntiles = NTILES[z]
663        (fromx, fromy) = get_tile_from_ll(envLL.maxy, envLL.minx, z)
664        (tox, toy) = get_tile_from_ll(envLL.miny, envLL.maxx, z)
665        for x in range(fromx, tox+1, ntiles):
666            for y in range(fromy, toy+1, ntiles):
667                renderer.render('render', z, x, y)
668
669def create_jpeg_tiles(envLL, minz, maxz, quality):
670    for z in range(minz, maxz+1):
671        (fromx, fromy) = get_tile_from_ll(envLL.maxy, envLL.minx, z)
672        (tox, toy) = get_tile_from_ll(envLL.miny, envLL.maxx, z)
673        for x in range(fromx, tox+1):
674            for y in range(fromy, toy+1):
675                create_jpeg_tile(z, x, y, quality)
676
Note: See TracBrowser for help on using the repository browser.