From 660adf9123920e9dba297d1ef1dd3c2f163a1db5 Mon Sep 17 00:00:00 2001 From: jomu Date: Fri, 15 Jul 2016 09:40:01 +0200 Subject: [PATCH] added renderpoly.py --- scripts/renderpoly.py | 549 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 549 insertions(+) create mode 100755 scripts/renderpoly.py diff --git a/scripts/renderpoly.py b/scripts/renderpoly.py new file mode 100755 index 0000000..4f8e6f8 --- /dev/null +++ b/scripts/renderpoly.py @@ -0,0 +1,549 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- + +# this file is based on generate_tiles_multiprocess.py +# run it without arguments to see options list + +from math import pi,cos,sin,log,exp,atan +from subprocess import call +import sys, os +#import mapnik2 as mapnik +import multiprocessing +import psycopg2 +from shapely.geometry import Polygon +from shapely.wkb import loads +import ogr +import sqlite3 +import getpass +import argparse +import mapnik + +#try: +# import mapnik2 as mapnik +#except: +# import mapnik + + +DEG_TO_RAD = pi/180 +RAD_TO_DEG = 180/pi +TILE_SIZE = 256 + + +def box(x1,y1,x2,y2): + return Polygon([(x1,y1), (x2,y1), (x2,y2), (x1,y2)]) + +def minmax (a,b,c): + a = max(a,b) + a = min(a,c) + return a + +class GoogleProjection: + def __init__(self,levels=19): + self.Bc = [] + self.Cc = [] + self.zc = [] + self.Ac = [] + c = 256 + for d in range(0,levels): + e = c/2; + self.Bc.append(c/360.0) + self.Cc.append(c/(2 * pi)) + self.zc.append((e,e)) + self.Ac.append(c) + c *= 2 + + def fromLLtoPixel(self,ll,zoom): + d = self.zc[zoom] + e = round(d[0] + ll[0] * self.Bc[zoom]) + f = minmax(sin(DEG_TO_RAD * ll[1]),-0.9999,0.9999) + g = round(d[1] + 0.5*log((1+f)/(1-f))*-self.Cc[zoom]) + return (e,g) + + def fromPixelToLL(self,px,zoom): + e = self.zc[zoom] + f = (px[0] - e[0])/self.Bc[zoom] + g = (px[1] - e[1])/-self.Cc[zoom] + h = RAD_TO_DEG * ( 2 * atan(exp(g)) - 0.5 * pi) + return (f,h) + + +class ListWriter: + def __init__(self, f): + self.f = f + self.force = False + + def __str__(self): + return "ListWriter({0})".format(self.f.name) + + def write_poly(self, poly): + self.f.write("BBox: {0}\n".format(poly.bounds)) + + def write(self, x, y, z): + self.f.write("{0}/{1}/{2}\n".format(z,x,y)) + + def exists(self, x, y, z): + return False + + def need_image(self): + return False + + def multithreading(self): + return False + + def close(self): + self.f.close() + +class FileWriter: + def __init__(self, tile_dir): + self.tile_dir = tile_dir + self.force = false + if not self.tile_dir.endswith('/'): + self.tile_dir = self.tile_dir + '/' + if not os.path.isdir(self.tile_dir): + os.mkdir(self.tile_dir) + + def __init__(self, tile_dir, force): + self.tile_dir = tile_dir + self.force = force + if not self.tile_dir.endswith('/'): + self.tile_dir = self.tile_dir + '/' + if not os.path.isdir(self.tile_dir): + os.mkdir(self.tile_dir) + + def __str__(self): + return "FileWriter({0})".format(self.tile_dir) + + def write_poly(self, poly): + pass + + def tile_uri(self, x, y, z): + return '{0}{1}/{2}/{3}.png'.format(self.tile_dir, z, x, y) + + def exists(self, x, y, z): + return os.path.isfile(self.tile_uri(x, y, z)) + + def write(self, x, y, z, image): + uri = self.tile_uri(x, y, z) + try: + os.makedirs(os.path.dirname(uri)) + except OSError: + pass + image.save(uri, 'png256') + + def need_image(self): + return True + + def multithreading(self): + return True + + def close(self): + pass + +class TMSWriter(FileWriter): + def tile_uri(self, x, y, z): + return '{0}{1}/{2}/{3}.png'.format(self.tile_dir, z, x, 2**z-1-y) + + def __str__(self): + return "TMSWriter({0})".format(self.tile_dir) + +# https://github.com/mapbox/mbutil/blob/master/mbutil/util.py +class MBTilesWriter: + def __init__(self, filename, setname, overlay=False, version=1, description=None): + self.filename = filename + if not self.filename.endswith('.mbtiles'): + self.filename = self.filename + '.mbtiles' + self.con = sqlite3.connect(self.filename) + self.cur = self.con.cursor() + self.cur.execute("""PRAGMA synchronous=0""") + self.cur.execute("""PRAGMA locking_mode=EXCLUSIVE""") + #self.cur.execute("""PRAGMA journal_mode=TRUNCATE""") + self.cur.execute("""create table if not exists tiles (zoom_level integer, tile_column integer, tile_row integer, tile_data blob);""") + self.cur.execute("""create table if not exists metadata (name text, value text);""") + self.cur.execute("""create unique index if not exists name on metadata (name);""") + self.cur.execute("""create unique index if not exists tile_index on tiles (zoom_level, tile_column, tile_row);""") + metadata = [ ('name', setname), ('format', 'png'), ('type', 'overlay' if overlay else 'baselayer'), ('version', version) ] + if description: + metadata.append(('description', description)) + for name, value in metadata: + self.cur.execute('insert or replace into metadata (name, value) values (?, ?)', (name, value)) + self.force = False + + def __str__(self): + return "MBTilesWriter({0})".format(self.filename) + + def write_poly(self, poly): + bbox = poly.bounds + self.cur.execute("""select value from metadata where name='bounds'""") + result = self.cur.fetchone + if result: + b = result['value'].split(',') + oldbbox = box(int(b[0]), int(b[1]), int(b[2]), int(b[3])) + bbox = bbox.union(oldbbox).bounds + self.cur.execute("""insert or replace into metadata (name, value) values ('bounds', ?)""", ','.join(bbox)) + + def exists(self, x, y, z): + self.cur.execute("""select 1 from tiles where zoom_level = ? and tile_column = ? and tile_row = ?""", (z, x, 2**z-1-y)) + return self.cur.fetchone() + + def write(self, x, y, z, image): + 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')))) + + def need_image(self): + return True + + def multithreading(self): + return False + + def close(self): + self.cur.execute("""ANALYZE;""") + self.cur.execute("""VACUUM;""") + self.cur.close() + self.con.close() + + +# todo: make queue-based writer +class ThreadedWriter: + def __init__(self, writer): + self.writer = writer + self.queue = multiprocessing.Queue(10) + + def __str__(self): + return "Threaded{0}".format(self.writer) + + def write_poly(self, poly): + pass + + def exists(self, x, y, z): + pass + + def write(self, x, y, z, image): + pass + + def need_image(self): + return writer.need_image() + + def multithreading(self): + return True + + def close(self): + writer.close() + + +class RenderThread: + def __init__(self, writer, mapfile, q, printLock, verbose=True): + self.writer = writer + self.q = q + self.mapfile = mapfile + self.printLock = printLock + self.verbose = verbose + + def render_tile(self, x, y, z): + # Calculate pixel positions of bottom-left & top-right + p0 = (x * TILE_SIZE, (y + 1) * TILE_SIZE) + p1 = ((x + 1) * TILE_SIZE, y * TILE_SIZE) + + # Convert to LatLong (EPSG:4326) + l0 = self.tileproj.fromPixelToLL(p0, z); + l1 = self.tileproj.fromPixelToLL(p1, z); + + # Convert to map projection (e.g. mercator co-ords EPSG:900913) + c0 = self.prj.forward(mapnik.Coord(l0[0],l0[1])) + c1 = self.prj.forward(mapnik.Coord(l1[0],l1[1])) + + # Bounding box for the tile + if hasattr(mapnik,'mapnik_version') and mapnik.mapnik_version() >= 800: + bbox = mapnik.Box2d(c0.x,c0.y, c1.x,c1.y) + else: + bbox = mapnik.Envelope(c0.x,c0.y, c1.x,c1.y) + render_size = TILE_SIZE + self.m.resize(render_size, render_size) + self.m.zoom_to_box(bbox) + self.m.buffer_size = 128 + + # Render image with default Agg renderer + im = mapnik.Image(render_size, render_size) + mapnik.render(self.m, im) + self.writer.write(x, y, z, im) + + + def loop(self): + self.m = mapnik.Map(TILE_SIZE, TILE_SIZE) + # Load style XML + mapnik.load_map(self.m, self.mapfile, True) + # Obtain projection + self.prj = mapnik.Projection(self.m.srs) + # Projects between tile pixel co-ordinates and LatLong (EPSG:4326) + self.tileproj = GoogleProjection() + + while True: + #Fetch a tile from the queue and render it + r = self.q.get() + if (r == None): + self.q.task_done() + break + else: + (x, y, z) = r + + exists= "" + if ((self.writer.exists(x, y, z)) & (self.writer.force == False)): + exists= "exists" + elif self.writer.need_image(): + self.render_tile(x, y, z) + else: + self.writer.write(x, y, z) + empty = '' + #if os.path.exists(tile_uri): + # bytes=os.stat(tile_uri)[6] + # empty= '' + # if bytes == 103: + # empty = " Empty Tile " + #else: + # empty = " Missing " + if self.verbose: + self.printLock.acquire() + print z, x, y, exists, empty + self.printLock.release() + self.q.task_done() + +class ListGenerator: + def __init__(self, f): + self.f = f + + def __str__(self): + return "ListGenerator({0})".format(self.f.name) + + def generate(self, queue): + import re + for line in self.f: + m = re.search(r"(\d{1,2})\D+(\d+)\D+(\d+)", line) + if m: + queue.put((int(m.group(2)), int(m.group(3)), int(m.group(1)))) + + +class PolyGenerator: + def __init__(self, poly, zooms): + self.poly = poly + self.zooms = zooms + self.zooms.sort() + + def __str__(self): + return "PolyGenerator({0}, {1})".format(self.poly.bounds, self.zooms) + + def generate(self, queue): + gprj = GoogleProjection(self.zooms[-1]+1) + + bbox = self.poly.bounds + ll0 = (bbox[0], bbox[3]) + ll1 = (bbox[2], bbox[1]) + + for z in self.zooms: + px0 = gprj.fromLLtoPixel(ll0, z) + px1 = gprj.fromLLtoPixel(ll1, z) + + for x in range(int(px0[0]/float(TILE_SIZE)), int(px1[0]/float(TILE_SIZE))+1): + # Validate x co-ordinate + if (x < 0) or (x >= 2**z): + continue + for y in range(int(px0[1]/float(TILE_SIZE)), int(px1[1]/float(TILE_SIZE))+1): + # Validate x co-ordinate + if (y < 0) or (y >= 2**z): + continue + + # Calculate pixel positions of bottom-left & top-right + tt_p0 = (x * TILE_SIZE, (y + 1) * TILE_SIZE) + tt_p1 = ((x + 1) * TILE_SIZE, y * TILE_SIZE) + + # Convert to LatLong (EPSG:4326) + tt_l0 = gprj.fromPixelToLL(tt_p0, z); + tt_l1 = gprj.fromPixelToLL(tt_p1, z); + + tt_p = box(tt_l0[0], tt_l1[1], tt_l1[0], tt_l0[1]) + if not self.poly.intersects(tt_p): + continue + + # Submit tile to be rendered into the queue + t = (x, y, z) + queue.put(t) + + +def render_tiles(generator, mapfile, writer, num_threads=1, verbose=True): + if verbose: + print "render_tiles(",generator, mapfile, writer, ")" + + # Launch rendering threads + queue = multiprocessing.JoinableQueue(32 if writer.multithreading() else 0) + printLock = multiprocessing.Lock() + renderers = {} + if writer.multithreading(): + for i in range(num_threads): + renderer = RenderThread(writer, mapfile, queue, printLock, verbose=verbose) + render_thread = multiprocessing.Process(target=renderer.loop) + render_thread.start() + #print "Started render thread %s" % render_thread.getName() + renderers[i] = render_thread + + generator.generate(queue) + + if writer.multithreading(): + # Signal render threads to exit by sending empty request to queue + for i in range(num_threads): + queue.put(None) + # wait for pending rendering jobs to complete + queue.join() + for i in range(num_threads): + renderers[i].join() + else: + renderer = RenderThread(writer, mapfile, queue, printLock, verbose=verbose) + queue.put(None) + renderer.loop() + + +def poly_parse(fp): + poly = [] + data = False + for l in fp: + l = l.strip() + if not l: continue + if l == 'END': break + if l == '1': + data = True + continue + if not data: continue + poly.append(map(lambda x: float(x.strip()), l.split()[:2])) + return poly + + +def project(geom, from_epsg=900913, to_epsg=4326): + # source: http://hackmap.blogspot.com/2008/03/ogr-python-projection.html + to_srs = ogr.osr.SpatialReference() + to_srs.ImportFromEPSG(to_epsg) + + from_srs = ogr.osr.SpatialReference() + from_srs.ImportFromEPSG(from_epsg) + + ogr_geom = ogr.CreateGeometryFromWkb(geom.wkb) + ogr_geom.AssignSpatialReference(from_srs) + + ogr_geom.TransformTo(to_srs) + return loads(ogr_geom.ExportToWkb()) + +def read_db(db, osm_id=0): + # Zero for DB bbox + cur = db.cursor() + if osm_id: + cur.execute("""SELECT way FROM planet_osm_polygon WHERE osm_id = %s;""", (osm_id,)) + else: + cur.execute("""SELECT ST_ConvexHull(ST_Collect(way)) FROM planet_osm_polygon;""") + way = cur.fetchone()[0] + cur.close() + poly = loads(way.decode('hex')) + return project(poly) + +def read_cities(db, osm_id=0): + cur = db.cursor() + if osm_id: + 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,)) + else: + cur.execute("""SELECT ST_Union(way) FROM planet_osm_polygon WHERE place IN ('town', 'city') AND ST_Area(way) < 500*1000*1000;""") + result = cur.fetchone() + poly = loads(result[0].decode('hex')) if result else Polygon() + if osm_id: + 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,)) + else: + 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));""") + result = cur.fetchone() + if result: + poly = poly.union(loads(result[0].decode('hex'))) + return project(poly) + + +if __name__ == "__main__": + try: + mapfile = os.environ['MAPNIK_MAP_FILE'] + except KeyError: + mapfile = os.getcwd() + '/osm.xml' + + default_user = getpass.getuser() + + parser = argparse.ArgumentParser(description='Generate mapnik tiles for OSM polygon') + apg_input = parser.add_argument_group('Input') + apg_input.add_argument("-b", "--bbox", nargs=4, type=float, metavar=('X1', 'Y1', 'X2', 'Y2'), help="generate tiles inside a bounding box") + apg_input.add_argument('-p', '--poly', type=argparse.FileType('r'), help='use a poly file for area') + 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") + apg_input.add_argument("-c", "--cities", type=int, metavar='OSM_ID', help='generate tiles for all towns inside a polygon') + apg_input.add_argument('-l', '--list', type=argparse.FileType('r'), metavar='TILES.LST', help='process tile list') + apg_input.add_argument('-f', '--force', action='store_true', help='force creation of tile, overwrite existing') + apg_output = parser.add_argument_group('Output') + apg_output.add_argument('-t', '--tiledir', metavar='DIR', help='output tiles to directory (default: {0}/tiles)'.format(os.getcwd())) + apg_output.add_argument('--tms', action='store_true', help='write files in TMS order', default=False) + apg_output.add_argument('-m', '--mbtiles', help='generate mbtiles file') + apg_output.add_argument('--name', help='name for mbtiles', default='Test MBTiles') + apg_output.add_argument('--overlay', action='store_true', help='if this layer is an overlay (for mbtiles metadata)', default=False) + apg_output.add_argument('-x', '--export', type=argparse.FileType('w'), metavar='TILES.LST', help='save tile list into file') + 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)) + apg_other = parser.add_argument_group('Settings') + apg_other.add_argument('-s', '--style', help='style file for mapnik (default: {0})'.format(mapfile), default=mapfile) + apg_other.add_argument('--threads', type=int, metavar='N', help='number of threads (default: 2)', default=2) + apg_other.add_argument('-q', '--quiet', dest='verbose', action='store_false', help='do not print any information', default=True) + apg_db = parser.add_argument_group('Database (for poly/cities)') + apg_db.add_argument('-d', '--dbname', metavar='DB', help='database (default: gis)', default='gis') + apg_db.add_argument('--host', help='database host', default='localhost') + apg_db.add_argument('--port', type=int, help='database port', default='5432') + apg_db.add_argument('-u', '--user', help='user name for db (default: {0})'.format(default_user), default=default_user) + apg_db.add_argument('-w', '--password', action='store_true', help='ask for password', default=False) + options = parser.parse_args() + + # check for required argument + if options.bbox == None and options.poly == None and options.cities == None and options.list == None and options.area == None: + parser.print_help() + sys.exit() + + # writer + if options.tiledir: + writer = FileWriter(options.tiledir, options.force) if not options.tms else TMSWriter(options.tiledir) + elif options.mbtiles: + writer = MBTilesWriter(options.mbtiles, options.name, overlay=options.overlay) + elif options.export: + writer = ListWriter(options.export) + else: + writer = FileWriter(os.getcwd() + '/tiles', options.force) if not options.tms else TMSWriter(os.getcwd() + '/tiles') + + # input and process + poly = None + if options.bbox: + b = options.bbox + tpoly = box(b[0], b[1], b[2], b[3]) + poly = tpoly if not poly else poly.intersection(tpoly) + if options.poly: + tpoly = Polygon(poly_parse(options.poly)) + poly = tpoly if not poly else poly.intersection(tpoly) + if options.area != None or options.cities != None: + passwd = "" + if options.password: + passwd = getpass.getpass("Please enter your password: ") + + try: + db = psycopg2.connect(database=options.dbname, user=options.user, password=passwd, host=options.host, port=options.port) + if options.area != None: + tpoly = read_db(db, options.area) + poly = tpoly if not poly else poly.intersection(tpoly) + if options.cities != None: + tpoly = read_cities(db, options.cities) + poly = tpoly if not poly else poly.intersection(tpoly) + db.close() + except Exception, e: + print "Error connecting to database: ", e.pgerror + sys.exit(1) + + if options.list: + generator = ListGenerator(options.list) + elif poly: + generator = PolyGenerator(poly, range(options.zooms[0], options.zooms[1] + 1)) + else: + print "Please specify a region for rendering." + sys.exit() + + render_tiles(generator, options.style, writer, num_threads=options.threads, verbose=options.verbose) + writer.close() + +