pythongoogle-mapsmapnik

Incorrectly render png tiles


I'm trying to generate png tiles to overlay them on Google Map. Genenriruyu using mnapnik. But the image of the tile after the 15 scale is not generated. Any ideas? Here is the script generator:

DEG_TO_RAD = pi/180
RAD_TO_DEG = 180/pi

NUM_THREADS = 1

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 = 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

class RenderThread:
    def __init__(self, tile_dir, q, maxZoom):
        self.tile_dir = tile_dir
        self.q = q

        # Create Map
        self.my_map = mapnik.Map(256, 256)

        # create style
        s = mapnik.Style()
        r = mapnik.Rule()

        # add polygon
        polygon_symbolizer = mapnik.PolygonSymbolizer()
        polygon_symbolizer.fill = mapnik.Color('#ffffff')
        r.symbols.append(polygon_symbolizer)

        line_symbolizer = mapnik.LineSymbolizer()
        line_symbolizer.stroke = mapnik.Color('#000000')
        line_symbolizer.stroke_width = 0.5

        r.symbols.append(line_symbolizer)
        s.rules.append(r)
        self.my_map.append_style('My Style', s)

        layer = mapnik.Layer('world')
        layer.styles.append('My Style')
        layer.datasource = mapnik.PostGIS(**polygon_db_params)
        self.my_map.layers.append(layer)

        # Obtain <Map> projection
        self.prj = mapnik.Projection(self.my_map.srs)
        print 'PRG ', self.prj
        # 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 * 256, (y + 1) * 256)
        p1 = ((x + 1) * 256, y * 256)

        # 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 = 256
        self.my_map.zoom_to_box(bbox)

        if self.my_map.buffer_size < 128:
            self.my_map.buffer_size = 128

        # Render image with default Agg renderer
        im = mapnik.Image(render_size, render_size)
        mapnik.render(self.my_map, im)
        im.save(tile_uri, 'png256')

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

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

    # Launch rendering threads
    queue = Queue(32)
    renderers = {}
    for i in range(num_threads):
        renderer = RenderThread(tile_dir, queue, maxZoom)
        render_thread = threading.Thread(target=renderer.loop)
        render_thread.start()
        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)
        print z, " >> convert to pixel ", px0, '  ', px1

        # 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]/256.0), int(px1[0]/256.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]/256.0), int(px1[1]/256.0)+1):
                # Validate x co-ordinate
                if (y < 0) or (y >= 2**z):
                    continue
                # flip y to match OSGEO TMS spec
                if tms_scheme:
                    str_y = "%s" % ((2**z-1) - y)
                else:
                    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)
                try:
                    queue.put(t)
                except KeyboardInterrupt:
                    raise SystemExit("Ctrl-c detected, exiting...")
    # 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']
    tile_dir = home + "/osm/tiles/"

    # Change the following for different bounding boxes and zoom levels
    bbox = (50.18117, 53.21847, 50.3, 53.22)
    render_tiles(bbox, tile_dir, 12, 15, "World")

the coordinates in the database are stored in the format of WGS84. The Render goes from scale 1 to scale 18. on the 18th scale there is a duplication of tiles. Example:

incorrect_image

correct_image

incorrect_image


Solution

  • set:

    self.my_map = mapnik.Map(256, 256)
    self.my_map.srs="+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"