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

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

Added TopOSM code and resources.

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