Render metatiles

generate_metatiles.py

#!/usr/bin/env python

from math import pi,cos,sin,log,exp,atan
import sys, os

DEG_TO_RAD = pi/180
RAD_TO_DEG = 180/pi

def minmax (a,b,c):
    a = max(a,b)
    a = min(a,c)
    return a

def mkdir_p(path):
    try:
        os.makedirs(path)
    except OSError as exc:
        if exc.errno == errno.EEXIST:
            pass
        else: raise

def xyz_to_path (x, y, z, base_dir, layer_name,file_extension):
    y_flipped = (1<<z) - 1 - y 
    return "%s/%s/%d/%d/%d.%s" % (base_dir, layer_name, z,x,y_flipped ,file_extension)

def xyz_to_tilecache_path (x, y, z, base_dir, layer_name,file_extension):
    y_flipped = (1<<z) - 1 - y 
    components = ( base_dir,
                   layer_name,
                   "%02d" % z,
                   "%03d" % int( x / 1000000),
                   "%03d" % (int( x / 1000) % 1000),
                   "%03d" % (int(x) % 1000),
                   "%03d" % int(y_flipped / 1000000),
                   "%03d" % (int(y_flipped / 1000) % 1000),
                   "%03d.%s" % (int(y_flipped) % 1000, file_extension)
                   )
    filename = os.path.join( *components )
    return filename


def quad(x,y,z):
    quad = ""
    for i in range(z,0,-1):
        mask = 1 << (i-1)
        cell = 0
        if x & mask != 0:
            cell +=1
        if y & mask != 0:
            cell +=2
        quad += "%d" % cell
    return quad

def xyz_to_hash (x,y,z):
    h = []
    for i in range(0,5):
        h.append(((x & 0x0f) << 4) | (y & 0x0f))
        x >>= 4
        y >>= 4
    h.reverse()
    h.insert(0,z)
    return h


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)


def clamp(x,zoom):
    if x < 0 : 
        return 0
    elif x > 2**zoom:
        return 2**zoom
    else:
        return x

from mapnik2 import *

def render_tiles(bbox, mapfile, base_dir, layer_name, minZoom=1,maxZoom=18):
    print "rendering  tiles (",bbox, mapfile, base_dir, minZoom,maxZoom,")"
    max_meta_step = 16
    tile_size = 256.0
    
    gprj = GoogleProjection(maxZoom+1) 
    
    m = Map(0,0)
    load_map(m,mapfile)
    m.buffer_size=128
    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 +nadgrids=@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)
        
        x0 = clamp(int(px0[0]/tile_size),z)
        x1 = clamp(int(px1[0]/tile_size),z)
        y0 = clamp(int(px0[1]/tile_size),z)
        y1 = clamp(int(px1[1]/tile_size),z)
        
        meta_step_x = min(x1-x0,max_meta_step)
        meta_step_y = min(y1-y0,max_meta_step)
        print>>sys.stderr,px0, px1
        print>>sys.stderr, x0,x1,y0,y1,"META STEP X:",meta_step_x, " Y:",meta_step_y
        m.width  = int((meta_step_x + 1) * tile_size)
        m.height = int((meta_step_y + 1) * tile_size)
        for x in range(x0,x1, meta_step_x):
            for y in range(y0, y1, meta_step_y):
                p0 = gprj.fromPixelToLL((x * tile_size, (y + meta_step_y) * tile_size),z)
                p1 = gprj.fromPixelToLL(((x + meta_step_x) * tile_size, y * tile_size),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]))
                
                #h = xyz_to_hash(x,y,z)
                
                missing_tiles=False
                for xx in range(0, meta_step_x):
                    for yy in range(0,meta_step_y):                
                        tile_path = xyz_to_path (x+xx, y+yy, z, base_dir, layer_name,'png')                            
                        if not os.path.isfile(tile_path):
                            missing_tiles=True
                            break
                if missing_tiles:
                    bbox = Box2d(c0.x,c0.y,c1.x,c1.y)
                    bbox.width(bbox.width() * (meta_step_x + 1.0)/meta_step_x)
                    bbox.height(bbox.height() * (meta_step_y + 1.0)/meta_step_y)
                    m.zoom_to_box(bbox)
                    im = Image(m.width, m.height)
                    render(m, im)
                    
                    for xx in range(0, meta_step_x):
                        for yy in range(0,meta_step_y):
                            
                            tile_path=xyz_to_path (x + xx, y + yy, z, base_dir, layer_name,'png')
    
                            dirname = os.path.dirname(tile_path)
                            if not os.path.isdir(dirname) :
                                mkdir_p(dirname)
                            
                            if not os.path.isfile(tile_path):
                                view = im.view(128 + xx * 256, 128 + yy * 256, 256, 256)
                                view.save(tile_path,'png256')

                    print>>sys.stderr, "[",minZoom,"-",maxZoom,"]: " ,z, x, y, "p:", p0 , p1
                else :
                    print>>sys.stderr, "z=%d x=%d y=%d exists" % (z,x,y)


if __name__ == "__main__":
    MIN_ZOOM=0
    MAX_ZOOM=9
    mapfile = "world.xml"
    base_dir = "/tmp/tiles/"
    layer_name = "world"
    render_tiles( (-180,-90,180,90), mapfile, base_dir, layer_name, 0, 7)
    render_tiles( (-127,45,-100,60), mapfile, base_dir, layer_name, 8, 10)