source: subversion/applications/rendering/mapnik/polytiles.py @ 27069

Last change on this file since 27069 was 27069, checked in by zverik, 8 years ago

polytiles.py. I hope it is allowed.

  • Property svn:executable set to *
File size: 20.3 KB
Line 
1#!/usr/bin/python
2# -*- coding: utf-8 -*-
3# This file is based on generate_tiles_multiprocess.py
4from math import pi,cos,sin,log,exp,atan
5from subprocess import call
6import sys, os
7from Queue import Queue
8import mapnik2 as mapnik
9import multiprocessing
10import psycopg2
11import shapely
12import ogr
13import sqlite3
14import getpass
15import argparse
16
17DEG_TO_RAD = pi/180
18RAD_TO_DEG = 180/pi
19TILE_SIZE = 256
20
21
22def box(x1,y1,x2,y2):
23    from shapely.geometry import Polygon
24    return Polygon([(x1,y1), (x2,y1), (x2,y2), (x1,y2)])
25
26def minmax (a,b,c):
27    a = max(a,b)
28    a = min(a,c)
29    return a
30
31class GoogleProjection:
32    def __init__(self,levels=18):
33        self.Bc = []
34        self.Cc = []
35        self.zc = []
36        self.Ac = []
37        c = 256
38        for d in range(0,levels):
39            e = c/2;
40            self.Bc.append(c/360.0)
41            self.Cc.append(c/(2 * pi))
42            self.zc.append((e,e))
43            self.Ac.append(c)
44            c *= 2
45               
46    def fromLLtoPixel(self,ll,zoom):
47         d = self.zc[zoom]
48         e = round(d[0] + ll[0] * self.Bc[zoom])
49         f = minmax(sin(DEG_TO_RAD * ll[1]),-0.9999,0.9999)
50         g = round(d[1] + 0.5*log((1+f)/(1-f))*-self.Cc[zoom])
51         return (e,g)
52     
53    def fromPixelToLL(self,px,zoom):
54         e = self.zc[zoom]
55         f = (px[0] - e[0])/self.Bc[zoom]
56         g = (px[1] - e[1])/-self.Cc[zoom]
57         h = RAD_TO_DEG * ( 2 * atan(exp(g)) - 0.5 * pi)
58         return (f,h)
59
60
61class ListWriter:
62    def __init__(self, f):
63        self.f = f
64
65    def __str__(self):
66        return "ListWriter({0})".format(self.f.name)
67
68    def write_poly(self, poly):
69        self.f.write("BBox: {0}\n".format(poly.bounds))
70
71    def write(self, x, y, z):
72        self.f.write("{0}/{1}/{2}\n".format(z,x,y))
73
74    def exists(self, x, y, z):
75        return False
76
77    def need_image(self):
78        return False
79
80    def multithreading(self):
81        return False
82
83    def close(self):
84        self.f.close()
85
86class FileWriter:
87    def __init__(self, tile_dir):
88        self.tile_dir = tile_dir
89        if not self.tile_dir.endswith('/'):
90            self.tile_dir = self.tile_dir + '/'
91        if not os.path.isdir(self.tile_dir):
92            os.mkdir(self.tile_dir)
93
94    def __str__(self):
95        return "FileWriter({0})".format(self.tile_dir)
96
97    def write_poly(self, poly):
98        pass
99
100    def tile_uri(self, x, y, z):
101        return '{0}{1}/{2}/{3}.png'.format(self.tile_dir, z, x, y)
102
103    def exists(self, x, y, z):
104        return os.path.isfile(self.tile_uri(x, y, z))
105
106    def write(self, x, y, z, image):
107        uri = self.tile_uri(x, y, z)
108        try:
109            os.makedirs(os.path.dirname(uri))
110        except OSError:
111            pass
112        image.save(uri, 'png256')
113
114    def need_image(self):
115        return True
116
117    def multithreading(self):
118        return True
119
120    def close(self):
121        pass
122
123class TMSWriter(FileWriter):
124    def tile_uri(self, x, y, z):
125        return '{0}{1}/{2}/{3}.png'.format(self.tile_dir, z, x, 2**z-1-y)
126
127    def __str__(self):
128        return "TMSWriter({0})".format(self.tile_dir)
129
130# https://github.com/mapbox/mbutil/blob/master/mbutil/util.py
131class MBTilesWriter:
132    def __init__(self, setname, filename, overlay=False, version=1, description=None):
133        self.filename = filename
134        if not self.filename.endswith('.mbtiles'):
135            self.filename = self.filename + '.mbtiles'
136        self.con = sqlite3.connect(self.filename)
137        self.cur = self.con.cursor()
138        self.cur.execute("""PRAGMA synchronous=0""")
139        self.cur.execute("""PRAGMA locking_mode=EXCLUSIVE""")
140        #self.cur.execute("""PRAGMA journal_mode=TRUNCATE""")
141        self.cur.execute("""create table if not exists tiles (zoom_level integer, tile_column integer, tile_row integer, tile_data blob);""")
142        self.cur.execute("""create table if not exists metadata (name text, value text);""")
143        self.cur.execute("""create unique index if not exists name on metadata (name);""")
144        self.cur.execute("""create unique index if not exists tile_index on tiles (zoom_level, tile_column, tile_row);""")
145        metadata = [ ('name', setname), ('format', 'png'), ('type', 'overlay' if overlay else 'baselayer'), ('version', version) ]
146        if description:
147            metadata.append(('description', description))
148        for name, value in metadata:
149            self.cur.execute('insert or replace into metadata (name, value) values (?, ?)', (name, value))
150
151    def __str__(self):
152        return "MBTilesWriter({0})".format(self.filename)
153
154    def write_poly(self, poly):
155        bbox = poly.bounds
156        self.cur.execute("""select value from metadata where name='bounds'""")
157        result = self.cur.fetchone
158        if result:
159            b = result['value'].split(',')
160            oldbbox = box(int(b[0]), int(b[1]), int(b[2]), int(b[3]))
161            bbox = bbox.union(oldbbox).bounds
162        self.cur.execute("""insert or replace into metadata (name, value) values ('bounds', ?)""", ','.join(bbox))
163
164    def exists(self, x, y, z):
165        self.cur.execute("""select 1 from tiles where zoom_level = ? and tile_column = ? and tile_row = ?""", (z, x, 2**z-1-y))
166        return self.cur.fetchone()
167
168    def write(self, x, y, z, image):
169        self.cur.execute("""insert or replace into tiles (zoom_level, tile_column, tile_row, tile_data) values (?, ?, ?, ?);""", (z, x, 2**z-1-y, sqlite3.Binary(image.tostring('png256'))))
170
171    def need_image(self):
172        return True
173
174    def multithreading(self):
175        return False
176
177    def close(self):
178        self.cur.execute("""ANALYZE;""")
179        self.cur.execute("""VACUUM;""")
180        self.cur.close()
181        self.con.close()
182
183
184# todo: make queue-based writer
185class ThreadedWriter:
186    def __init__(self, writer):
187        self.writer = writer
188        self.queue = Queue(10)
189
190    def __str__(self):
191        return "Threaded{0}".format(self.writer)
192
193    def write_poly(self, poly):
194        pass
195
196    def exists(self, x, y, z):
197        pass
198
199    def write(self, x, y, z, image):
200        pass
201
202    def need_image(self):
203        return writer.need_image()
204
205    def multithreading(self):
206        return True
207
208    def close(self):
209        writer.close()
210
211
212class RenderThread:
213    def __init__(self, writer, mapfile, q, printLock, verbose=True):
214        self.writer = writer
215        self.q = q
216        self.mapfile = mapfile
217        self.printLock = printLock
218        self.verbose = verbose
219
220    def render_tile(self, x, y, z):
221        # Calculate pixel positions of bottom-left & top-right
222        p0 = (x * TILE_SIZE, (y + 1) * TILE_SIZE)
223        p1 = ((x + 1) * TILE_SIZE, y * TILE_SIZE)
224
225        # Convert to LatLong (EPSG:4326)
226        l0 = self.tileproj.fromPixelToLL(p0, z);
227        l1 = self.tileproj.fromPixelToLL(p1, z);
228
229        # Convert to map projection (e.g. mercator co-ords EPSG:900913)
230        c0 = self.prj.forward(mapnik.Coord(l0[0],l0[1]))
231        c1 = self.prj.forward(mapnik.Coord(l1[0],l1[1]))
232
233        # Bounding box for the tile
234        if hasattr(mapnik,'mapnik_version') and mapnik.mapnik_version() >= 800:
235            bbox = mapnik.Box2d(c0.x,c0.y, c1.x,c1.y)
236        else:
237            bbox = mapnik.Envelope(c0.x,c0.y, c1.x,c1.y)
238        render_size = TILE_SIZE
239        self.m.resize(render_size, render_size)
240        self.m.zoom_to_box(bbox)
241        self.m.buffer_size = 128
242
243        # Render image with default Agg renderer
244        im = mapnik.Image(render_size, render_size)
245        mapnik.render(self.m, im)
246        self.writer.write(x, y, z, im)
247
248
249    def loop(self):
250        self.m = mapnik.Map(TILE_SIZE, TILE_SIZE)
251        # Load style XML
252        mapnik.load_map(self.m, self.mapfile, True)
253        # Obtain <Map> projection
254        self.prj = mapnik.Projection(self.m.srs)
255        # Projects between tile pixel co-ordinates and LatLong (EPSG:4326)
256        self.tileproj = GoogleProjection()
257
258        while True:
259            #Fetch a tile from the queue and render it
260            r = self.q.get()
261            if (r == None):
262                self.q.task_done()
263                break
264            else:
265                (x, y, z) = r
266
267            exists= ""
268            if self.writer.exists(x, y, z):
269                exists= "exists"
270            elif self.writer.need_image():
271                self.render_tile(x, y, z)
272            else:
273                self.writer.write(x, y, z)
274            empty = ''
275            #if os.path.exists(tile_uri):
276            #    bytes=os.stat(tile_uri)[6]
277            #    empty= ''
278            #    if bytes == 103:
279            #        empty = " Empty Tile "
280            #else:
281            #    empty = " Missing "
282            if self.verbose:
283                self.printLock.acquire()
284                print z, x, y, exists, empty
285                self.printLock.release()
286            self.q.task_done()
287
288class ListGenerator:
289    def __init__(self, f):
290        self.f = f
291
292    def __str__(self):
293        return "ListGenerator({0})".format(self.f.name)
294
295    def generate(self, queue):
296        import re
297        for line in self.f:
298            m = re.search(r"(\d{1,2})\D+(\d+)\D+(\d+)", line)
299            if m:
300                queue.put((int(m.group(2)), int(m.group(3)), int(m.group(1))))
301
302
303class PolyGenerator:
304    def __init__(self, poly, zooms):
305        self.poly = poly
306        self.zooms = zooms
307        self.zooms.sort()
308
309    def __str__(self):
310        return "PolyGenerator({0}, {1})".format(self.poly.bounds, self.zooms)
311
312    def generate(self, queue):
313        gprj = GoogleProjection(self.zooms[-1]+1) 
314
315        bbox = self.poly.bounds
316        ll0 = (bbox[0], bbox[3])
317        ll1 = (bbox[2], bbox[1])
318
319        for z in self.zooms:
320            px0 = gprj.fromLLtoPixel(ll0, z)
321            px1 = gprj.fromLLtoPixel(ll1, z)
322
323            for x in range(int(px0[0]/float(TILE_SIZE)), int(px1[0]/float(TILE_SIZE))+1):
324                # Validate x co-ordinate
325                if (x < 0) or (x >= 2**z):
326                    continue
327                for y in range(int(px0[1]/float(TILE_SIZE)), int(px1[1]/float(TILE_SIZE))+1):
328                    # Validate x co-ordinate
329                    if (y < 0) or (y >= 2**z):
330                        continue
331
332                    # Calculate pixel positions of bottom-left & top-right
333                    tt_p0 = (x * TILE_SIZE, (y + 1) * TILE_SIZE)
334                    tt_p1 = ((x + 1) * TILE_SIZE, y * TILE_SIZE)
335
336                    # Convert to LatLong (EPSG:4326)
337                    tt_l0 = gprj.fromPixelToLL(tt_p0, z);
338                    tt_l1 = gprj.fromPixelToLL(tt_p1, z);
339
340                    tt_p = box(tt_l0[0], tt_l1[1], tt_l1[0], tt_l0[1])
341                    if not self.poly.intersects(tt_p):
342                        continue
343
344                    # Submit tile to be rendered into the queue
345                    t = (x, y, z)
346                    queue.put(t)
347
348
349def render_tiles(generator, mapfile, writer, num_threads=1, verbose=True):
350    if verbose:
351        print "render_tiles(",generator, mapfile, writer, ")"
352
353    # Launch rendering threads
354    queue = multiprocessing.JoinableQueue(32 if writer.multithreading() else 0)
355    printLock = multiprocessing.Lock()
356    renderers = {}
357    if writer.multithreading():
358        for i in range(num_threads):
359            renderer = RenderThread(writer, mapfile, queue, printLock, verbose=verbose)
360            render_thread = multiprocessing.Process(target=renderer.loop)
361            render_thread.start()
362            #print "Started render thread %s" % render_thread.getName()
363            renderers[i] = render_thread
364
365    generator.generate(queue)
366
367    if writer.multithreading():
368        # Signal render threads to exit by sending empty request to queue
369        for i in range(num_threads):
370            queue.put(None)
371        # wait for pending rendering jobs to complete
372        queue.join()
373        for i in range(num_threads):
374            renderers[i].join()
375    else:
376        renderer = RenderThread(writer, mapfile, queue, printLock, verbose=verbose)
377        queue.put(None)
378        renderer.loop()
379
380
381def poly_parse(fp):
382    poly = []
383    data = False
384    for l in fp:
385        l = l.strip()
386        if not l: continue
387        if l == 'END': break
388        if l == '1':
389            data = True
390            continue
391        if not data: continue
392        poly.append(map(lambda x: float(x.strip()), l.split()[:2]))
393    return poly
394
395
396def read_poly(filename):
397    from shapely.geometry import Polygon
398    poly = Polygon(poly_parse(open(filename)))
399    return poly
400
401def project(geom, from_epsg=900913, to_epsg=4326):
402    # source: http://hackmap.blogspot.com/2008/03/ogr-python-projection.html
403    to_srs = ogr.osr.SpatialReference()
404    to_srs.ImportFromEPSG(to_epsg)
405
406    from_srs = ogr.osr.SpatialReference()
407    from_srs.ImportFromEPSG(from_epsg)
408
409    ogr_geom = ogr.CreateGeometryFromWkb(geom.wkb)
410    ogr_geom.AssignSpatialReference(from_srs)
411
412    ogr_geom.TransformTo(to_srs)
413    from shapely.wkb import loads
414    return loads(ogr_geom.ExportToWkb())
415
416def read_db(db, osm_id=0):
417    # Zero for DB bbox
418    cur = db.cursor()
419    if osm_id:
420        cur.execute("""SELECT way FROM planet_osm_polygon WHERE osm_id = %s;""", (osm_id,))
421    else:
422        cur.execute("""SELECT ST_ConvexHull(ST_Collect(way)) FROM planet_osm_polygon;""")
423    way = cur.fetchone()[0]
424    cur.close()
425    from shapely.wkb import loads
426    poly = loads(way.decode('hex'))
427    return project(poly)
428
429def read_cities(db, osm_id=0):
430    from shapely.wkb import loads
431    cur = db.cursor()
432    if osm_id:
433        cur.execute("""SELECT ST_Union(pl.way) FROM planet_osm_polygon pl, planet_osm_polygon b WHERE b.osm_id = %s AND pl.place IN ('town', 'city') AND AT_Area(pl.way) < 500*1000*1000 AND ST_Contains(b.way, pl.way);""", (osm_id,))
434    else:
435        cur.execute("""SELECT ST_Union(way) FROM planet_osm_polygon WHERE place IN ('town', 'city') AND ST_Area(way) < 500*1000*1000;""")
436    result = cur.fetchone()
437    way = result[0] if result else shapely.geometry.Polygon()
438    poly = loads(way.decode('hex'))
439    if osm_id:
440        cur.execute("""SELECT ST_Union(ST_Buffer(p.way, 5000)) FROM planet_osm_point p, planet_osm_polygon b WHERE b.osm_id=%s AND ST_Contains(b.way, p.way) AND p.place IN ('town', 'city') AND NOT EXISTS(SELECT 1 FROM planet_osm_polygon pp WHERE pp.name=p.name AND ST_Contains(pp.way, p.way));""", (osm_id,))
441    else:
442        cur.execute("""SELECT ST_Union(ST_Buffer(p.way, 5000)) FROM planet_osm_point p WHERE p.place in ('town', 'city') AND NOT EXISTS(SELECT 1 FROM planet_osm_polygon pp WHERE pp.name=p.name AND ST_Contains(pp.way, p.way));""")
443    result = cur.fetchone()
444    if result:
445        poly = poly.union(loads(result[0].decode('hex')))
446    return project(poly)
447
448
449if __name__ == "__main__":
450    try:
451        mapfile = os.environ['MAPNIK_MAP_FILE']
452    except KeyError:
453        mapfile = os.getcwd() + '/osm.xml'
454
455    default_user = getpass.getuser()
456
457    parser = argparse.ArgumentParser(description='Generate mapnik tiles for OSM polygon')
458    apg_input = parser.add_argument_group('Input')
459    apg_input.add_argument("-b", "--bbox", nargs=4, type=float, metavar=('X1', 'Y1', 'X2', 'Y2'), help="generate tiles inside a bounding box")
460    apg_input.add_argument('-p', '--poly', type=argparse.FileType('r'), help='use a poly file for area')
461    apg_input.add_argument("-a", "--area", type=int, metavar='OSM_ID', help="generate tiles inside an OSM polygon: positive for polygons, negative for relations, 0 for whole database")
462    apg_input.add_argument("-c", "--cities", type=int, metavar='OSM_ID', help='generate tiles for all towns inside a polygon')
463    apg_input.add_argument('-l', '--list', type=argparse.FileType('r'), metavar='TILES.LST', help='process tile list')
464    apg_output = parser.add_argument_group('Output')
465    apg_output.add_argument('-t', '--tiledir', metavar='DIR', help='output tiles to directory (default: {0}/tiles)'.format(os.getcwd()))
466    apg_output.add_argument('--tms', action='store_true', help='write files in TMS order', default=False)
467    apg_output.add_argument('-m', '--mbtiles', help='generate mbtiles file')
468    apg_output.add_argument('-x', '--export', type=argparse.FileType('w'), metavar='TILES.LST', help='save tile list into file')
469    apg_output.add_argument('-z', '--zooms', type=int, nargs=2, metavar=('ZMIN', 'ZMAX'), help='range of zoom levels to render (default: 6 12)', default=(6, 12))
470    apg_other = parser.add_argument_group('Settings')
471    apg_other.add_argument('-s', '--style', help='style file for mapnik (default: {0})'.format(mapfile), default=mapfile)
472    apg_other.add_argument('--threads', type=int, metavar='N', help='number of threads (default: 2)', default=2)
473    apg_other.add_argument('-q', '--quiet', dest='verbose', action='store_false', help='do not print any information',  default=True)
474    apg_db = parser.add_argument_group('Database (for poly/cities)')
475    apg_db.add_argument('-d', '--dbname', metavar='DB', help='database (default: gis)', default='gis')
476    apg_db.add_argument('--host', help='database host', default='localhost')
477    apg_db.add_argument('--port', type=int, help='database port', default='5432')
478    apg_db.add_argument('-u', '--user', help='user name for db (default: {0})'.format(default_user), default=default_user)
479    apg_db.add_argument('-w', '--password', action='store_true', help='ask for password', default=False)
480    options = parser.parse_args()
481
482    # check for required argument
483    if options.bbox == None and options.poly == None and options.cities == None and options.list == None and options.area == None:
484        parser.print_help()
485        sys.exit()
486
487    # writer
488    if options.tiledir:
489        writer = FileWriter(options.tiledir) if not options.tms else TMSWriter(options.tiledir)
490    elif options.mbtiles:
491        writer = MBTilesWriter(options.mbtiles)
492    elif options.export:
493        writer = ListWriter(options.export)
494    else:
495        writer = FileWriter(os.getcwd() + '/tiles') if not options.tms else TMSWriter(os.getcwd() + '/tiles')
496
497    # input and process
498    poly = None
499    if options.bbox:
500        b = options.bbox
501        tpoly = box(b[0], b[1], b[2], b[3])
502        poly = tpoly if not poly else poly.intersection(tpoly)
503    if options.poly:
504        tpoly = poly_parse(options.poly)
505        poly = tpoly if not poly else poly.intersection(tpoly)
506    db = None
507    if options.area != None or options.cities != None:
508        passwd = ""
509        if options.password:
510            passwd = getpass.getpass("Please enter your password: ")
511
512        try:
513            db = psycopg2.connect(database=options.dbname, user=options.user, password=passwd, host=options.host, port=options.port)
514        except Exception, e:
515            print "Error connecting to database: ", e.pgerror
516            sys.exit(1)
517    if options.area != None:
518        tpoly = read_db(db, options.area)
519        poly = tpoly if not poly else poly.intersection(tpoly)
520    if options.cities != None:
521        tpoly = read_cities(db, options.poly)
522        poly = tpoly if not poly else poly.intersection(tpoly)
523
524    if db:
525        db.close()
526
527    if options.list:
528        generator = ListGenerator(options.list)
529    elif poly:
530        generator = PolyGenerator(poly, range(options.zooms[0], options.zooms[1] + 1))
531    else:
532        print "Please specify a region for rendering."
533        sys.exit()
534
535    render_tiles(generator, options.style, writer, num_threads=options.threads, verbose=options.verbose)
536    writer.close()
537    sys.exit()
538
539    #vas = read_db(db, -1114252)
540    #vyb = read_db(db, -1138534)
541    writer = MBTilesWriter('Vasileostrovsky', 'vo')
542    #render_tiles(vas, mapfile, FileWriter(tile_dir), 12, 13, "spb")
543    render_tiles(vas, mapfile, writer, [10,11,12, 13,14], "spb")
544    render_tiles(vyb, mapfile, writer, [10,11,12,13], "vyb")
545    writer.close()
546    sys.exit()
547    #-------------------------------------------------------------------------
548    #
549    # Change the following for different bounding boxes and zoom levels
550    #
551    # Start with an overview
552    # World
553
554    ######
555   
556    # Обзорная часть, по bbox
557    bbox = box(19.364370,41.116250,180.000000,77.989470)
558    render_tiles(bbox, mapfile, tile_dir, [2, 6], "Россия 2-6")
559
560    #sys.exit()
561
562    # Калининградская область по границам
563    # bbox = (19.364370,54.307900,22.916120,55.412780)
564    poly = Polygon(poly_parse(open('/home/maks/OSM/russia-kaliningrad-mapnik.poly' )))
565    render_tiles(poly, mapfile, tile_dir, [7, 12], "Калиниград 7-13")
566
567    # Материковая часть страны по границе
568    bbox = box(26.574820,41.116250,100.000000,77.989470)
569    poly = Polygon(poly_parse(open('/home/maks/OSM/russia-mapnik.poly')))
570    render_tiles(poly.intersection(bbox), mapfile, tile_dir, [7, 10], "Россия 7-13")
571
Note: See TracBrowser for help on using the repository browser.