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: 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")
    

