I'm trying to generate png tiles to overlay them on Google Map. Genenriruyu using mnapnik. But the image of the tile after the 15 scale is not generated. Any ideas? Here is the script generator:
DEG_TO_RAD = pi/180
RAD_TO_DEG = 180/pi
NUM_THREADS = 1
def minmax(a, b, c):
a = max(a, b)
a = min(a, c)
return a
class GoogleProjection:
def __init__(self, levels=18):
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 RenderThread:
def __init__(self, tile_dir, q, maxZoom):
self.tile_dir = tile_dir
self.q = q
# Create Map
self.my_map = mapnik.Map(256, 256)
# create style
s = mapnik.Style()
r = mapnik.Rule()
# add polygon
polygon_symbolizer = mapnik.PolygonSymbolizer()
polygon_symbolizer.fill = mapnik.Color('#ffffff')
r.symbols.append(polygon_symbolizer)
line_symbolizer = mapnik.LineSymbolizer()
line_symbolizer.stroke = mapnik.Color('#000000')
line_symbolizer.stroke_width = 0.5
r.symbols.append(line_symbolizer)
s.rules.append(r)
self.my_map.append_style('My Style', s)
layer = mapnik.Layer('world')
layer.styles.append('My Style')
layer.datasource = mapnik.PostGIS(**polygon_db_params)
self.my_map.layers.append(layer)
# Obtain <Map> projection
self.prj = mapnik.Projection(self.my_map.srs)
print 'PRG ', self.prj
# Projects between tile pixel co-ordinates and LatLong (EPSG:4326)
self.tileproj = GoogleProjection(maxZoom+1)
def render_tile(self, tile_uri, x, y, z):
# Calculate pixel positions of bottom-left & top-right
p0 = (x * 256, (y + 1) * 256)
p1 = ((x + 1) * 256, y * 256)
# 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 = 256
self.my_map.zoom_to_box(bbox)
if self.my_map.buffer_size < 128:
self.my_map.buffer_size = 128
# Render image with default Agg renderer
im = mapnik.Image(render_size, render_size)
mapnik.render(self.my_map, im)
im.save(tile_uri, 'png256')
def loop(self):
while True:
#Fetch a tile from the queue and render it
r = self.q.get()
if r is None:
self.q.task_done()
break
else:
(name, tile_uri, x, y, z) = r
if os.path.isfile(tile_uri):
pass
else:
self.render_tile(tile_uri, x, y, z)
self.q.task_done()
def render_tiles(bbox, tile_dir, minZoom=1, maxZoom=18,
name="unknown", num_threads=NUM_THREADS, tms_scheme=False):
# Launch rendering threads
queue = Queue(32)
renderers = {}
for i in range(num_threads):
renderer = RenderThread(tile_dir, queue, maxZoom)
render_thread = threading.Thread(target=renderer.loop)
render_thread.start()
renderers[i] = render_thread
if not os.path.isdir(tile_dir):
os.mkdir(tile_dir)
gprj = GoogleProjection(maxZoom+1)
ll0 = (bbox[0], bbox[3])
ll1 = (bbox[2], bbox[1])
for z in range(minZoom, maxZoom + 1):
px0 = gprj.fromLLtoPixel(ll0, z)
px1 = gprj.fromLLtoPixel(ll1, z)
print z, " >> convert to pixel ", px0, ' ', px1
# check if we have directories in place
zoom = "%s" % z
if not os.path.isdir(tile_dir + zoom):
os.mkdir(tile_dir + zoom)
for x in range(int(px0[0]/256.0), int(px1[0]/256.0)+1):
# Validate x co-ordinate
if (x < 0) or (x >= 2**z):
continue
# check if we have directories in place
str_x = "%s" % x
if not os.path.isdir(tile_dir + zoom + '/' + str_x):
os.mkdir(tile_dir + zoom + '/' + str_x)
for y in range(int(px0[1]/256.0), int(px1[1]/256.0)+1):
# Validate x co-ordinate
if (y < 0) or (y >= 2**z):
continue
# flip y to match OSGEO TMS spec
if tms_scheme:
str_y = "%s" % ((2**z-1) - y)
else:
str_y = "%s" % y
tile_uri = tile_dir + zoom + '/' + str_x + '/' + str_y + '.png'
# Submit tile to be rendered into the queue
t = (name, tile_uri, x, y, z)
try:
queue.put(t)
except KeyboardInterrupt:
raise SystemExit("Ctrl-c detected, exiting...")
# 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()
if __name__ == "__main__":
home = os.environ['HOME']
tile_dir = home + "/osm/tiles/"
# Change the following for different bounding boxes and zoom levels
bbox = (50.18117, 53.21847, 50.3, 53.22)
render_tiles(bbox, tile_dir, 12, 15, "World")
the coordinates in the database are stored in the format of WGS84. The Render goes from scale 1 to scale 18. on the 18th scale there is a duplication of tiles. Example:
set:
self.my_map = mapnik.Map(256, 256)
self.my_map.srs="+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"