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

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

Work on highway shields. Minor style fixes.

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