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

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

Applied Phil Gold's metatiles and generate-all-color-relief patches.
Added power line rendering support (partly based on Phil Gold's patch).

File size: 26.7 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 or not allConstituentTopoSubtilesExist(z, x, y):
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 allConstituentTopoSubtilesExist(z, x, y):
281    ntiles = NTILES[z]
282    return allSubtilesExist('color-relief', z+1, 2*x, 2*(x+ntiles)-1, 2*y, 2*(y+ntiles)-1, '.jpg')
283
284def allFinalTopoSubtilesExist(z, x, y):
285    ntiles = NTILES[z]
286    return allSubtilesExist('color-relief', z, x, x+ntiles-1, y, y+ntiles-1, '.jpg')
287
288def allFinalMapnikSubtilesExist(z, x, y):
289    ntiles = NTILES[z]
290    return \
291        allSubtilesExist('contours', z, x, x+ntiles-1, y, y+ntiles-1, '.png') and \
292        allSubtilesExist('features', z, x, x+ntiles-1, y, y+ntiles-1, '.png')
293
294def ensureDirExists(path):
295    fslock.acquire()
296    try:
297        if not os.path.isdir(path):
298            os.makedirs(path)
299    finally:
300        fslock.release()
301
302def tryRemove(filename):
303    fslock.acquire()
304    try:
305            if path.isfile(filename):
306                    os.remove(filename)
307    finally:
308        fslock.release()
309
310def writeEmpty(filename):
311        "Overwrites the specified filename with a new empty file."
312        fslock.acquire()
313        try:
314            open(filename, 'w').close();
315        finally:
316            fslock.release()
317       
318def render_mapnik_tile(z, x, y, ntiles, mapname, map):
319    bboxMerc = get_bbox_merc(z, x, y, ntiles, True)
320    to_x = x + ntiles
321    to_y = y + ntiles
322    if tileExists(mapname, z, x, y, '.png'):
323        pass
324    else:
325        map.zoom_to_box(bboxMerc)
326        tilesize = getTileSize(ntiles, True)
327        image = Image(tilesize, tilesize)
328        render(map, image)
329        view = image.view(BORDER_WIDTH, BORDER_WIDTH, tilesize, tilesize)
330        ensureDirExists(getTileDir(mapname, z))
331        view.save(getTilePath(mapname, z, x, y, '.png'))
332
333def get_bbox_ll(z, x, y, ntiles, includeBorder = True):
334    "Returns the lat/lon bbox for the specified tile"
335    border = 0
336    if includeBorder:
337        border = BORDER_WIDTH
338    p0 = GOOGLE_PROJECTION.fromPixelToLL(
339        (x * SUBTILE_SIZE - border,
340         (y + ntiles) * SUBTILE_SIZE + border), z)
341    p1 = GOOGLE_PROJECTION.fromPixelToLL(
342        ((x + ntiles) * SUBTILE_SIZE + border,
343         y * SUBTILE_SIZE - border), z)
344    return Envelope(p0[0], p0[1], p1[0], p1[1])
345                                       
346def get_bbox_merc(z, x, y, ntiles, includeBorder = True):
347    "Returns the mercator bbox for the specified tile"
348    border = 0
349    if includeBorder:
350        border = BORDER_WIDTH
351    p0 = GOOGLE_PROJECTION.fromPixelToLL(
352        (x * SUBTILE_SIZE - border,
353         (y + ntiles) * SUBTILE_SIZE + border), z)
354    p1 = GOOGLE_PROJECTION.fromPixelToLL(
355        ((x + ntiles) * SUBTILE_SIZE + border,
356         y * SUBTILE_SIZE - border), z)
357    c0 = MERCATOR_PROJECTION.forward(Coord(p0[0], p0[1]))
358    c1 = MERCATOR_PROJECTION.forward(Coord(p1[0], p1[1]))           
359    return Envelope(c0.x, c0.y, c1.x, c1.y)
360
361def get_tile_from_ll(lat, lon, z):
362    "Returns the subtile number (x, y) at the specified lat/lon/zoom"
363    (px,py) = GOOGLE_PROJECTION.fromLLtoPixel((lon,lat),z)
364    return ((int)(px/SUBTILE_SIZE), (int)(py/SUBTILE_SIZE))
365
366def get_ned13_tilepath(basename):
367    return os.path.join(NED13DIR, basename, 'grd' + basename + '_13', 'w001001.adf')
368
369def get_ned13_tiles(envLL):
370    "Gets the (basename, envelope) of all (existing) 1/3 NED tiles to cover the specified envelope"
371    fromx = int(floor(envLL.minx))
372    tox = int(floor(envLL.maxx))
373    fromy = int(ceil(envLL.miny))
374    toy = int(ceil(envLL.maxy))
375    tiles = []
376    for y in range(fromy, toy+1):
377        for x in range(fromx, tox+1):
378            basename = 'n%02dw%03d' % (y, -x)
379            tilepath = get_ned13_tilepath(basename)
380            if os.path.isfile(tilepath):
381                tiles.append((basename, Envelope(x, y-1, x+1, y)))
382    return tiles
383
384def get_ned13_slice(prefix, x, y, suffix = '.tif'):
385    filename = prefix + '_%.1f_%.1f%s' % (float(x), float(y), suffix)
386    return path.join(TEMPDIR, filename)
387
388def get_ned13_slices(prefix, envLL, suffix = '.tif'):
389    fromx = floor(envLL.minx/STEP)*STEP
390    fromy = floor(envLL.miny/STEP)*STEP
391    tox = ceil(envLL.maxx/STEP)*STEP
392    toy = ceil(envLL.maxy/STEP)*STEP
393    slices = []
394    for y in numpy.arange(fromy, toy, STEP):
395        for x in numpy.arange(fromx, tox, STEP):
396            slicefile = get_ned13_slice(prefix, x, y, suffix)
397            if os.path.isfile(slicefile):
398                slices.append(slicefile)
399    return slices
400
401def get_ned13_tilecoords(lat, lon, step = 1):
402    return (int(ceil(lat/float(step)))*float(step), int(floor(lon/float(step)))*float(step))
403
404def get_ned13_tilename(lat, lon, step = 1):
405    (y, x) = get_ned13_tilecoords(lat, lon, step)
406    if step == 1:
407        return 'n%02dw%03d' % (y, -x)
408    return 'n%02.2fw%03.2f' % (y, -x)
409                       
410def prep_ned13_data_file(basename, env):
411    print 'Preparing NED 1/3" data file:', basename
412    print '  Converting to GeoTIFF...'
413    sourcepath = get_ned13_tilepath(basename)
414    geotiff = path.join(TEMPDIR, basename + '.tif')
415    if not path.isfile(geotiff):
416        cmd = 'gdal_translate "%s" "%s"' % (sourcepath, geotiff)
417        call(cmd, shell = True)
418
419    print '  Generating contour lines...'
420    # split the GeoTIFF, since it's often too large otherwise
421    for y in numpy.arange(env.miny, env.maxy, STEP):
422        for x in numpy.arange(env.minx, env.maxx, STEP):
423            print '  Cutting geotiff slice...'
424            nedslice = get_ned13_slice('ned', x, y)
425            if not path.isfile(nedslice):
426                cmd = 'gdalwarp -q -te %f %f %f %f "%s" "%s"' % \
427                    (x, y, x+STEP, y+STEP, geotiff, nedslice)
428                call(cmd, shell=True)
429           
430            print '  Generating contour lines...'
431            contourbasefile = path.join(TEMPDIR, 'contours_' + str(x) + '_' + str(y))
432            contourfile = contourbasefile + '.shp'
433            if not path.isfile(contourfile):
434                cmd = 'gdal_contour -i %f -snodata 32767 -a height "%s" "%s"' % \
435                    (CONTOUR_INTERVAL, nedslice, contourfile)
436                call(cmd, shell=True)
437
438                print '  Importing contour lines...'
439                # NOTE: this assumes that the table is already set up
440                cmd = 'shp2pgsql -a -g way "%s" "%s" | psql -q "%s"' % \
441                    (contourfile, CONTOURS_TABLE, DATABASE)
442                call(cmd, shell=True)
443               
444                # Clear contents (but keep file to prevent us from importing
445                # these contours again).
446                writeEmpty(contourfile)
447               
448                # Remove shape index and attribute files.
449                tryRemove(contourbasefile + ".shx")
450                tryRemove(contourbasefile + ".dbf")
451
452            print '  Generating hillshade slice...'
453            hillshadeslice = get_ned13_slice('hillshade', x, y)
454            hillshadeslicePng = get_ned13_slice('hillshade', x, y, '.png')
455            if not path.isfile(hillshadeslicePng):
456                cmd = '"%s" "%s" "%s" -z 0.00001' % \
457                    (HILLSHADE, nedslice, hillshadeslice)
458                call(cmd, shell = True)
459                # convert to PNG + world file to save space
460                cmd = 'gdal_translate -quiet -of PNG -co WORLDFILE=YES "%s" "%s"' % \
461                    (hillshadeslice, hillshadeslicePng)
462                call(cmd, shell = True)
463                tryRemove(hillshadeslice)
464
465            print '  Generating colormap slice...'
466            colormapslice = get_ned13_slice('colormap', x, y)
467            colormapslicePng = get_ned13_slice('colormap', x, y, '.png')
468            if not path.isfile(colormapslicePng):
469                cmd = '"%s" "%s" "%s" "%s"' % \
470                    (COLORRELIEF, nedslice, COLORFILE, colormapslice)
471                call(cmd, shell = True)
472                # convert to PNG + world file to save space
473                cmd = 'gdal_translate -quiet -of PNG -co WORLDFILE=YES "%s" "%s"' % \
474                    (colormapslice, colormapslicePng)
475                call(cmd, shell = True)
476                tryRemove(colormapslice)
477
478            writeEmpty(nedslice) # don't need the raw slice anymore.
479               
480    writeEmpty(geotiff) # done with this GeoTIFF.
481
482def pad_envelope(envLL, z, ntiles):
483    """Extends the specified area with the width/height of a subtile toward
484    increasing x/y, plus one unit of padding on all sides. Used to ensure that
485    enough data is included for rendering a set of tiles."""
486    (xmin, ymax) = GOOGLE_PROJECTION.fromLLtoPixel((envLL.minx, envLL.miny), z)
487    (xmax, ymin) = GOOGLE_PROJECTION.fromLLtoPixel((envLL.maxx, envLL.maxy), z)
488    xmin -= BORDER_WIDTH
489    ymin -= BORDER_WIDTH
490    xmax += BORDER_WIDTH + SUBTILE_SIZE
491    ymax += BORDER_WIDTH + SUBTILE_SIZE
492    (lonmin, latmin) = GOOGLE_PROJECTION.fromPixelToLL((xmin, ymax), z)
493    (lonmax, latmax) = GOOGLE_PROJECTION.fromPixelToLL((xmax, ymin), z)
494    return Envelope(lonmin, latmin, lonmax, latmax)
495
496# TODO: rename method... e.g. render_topo_tile
497def render_geotiff_tile(z, x, y, ntiles, mapname):
498    bbox = get_bbox_merc(z, x, y, ntiles, False)
499    bboxLL = get_bbox_ll(z, x, y, ntiles, False)
500    destdir = getTileDir(mapname, z)
501    # NOTE: gdalwarp won't save as png directly, hence this "hack"
502    destTilePath = getTilePath(mapname, z, x, y, '.tif')
503    finalTilePath = getTilePath(mapname, z, x, y, '.png')
504    if os.path.isfile(destTilePath):
505        pass
506    else:
507        ensureDirExists(destdir)
508        nedSlices = get_ned13_slices(mapname, bboxLL, '.png')
509        tilesize = getTileSize(ntiles, False)
510        if len(nedSlices) > 0:
511            cmd = 'gdalwarp -q -t_srs "%s" -te %f %f %f %f -ts %d %d -r lanczos -dstnodata 255 ' % \
512                (OSM_SRS, bbox.minx, bbox.miny, bbox.maxx, bbox.maxy, tilesize, tilesize)
513            for slice in nedSlices:
514                cmd = cmd + '"' + slice + '" '
515            cmd = cmd + '"' + destTilePath + '"'
516            print cmd
517            call(cmd, shell=True)
518            # convert to png and remove tif (to conserve space)
519            cmd = 'convert "%s" "%s" && rm "%s"' % (destTilePath, finalTilePath, destTilePath)
520            call(cmd, shell=True)
521        else:
522            return False
523    return True
524
525def slice_tile(z, x, y, ntiles, mapname, destSuffix = '.png'):
526    srcTilePath = getTilePath(mapname, z, x, y, '.png')
527    for dx in range(0, ntiles):
528        destDir = path.join(getTileDir(mapname, z), str(x+dx))
529        ensureDirExists(destDir)
530        for dy in range(0, ntiles):
531            destfile = path.join(destDir, str(y+dy) + destSuffix)
532            if not path.isfile(destfile):
533                offsetx = dx * SUBTILE_SIZE
534                offsety = dy * SUBTILE_SIZE
535                cmd = 'convert "%s" -crop %dx%d+%d+%d +repage "%s"' % \
536                    (srcTilePath, SUBTILE_SIZE, SUBTILE_SIZE, offsetx, offsety,
537                     path.join(destDir, str(y+dy) + destSuffix))
538                call(cmd, shell=True)
539            else:
540                pass
541    tryRemove(srcTilePath)
542
543
544def merge_subtiles(z, x, y, mapname, suffix = '.jpg'):
545    """Merges (up to) four subtiles from the next higher
546    zoom level into one subtile at the specified location"""
547#    print 'Merging tile ', z, x, y, mapname
548    cmd = 'convert -size 512x512 xc:white'
549    for dx in range(0,2):
550        for dy in range(0,2):
551            srcx = x*2 + dx
552            srcy = y*2 + dy
553            srcpath = getSubtilePath(mapname, z+1, srcx, srcy, suffix)
554            if os.path.isfile(srcpath):
555                cmd = cmd + ' "' + srcpath + '"'
556                cmd = cmd + ' -geometry +' + str(dx*256) + '+' + str(dy*256)
557                cmd = cmd + ' -composite'
558    cmd = cmd + ' -scale 256x256'
559    ensureDirExists(getSubtileDir(mapname, z, x))
560    destpath = getSubtilePath(mapname, z, x, y, suffix)
561    cmd = cmd + ' "' + destpath + '"'
562    call(cmd, shell=True)
563
564def create_jpeg_tile(z, x, y, quality):
565    print 'Creating jpeg tile at', z, x, y, '...',
566    colorreliefsrc = getSubtilePath('color-relief', z, x, y, '.jpg')
567    contourssrc = getSubtilePath('contours', z, x, y, '.png')
568    featuressrc = getSubtilePath('features', z, x, y, '.png')
569    desttile = getSubtilePath('jpeg' + str(quality), z, x, y, '.jpg')
570    if path.isfile(colorreliefsrc) and path.isfile(featuressrc):
571        ensureDirExists(path.dirname(desttile))
572
573        #tile = Image.open(colorreliefsrc)
574        #if path.isfile(contourssrc):
575        #    contours = Image.open(contourssrc)
576        #    tile.paste(contours, None)
577        #features = Image.open(featuressrc)
578        #tile.paste(features, None)
579        #ensure_dir_exists(path.dirname(desttile))
580        #tile.save(desttile, 'JPEG', quality=quality, optimize=True)
581       
582        # PIL generates internal errors opening the JPEG
583        # tiles so it's back to ImageMagick for now...
584        cmd = "convert " + colorreliefsrc;
585        if path.isfile(contourssrc):
586            cmd = cmd + " " + contourssrc + " -composite"
587        cmd = cmd + " " + featuressrc + " -composite"
588        cmd = cmd + " -quality " + str(quality) + " -strip " + desttile
589        call(cmd, shell=True)
590       
591        print 'done.'
592    else:
593        print 'source tiles not found. skipping.', colorreliefsrc, featuressrc
594
595
596# public methods
597
598def prepare_data_single(envLL, minz, maxz):
599    ntiles = NTILES[maxz]
600    envLL = pad_envelope(envLL, minz, ntiles)
601    tiles = get_ned13_tiles(envLL)
602    for tile in tiles:
603        print tile
604        prep_ned13_data_file(tile[0], tile[1])
605    print '  Converting meters to feet...'
606    cmd = 'echo "UPDATE %s SET height_ft = CAST(height * 3.28085 AS INT) WHERE height_ft IS NULL;" | psql -q "%s"' % (CONTOURS_TABLE, DATABASE)
607    call(cmd, shell=True)
608
609def prepare_data(envLL, minz, maxz):
610    ntiles = NTILES[maxz]
611    envLL = pad_envelope(envLL, minz, ntiles)
612    tiles = get_ned13_tiles(envLL)
613    queue = Queue(32)
614    renderers = {}
615    for i in range(NUM_THREADS):
616        renderer = RenderThread(queue, maxz, i)
617        render_thread = threading.Thread(target=renderer.render_loop)
618        render_thread.start()
619        renderers[i] = render_thread
620    for tile in tiles:
621        queue.put(('prepare_data', minz, tile[0], tile[1]))
622    for i in range(NUM_THREADS):
623        queue.put(None)
624    queue.join()
625    for i in range(NUM_THREADS):
626        renderers[i].join()
627    print 'Converting meters to feet...'
628    cmd = 'echo "UPDATE %s SET height_ft = CAST(height * 3.28085 AS INT) WHERE height_ft IS NULL;" | psql -q "%s"' % (CONTOURS_TABLE, DATABASE)
629    call(cmd, shell=True)
630    print 'Done.'
631
632def render_tiles(envLL, minz, maxz):
633    print 'Using mapnik version:', mapnik_version()
634    # Create mapnik rendering threads.
635    queue = Queue(32)
636    renderers = {}
637    for i in range(NUM_THREADS):
638        renderer = RenderThread(queue, maxz, i)
639        render_thread = threading.Thread(target=renderer.render_loop)
640        render_thread.start()
641        renderers[i] = render_thread
642    # Queue up jobs. High-to-low zoom, so we can merge topo-subtiles rather
643    # than render from scratch at lower zoom levels.
644    for z in range(maxz, minz-1, -1):
645        ntiles = NTILES[z]
646        (fromx, fromy) = get_tile_from_ll(envLL.maxy, envLL.minx, z)
647        (tox, toy) = get_tile_from_ll(envLL.miny, envLL.maxx, z)
648       
649        # Phil Gold: This causes all metatiles to be aligned to multiples of
650        # their zoom levels' NTILES value.  This lets me run renders with
651        # different windows into the same set of directories while reusing as
652        # much of the previous renders' work as possible.  It also lets me
653        # easily expire metatiles, since I always know which metatiles
654        # correspond to expired tiles.
655        for x in range(fromx - fromx%ntiles, tox+1, ntiles):
656            for y in range(fromy - fromy%ntiles, toy+1, ntiles):
657                queue.put(('render', z, x, y))
658        # Join threads after each completed level, because
659        # the color-relief layer of a lower zoom level depends
660        # on that of the next higher level being finished.
661        queue.join()
662    # Signal render threads to exit and join threads
663    for i in range(NUM_THREADS):
664        queue.put(None)
665    queue.join()
666    for i in range(NUM_THREADS):
667        renderers[i].join()
668
669def render_tiles_single(envLL, minz, maxz):
670    print 'Using mapnik version:', mapnik_version()
671    renderer = RenderThread(None, maxz, 1)
672    for z in range(maxz, minz-1, -1):
673        ntiles = NTILES[z]
674        (fromx, fromy) = get_tile_from_ll(envLL.maxy, envLL.minx, z)
675        (tox, toy) = get_tile_from_ll(envLL.miny, envLL.maxx, z)
676        for x in range(fromx, tox+1, ntiles):
677            for y in range(fromy, toy+1, ntiles):
678                renderer.render('render', z, x, y)
679
680def create_jpeg_tiles(envLL, minz, maxz, quality):
681    for z in range(minz, maxz+1):
682        (fromx, fromy) = get_tile_from_ll(envLL.maxy, envLL.minx, z)
683        (tox, toy) = get_tile_from_ll(envLL.miny, envLL.maxx, z)
684        for x in range(fromx, tox+1):
685            for y in range(fromy, toy+1):
686                create_jpeg_tile(z, x, y, quality)
687
Note: See TracBrowser for help on using the repository browser.