550 lines
20 KiB
Python
Executable File
550 lines
20 KiB
Python
Executable File
#!/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 <Map> 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()
|
|
|
|
|