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

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

Added support for creating composite JPEG tiles. Minor cleanups.

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