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

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

bugfixes in polytiles.py

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