Hello everyone.
First of all, I want to share with my two scripts for Mapnik rendering. First
is a modified version of generate_tiles.py to use multithread rendering. Second
is a modified generate_image.py to generate large bitmap posters like [1,2]. It
uses tiled rendering so it's possible to render very large areas with high
DPI. Only limitation is that whole poster must fit into computer's memory
during merging and texts drawing because I don't know any image rendering
library which can draw to memory mapped files or such. See scripts' code and
comments or ask for detailed information.
Seconly, I would like to ask how should I properly extend an OpenStreetMap
attribution on the [1,2] to include my name and to mention an OpenTrackMap
project [3]. Could you give me an advice about this?
Thank you in advance.
[1] - http://blackhex.no-ip.org/raw-attachment/wiki/Files/poster_preview.png -
OpenTrackMap poster, scaled down to 25% (13.4 MB)
[2] - http://blackhex.no-ip.org/raw-
attachment/wiki/Files/poster_small_preview.png - OpenTrackMap poster, scaled
down to 10% (2.4 MB)
[3] - http://opentrackmap.no-ip.org/ - OpenTrackMap project homepage.
--
Ing. Radek Bartoň
Faculty of Information Technology
Department of Computer Graphics and Multimedia
Brno University of Technology
E-mail: [email protected]
Web: http://blackhex.no-ip.org
Jabber: [email protected]
#!/usr/bin/python
# -*- coding: utf-8 -*-
#
# Generates a single large PNG image for a UK bounding box
# Tweak the lat/lon bounding box (ll) and image dimensions
# to get an image of arbitrary size.
#
# To use this script you must first have installed mapnik
# and imported a planet file into a Postgres DB using
# osm2pgsql.
#
# Note that mapnik renders data differently depending on
# the size of image. More detail appears as the image size
# increases but note that the text is rendered at a constant
# pixel size so will appear smaller on a large image.
from math import pi,cos,sin,log,exp,atan
from mapnik import *
import sys, os
if __name__ == "__main__":
# Use Mapnik template from MAPNIK_MAP_FILE environment variable or osm.xml
# file.
try:
map_file = os.environ['MAPNIK_MAP_FILE']
except KeyError:
map_file = "osm.xml"
# Configuration options
poster_filename = "poster.png" # Name of output PNG file.
lonlat_bottom_left = Coord(12.10, 48.55) # Bottom-left corner of area of interest in Lon/Lat coordinates.
lonlat_top_right = Coord(18.86, 51.06) # Top-right corner of area of interest in Lon/Lat coordinates.
x_size = 14040 # A0 width at 300 DPI
y_size = 9930 # A0 height at 300 DPI
#x_size = 28080 # A0 width at 600 DPI
#y_size = 19860 # A0 height at 600 DPI
x_parts = 8 # Number of tiles in x-direction.
y_parts = 8 # Number of tile in y-direction
overlap = 32 # Number of pixels of tile's overlapping.
borders = (400, 600, 400, 800) # Left, bottom, right and top white border of the poster.
frame_line_width = 10 # Line width of frame around map.
title = "OpenStreetMap - CTC Hiking Tracks Network - 30 September 2009" # Text of a title on the top.
title_font_size = 300 # Font size of the title.
title_offset = 700 # Offset of title bottom from the top of the poster.
attribution = "\302\251 OpenStreetMap & contributors (http://www.openstreetmap.org/), CC-BY-SA" # Text of an attribtuion on the bottom.
attribution_font_size = 100 # Font size of the attribution,
attribution_offset = 400 # Offset of attribution bottom from the bottom of the poster.
# Compute image size for map.
x_size = x_size - borders[0] - borders[2]
y_size = y_size - borders[1] - borders[3]
# Compute rendered area size in target projection.
projection = Projection("+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 +nadgri...@null +no_defs +over")
target_bottom_left = projection.forward(lonlat_bottom_left)
target_top_right = projection.forward(lonlat_top_right)
# Modify area to respect aspect ratio of image.
image_aspect = x_size / float(y_size)
target_center = Coord((target_bottom_left.x + target_top_right.x) * 0.5,
(target_bottom_left.y + target_top_right.y) * 0.5)
target_width = target_top_right.x - target_bottom_left.x
target_height = target_top_right.y - target_bottom_left.y
target_aspect = target_width / float(target_height)
if image_aspect > target_aspect:
target_bottom_left.x = target_center.x - (target_height * image_aspect * 0.5)
target_top_right.x = target_center.x + (target_height * image_aspect * 0.5)
else:
target_bottom_left.y = target_center.y - (target_width / image_aspect * 0.5)
target_top_right.y = target_center.y + (target_width / image_aspect * 0.5)
# Update target dimensions.
target_width = target_top_right.x - target_bottom_left.x
target_height = target_top_right.y - target_bottom_left.y
target_aspect = target_width / float(target_height)
# Inform about visible area and load map data.
lonlat_bottom_left = projection.inverse(target_bottom_left)
lonlat_top_right = projection.inverse(target_top_right)
print "Loading map of area (%f lon %f lat) - (%f lon %f lat)..." % (
lonlat_bottom_left.x, lonlat_bottom_left.y, lonlat_top_right.x,
lonlat_top_right.y),; sys.stdout.flush()
part_x_size = x_size / x_parts
part_y_size = y_size / y_parts
m = Map(part_x_size + 2 * overlap, part_y_size + 2 * overlap)
load_map(m, map_file)
print "done."
# Render tiles.
start_x = target_bottom_left.x
step_x = target_width / x_parts
start_y = target_bottom_left.y
step_y = target_height / y_parts
target_overlap = overlap * (target_height / y_size)
for x in range(x_parts):
for y in range(y_parts):
c0 = Coord(start_x + (x * step_x) - target_overlap, start_y + (y *
step_y) - target_overlap)
c1 = Coord(start_x + ((x + 1) * step_x) + target_overlap, start_y +
((y + 1) * step_y) + target_overlap)
bbox = Envelope(c0.x, c0.y, c1.x, c1.y)
m.zoom_to_box(bbox)
# Render part to image.
print "Rendering part (%d, %d)..." % (x, y),; sys.stdout.flush()
image = Image(part_x_size + 2 * overlap, part_y_size + 2 * overlap)
render(m, image)
view = image.view(overlap, overlap, part_x_size, part_y_size)
view.save('part%.2d%.2d.png' % (x, y), 'png')
print "done."
from cairo import ImageSurface, Context, FORMAT_RGB24, FONT_SLANT_NORMAL, \
FONT_WEIGHT_BOLD
# Join them to single image.
print "Merging parts...",; sys.stdout.flush()
x_size = part_x_size * x_parts + borders[0] + borders[2]
y_size = part_y_size * y_parts + borders[1] + borders[3]
surface = ImageSurface(FORMAT_RGB24, x_size, y_size)
context = Context(surface)
context.set_source_rgb(255, 255, 255)
context.paint()
for x in range(x_parts):
for y in range(y_parts):
image_surface = ImageSurface.create_from_png("part%.2d%.2d.png" % (x, y))
context.set_source_surface(image_surface, x * part_x_size + borders[0],
(y_parts - y - 1) * part_y_size + borders[3])
context.paint()
os.remove("part%.2d%.2d.png" % (x, y))
# Draw border
context.set_source_rgb(0, 0, 0)
context.set_line_width(frame_line_width)
context.rectangle(borders[0], borders[3], x_size - borders[0] - borders[2],
y_size - borders[1] - borders[3])
context.stroke()
# Draw title.
context.select_font_face("Nimbus Sans", FONT_SLANT_NORMAL, FONT_WEIGHT_BOLD)
context.set_font_size(title_font_size)
_, _, width, height, _, _ = context.text_extents(title)
context.move_to((x_size - width) * 0.5, title_offset)
context.show_text(title)
# Draw attribution.
context.set_font_size(attribution_font_size)
_, _, width, height, _, _ = context.text_extents(attribution)
context.move_to((x_size - width) * 0.5, y_size - attribution_offset)
context.show_text(attribution)
# Save to target file.
surface.write_to_png(poster_filename)
print "done."
#!/usr/bin/python
# -*- coding: utf-8 -*-
from math import pi,cos,sin,log,exp,atan
from subprocess import call
from threading import Thread, Lock
from datetime import datetime
from time import sleep
import sys, os
os.environ['LD_LIBRARY_PATH'] = '/home/blackhex/local/Mapnik/lib'
path = ['/home/blackhex/local/Mapnik/lib/python2.6/site-packages']
path.extend(sys.path)
sys.path = path
DEG_TO_RAD = pi/180
RAD_TO_DEG = 180/pi
fsLock = Lock()
stdoutLock = Lock()
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)
from mapnik import *
def render_tiles(mapFile, threadId, bbox, tileSQDir, tileHQDir, minZoom = 1,
maxZoom = 18, name = "unknown"):
global fsLock, stdoutLock, finish
print "Rendering: ", bbox, minZoom, maxZoom, name
# Setup Mapnik
gprj = GoogleProjection(maxZoom + 1)
m = Map(2 * 256, 2 * 256)
load_map(m, mapFile)
prj = Projection("+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 +nadgri...@null +no_defs +over")
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)
# Check if we have directories in place.
zoom = "%s" % z
fsLock.acquire()
if not os.path.exists(os.path.join(tileSQDir, zoom)):
os.mkdir(os.path.join(tileSQDir, zoom))
if not os.path.exists(os.path.join(tileHQDir, zoom)):
os.mkdir(os.path.join(tileHQDir, zoom))
fsLock.release()
for x in range(int(px0[0]/256.0),int(px1[0]/256.0)+1):
# Check if we have directories in place
str_x = "%s" % x
fsLock.acquire()
if not os.path.exists(os.path.join(tileSQDir, zoom, str_x)):
os.mkdir(os.path.join(tileSQDir, zoom, str_x))
if not os.path.exists(os.path.join(tileHQDir, zoom, str_x)):
os.mkdir(os.path.join(tileHQDir, zoom, str_x))
fsLock.release()
for y in range(int(px0[1]/256.0),int(px1[1]/256.0)+1):
start_time = datetime.now()
p0 = gprj.fromPixelToLL((x * 256.0, (y+1) * 256.0),z)
p1 = gprj.fromPixelToLL(((x+1) * 256.0, y * 256.0),z)
# render a new tile and store it on filesystem
c0 = prj.forward(Coord(p0[0],p0[1]))
c1 = prj.forward(Coord(p1[0],p1[1]))
bbox = Envelope(c0.x,c0.y,c1.x,c1.y)
bbox.width(bbox.width() * 2)
bbox.height(bbox.height() * 2)
m.zoom_to_box(bbox)
str_y = "%s" % y
tileHQFile = os.path.join(tileHQDir, zoom, str_x, str_y + '.png')
tileSQFile = os.path.join(tileSQDir, zoom, str_x, str_y + '.png')
exists = ""
if os.path.exists(tileHQFile):
exists = "exists"
else:
im = Image(512, 512)
render(m, im)
view = im.view(128,128,256,256) # x,y,width,height
fsLock.acquire()
# We should check file existance aggain for atomicity.
if not os.path.exists(tileHQFile):
view.save(tileHQFile,'png')
fsLock.release()
# This is not thread safe.
if not os.path.exists(tileSQFile):
call("convert -colors 255 %s %s" % (tileHQFile,
tileSQFile), shell = True)
bytes = os.stat(tileHQFile)[6]
empty = ""
if bytes == 137:
empty = "empty tile"
# Determine time of rendering.
end_time = datetime.now()
time_delta = end_time - start_time
stdoutLock.acquire()
print "Thread #" + str(threadId), time_delta, name, "[", \
minZoom, "-", maxZoom, "]: ", z, x, y, "p:", p0, p1, exists, \
empty
stdoutLock.release()
# Terminate thread if requested to finish.
if finish:
return
class RenderThread(Thread):
def __init__(self, mapFile, threadId, bbox, tileSQDir, tileHQDir, minZoom,
maxZoom):
Thread.__init__(self)
self.mapFile = mapFile
self.threadId = threadId
self.bbox = bbox
self.tileSQDir = tileSQDir
self.tileHQDir = tileHQDir
self.minZoom = minZoom
self.maxZoom = maxZoom
def run(self):
render_tiles(self.mapFile, self.threadId, self.bbox, self.tileSQDir,
self.tileHQDir, self.minZoom, self.maxZoom, "OpenStreetMap")
if __name__ == "__main__":
try:
mapFile = os.environ['MAPNIK_MAP_FILE']
except KeyError:
sys.stderr.write('MAPNIK_MAP_FILE environment variable is not set!')
exit(1)
try:
tileSQDir = os.environ['MAPNIK_TILE_SQ_DIR']
except KeyError:
sys.stderr.write('MAPNIK_TILE_SQ_DIR environment variable is not set!')
exit(1)
try:
tileHQDir = os.environ['MAPNIK_TILE_HQ_DIR']
except KeyError:
sys.stderr.write('MAPNIK_TILE_HQ_DIR environment variable is not set!')
exit(1)
try:
numThreads = int(os.environ['MAPNIK_NUM_THREADS'])
except KeyError:
numThreads = 1
# Czech Republic
minZoom = 0
maxZoom = 14
bbox = (12.10, 48.55, 18.86, 51.06)
# Create target directories.
if not os.path.exists(tileSQDir):
os.mkdir(tileSQDir)
if not os.path.exists(tileHQDir):
os.mkdir(tileHQDir)
# Create and setup threads.
finish = False
threads = []
threadLon = (bbox[2] - bbox[0]) / numThreads
for threadId in range(numThreads):
threadBbox = (bbox[0] + threadLon * threadId, bbox[1], bbox[0] +
threadLon * (threadId + 1), bbox[3])
thread = RenderThread(mapFile, threadId, threadBbox, tileSQDir, tileHQDir,
minZoom, maxZoom)
threads.append(thread)
# Start the threads,
for thread in threads:
thread.start()
# Wait for threads to finish.
try:
while not finish:
for thread in threads:
thread.join(1)
my_finish = True
for thread in threads:
if thread.isAlive():
my_finish = False
break;
finish = my_finish
except KeyboardInterrupt:
finish = True
for thread in threads:
thread.join()
_______________________________________________
talk mailing list
[email protected]
http://lists.openstreetmap.org/listinfo/talk
