OpenStreetMap generate georeferenced image

Question

I'm new to Openstreetmap and mapnick,

I'm trying to export map image which will be geo-referenced (So it can be used in other applications)

I've installed osm and mapnik inside ubuntu virtual machine

I've tried using generate_image.py script, but generated image is not equal to the bounding box. My python knowledge is not good enough for me to fix the script.

I've also tried using nik2img.py script using verbose mode, for example:

nik2img.py osm.xml sarajevo.png --srs 900913 --bbox 18.227 43.93 18.511 43.765 --dimensions 10000 10000

and tried using the log bounding box

Step: 11 // --> Map long/lat bbox: Envelope(18.2164733537,43.765,18.5215266463,43.93)

Unfortunately generated image is not equal to the bounding box :(

How can I change scripts so I can georeference generated image? Or do you know an easier way to accomplish this task? Image i'm getting using the http://www.openstreetmap.org/ export is nicely geo-referenced, but it's not big enough :(

Answer 1

I've managed to change generate_tiles.py to generate 1024x1024 images together with correct bounding box

Changed script is available bellow

#!/usr/bin/python
from math import pi,cos,sin,log,exp,atan
from subprocess import call
import sys, os
from Queue import Queue
import mapnik
import threading

DEG_TO_RAD = pi/180
RAD_TO_DEG = 180/pi

# Default number of rendering threads to spawn, should be roughly equal to number of CPU cores available
NUM_THREADS = 4


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 = 1024
        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, mapfile, q, printLock, maxZoom):
        self.tile_dir = tile_dir
        self.q = q
        self.m = mapnik.Map(1024, 1024)
        self.printLock = printLock
        # Load style XML
        mapnik.load_map(self.m, 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(maxZoom+1)

    def render_tile(self, tile_uri, x, y, z):
        # Calculate pixel positions of bottom-left & top-right
        p0 = (x * 1024, (y + 1) * 1024)
        p1 = ((x + 1) * 1024, y * 1024)

        # 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 = 1024
        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)
        im.save(tile_uri, 'png256')
    print "Rendered: ", tile_uri, "; ", l0 , "; ", l1


    # Write geo coding informations 
    file = open(tile_uri[:-4] + ".tab", 'w')
    file.write("!table\n")
    file.write("!version 300\n")
    file.write("!charset WindowsLatin2\n")

    file.write("Definition Table\n")
    file.write("  File \""+tile_uri[:-4]+".jpg\"\n")
    file.write("  Type \"RASTER\"\n")
    file.write("  ("+str(l0[0])+","+str(l1[1])+") (0,0) Label \"Pt 1\",\n")
    file.write("  ("+str(l1[0])+","+str(l1[1])+") (1023,0) Label \"Pt 2\",\n")
    file.write("  ("+str(l1[0])+","+str(l0[1])+") (1023,1023) Label \"Pt 3\",\n")
    file.write("  ("+str(l0[0])+","+str(l0[1])+") (0,1023) Label \"Pt 4\"\n")
    file.write("  CoordSys Earth Projection 1, 104\n")
    file.write("  Units \"degree\"\n")

    file.close()



    def loop(self):
        while True:
            #Fetch a tile from the queue and render it
            r = self.q.get()
            if (r == None):
                self.q.task_done()
                break
            else:
                (name, tile_uri, x, y, z) = r

            exists= ""
            if os.path.isfile(tile_uri):
                exists= "exists"
            else:
                self.render_tile(tile_uri, x, y, z)
            bytes=os.stat(tile_uri)[6]
            empty= ''
            if bytes == 103:
                empty = " Empty Tile "
            self.printLock.acquire()
            print name, ":", z, x, y, exists, empty
            self.printLock.release()
            self.q.task_done()



def render_tiles(bbox, mapfile, tile_dir, minZoom=1,maxZoom=18, name="unknown", num_threads=NUM_THREADS):
    print "render_tiles(",bbox, mapfile, tile_dir, minZoom,maxZoom, name,")"

    # Launch rendering threads
    queue = Queue(32)
    printLock = threading.Lock()
    renderers = {}
    for i in range(num_threads):
        renderer = RenderThread(tile_dir, mapfile, queue, printLock, maxZoom)
        render_thread = threading.Thread(target=renderer.loop)
        render_thread.start()
        #print "Started render thread %s" % render_thread.getName()
        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)

        # 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]/1024.0),int(px1[0]/1024.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]/1024.0),int(px1[1]/1024.0)+1):
                # Validate x co-ordinate
                if (y < 0) or (y >= 2**z):
                    continue
                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)
                queue.put(t)

    # 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']
    try:
        mapfile = "/home/emir/bin/mapnik/osm.xml" #os.environ['MAPNIK_MAP_FILE']
    except KeyError:
        mapfile = "/home/emir/bin/mapnik/osm.xml"
    try:
        tile_dir = os.environ['MAPNIK_TILE_DIR']
    except KeyError:
        tile_dir = home + "/osm/tiles/"

    if not tile_dir.endswith('/'):
        tile_dir = tile_dir + '/'

    #-------------------------------------------------------------------------
    #
    # Change the following for different bounding boxes and zoom levels
    #

    #render sarajevo at 16 zoom level
    bbox = (18.256, 43.785, 18.485, 43.907)
    render_tiles(bbox, mapfile, tile_dir, 16, 16, "World")

Answer 2

Try Maperitive's export-bitmap command, it generates various georeferencing sidecar files (worldfile, KML, OziExplorer .MAP file).