2011-02-15 71 views
0

我是新來的OpenStreetMap和mapnick,OpenStreetMap的生成地理參考圖像

我試圖出口將被地理參考 地圖圖像(所以它可以在其他應用程序中使用)

我已經在Ubuntu虛擬機中安裝了osm和mapnik

我試過使用generate_image.py腳本,但生成的圖像不等於邊界框。我的Python知識不足以修復腳本。

我使用nik2img.py腳本中使用詳細模式,例如也嘗試:

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

,並使用日誌邊框

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

不幸的是生成的圖像試過不等於邊界框:(

如何更改腳本,以便我可以對生成的圖像進行地理參考? 或者您是否知道一個更簡單的方法來完成此任務? 圖片我使用http://www.openstreetmap.org/出口越來越是很好的地理參照,但它不夠大:(

回答

0

我已經成功地改變generate_tiles.py生成1024×1024的圖像加上正確的邊框

更改腳本可用波紋管

#!/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") 
0

嘗試Maperitive's export-bitmap命令,它會產生各種地理參考附屬文件 (worldfile,KML,OziExplorer .MAP文件)。