generate_tiles.py

   1 #!/usr/bin/python
   2 from math import pi,cos,sin,log,exp,atan
   3 from subprocess import call
   4 import sys, os
   5
   6 DEG_TO_RAD = pi/180
   7 RAD_TO_DEG = 180/pi
   8
   9 def minmax (a,b,c):
  10     a = max(a,b)
  11     a = min(a,c)
  12     return a
  13
  14 class GoogleProjection:
  15     def __init__(self,levels=18):
  16         self.Bc = []
  17         self.Cc = []
  18         self.zc = []
  19         self.Ac = []
  20         c = 256
  21         for d in range(0,levels):
  22             e = c/2;
  23             self.Bc.append(c/360.0)
  24             self.Cc.append(c/(2 * pi))
  25             self.zc.append((e,e))
  26             self.Ac.append(c)
  27             c *= 2
  28
  29     def fromLLtoPixel(self,ll,zoom):
  30          d = self.zc[zoom]
  31          e = round(d[0] + ll[0] * self.Bc[zoom])
  32          f = minmax(sin(DEG_TO_RAD * ll[1]),-0.9999,0.9999)
  33          g = round(d[1] + 0.5*log((1+f)/(1-f))*-self.Cc[zoom])
  34          return (e,g)
  35
  36     def fromPixelToLL(self,px,zoom):
  37          e = self.zc[zoom]
  38          f = (px[0] - e[0])/self.Bc[zoom]
  39          g = (px[1] - e[1])/-self.Cc[zoom]
  40          h = RAD_TO_DEG * ( 2 * atan(exp(g)) - 0.5 * pi)
  41          return (f,h)
  42
  43 from mapnik import *
  44
  45 def render_tiles(bbox, mapfile, tile_dir, minZoom=1,maxZoom=18, name="unknown"):
  46     print "render_tiles(",bbox, mapfile, tile_dir, minZoom,maxZoom, name,")"
  47
  48     if not os.path.isdir(tile_dir):
  49          os.mkdir(tile_dir)
  50
  51     gprj = GoogleProjection(maxZoom+1)
  52     m = Map(2 * 256,2 * 256)
  53     load_map(m,mapfile)
  54     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")
  55
  56     ll0 = (bbox[0],bbox[3])
  57     ll1 = (bbox[2],bbox[1])
  58
  59     for z in range(minZoom,maxZoom + 1):
  60         px0 = gprj.fromLLtoPixel(ll0,z)
  61         px1 = gprj.fromLLtoPixel(ll1,z)
  62
  63         # check if we have directories in place
  64         zoom = "%s" % z
  65         if not os.path.isdir(tile_dir + zoom):
  66             os.mkdir(tile_dir + zoom)
  67         for x in range(int(px0[0]/256.0),int(px1[0]/256.0)+1):
  68             # check if we have directories in place
  69             str_x = "%s" % x
  70             if not os.path.isdir(tile_dir + zoom + '/' + str_x):
  71                 os.mkdir(tile_dir + zoom + '/' + str_x)
  72             for y in range(int(px0[1]/256.0),int(px1[1]/256.0)+1):
  73                 p0 = gprj.fromPixelToLL((x * 256.0, (y+1) * 256.0),z)
  74                 p1 = gprj.fromPixelToLL(((x+1) * 256.0, y * 256.0),z)
  75
  76                 # render a new tile and store it on filesystem
  77                 c0 = prj.forward(Coord(p0[0],p0[1]))
  78                 c1 = prj.forward(Coord(p1[0],p1[1]))
  79
  80                 bbox = Envelope(c0.x,c0.y,c1.x,c1.y)
  81                 bbox.width(bbox.width() * 2)
  82                 bbox.height(bbox.height() * 2)
  83                 m.zoom_to_box(bbox)
  84
  85                 str_y = "%s" % y
  86
  87                 tile_uri = tile_dir + zoom + '/' + str_x + '/' + str_y + '.png'
  88
  89                 exists= ""
  90                 if os.path.isfile(tile_uri):
  91                     exists= "exists"
  92                 else:
  93                     im = Image(512, 512)
  94                     render(m, im)
  95                     view = im.view(128,128,256,256) # x,y,width,height
  96                     view.save(tile_uri,'png')
  97                     command = "convert  -colors 255 %s %s" % (tile_uri,tile_uri)
  98                     call(command, shell=True)
  99
 100                 bytes=os.stat(tile_uri)[6]
 101                 empty= ''
 102                 if bytes == 137:
 103                     empty = " Empty Tile "
 104
 105                 print name,"[",minZoom,"-",maxZoom,"]: " ,z,x,y,"p:",p0,p1,exists, empty
 106
 107 if __name__ == "__main__":
 108     home = os.environ['HOME']
 109     try:
 110         mapfile = os.environ['MAPNIK_MAP_FILE']
 111     except KeyError:
 112         mapfile = home + "/svn.openstreetmap.org/applications/rendering/mapnik/osm-local.xml"
 113     try:
 114         tile_dir = os.environ['MAPNIK_TILE_DIR']
 115     except KeyError:
 116         tile_dir = home + "/osm/tiles/"
 117
 118     #-------------------------------------------------------------------------
 119     #
 120     # Change the following for different bounding boxes and zoom levels
 121     #
 122     # Start with an overview
 123     # World
 124     #bbox = (-180.0,-90.0, 180.0,90.0)
 125
 126     #render_tiles(bbox, mapfile, tile_dir, 0, 5, "World")
 127
 128     #minZoom = 10
 129     #maxZoom = 16
 130     #bbox = (-2, 50.0,1.0,52.0)
 131     #render_tiles(bbox, mapfile, tile_dir, minZoom, maxZoom)
 132
 133     # Caracal
 134     minZoom = 13
 135     maxZoom = 17
 136     bbox = (24.29, 44.10, 24.40, 44.14)
 137     render_tiles(bbox, mapfile, tile_dir, minZoom, maxZoom, "Caracal")
 138
 139     ## Muenchen
 140     #bbox = (11.4,48.07, 11.7,48.22)
 141     #render_tiles(bbox, mapfile, tile_dir, 1, 12 , "Muenchen")
 142
 143     ## Muenchen+
 144     #bbox = (11.3,48.01, 12.15,48.44)
 145     #render_tiles(bbox, mapfile, tile_dir, 7, 12 , "Muenchen+")
 146
 147     ## Muenchen++
 148     #bbox = (10.92,47.7, 12.24,48.61)
 149     #render_tiles(bbox, mapfile, tile_dir, 7, 12 , "Muenchen++")
 150
 151     ## Nuernberg
 152     #bbox=(10.903198,49.560441,49.633534,11.038085)
 153     #render_tiles(bbox, mapfile, tile_dir, 10, 16, "Nuernberg")
 154
 155     ## Karlsruhe
 156     #bbox=(8.179113,48.933617,8.489252,49.081707)
 157     #render_tiles(bbox, mapfile, tile_dir, 10, 16, "Karlsruhe")
 158
 159     ## Karlsruhe+
 160     #bbox = (8.3,48.95,8.5,49.05)
 161     #render_tiles(bbox, mapfile, tile_dir, 1, 16, "Karlsruhe+")
 162
 163     ## Augsburg
 164     #bbox = (8.3,48.95,8.5,49.05)
 165     #render_tiles(bbox, mapfile, tile_dir, 1, 16, "Augsburg")
 166
 167     ## Augsburg+
 168     #bbox=(10.773251,48.369594,10.883834,48.438577)
 169     #render_tiles(bbox, mapfile, tile_dir, 10, 14, "Augsburg+")
 170
 171     ## Europe+
 172     #bbox = (1.0,10.0, 20.6,50.0)
 173     #render_tiles(bbox, mapfile, tile_dir, 1, 11 , "Europe+")
 174