Tile script (Mapnik) with OpenLayers - Problem with projection (900913)

Question

I use Mapnik python bindings in a geodjango application to render tiles in openlayers with datasource stored in postGIS in 900913 projection. The feature seems to be rendered correctly in the spherical mercator projection, but at a very high zoomLevel in the very left corner of the world map so the features seems really tiny.

EDIT: After several test the problem seems to be in the first lines of the script where minLong, minLat, maxLong, and maxLat are defined. Someone has a similar tiling script in 900913?

This is weird because the bounds of the features are:

-8879435.28405854, 5620161.13498339, -6337967.06663213, 8998594.89464899

but they are rendered around this position in the map:

lat: -19955850.825909 lon: -19952809.28254

There is the script

MAX_ZOOM_LEVEL = 20
TILE_WIDTH = 256
TILE_HEIGHT = 256

dbSettings = settings.DATABASES['default']

def tile(request, version, shapefile_id, zoom, x, y):

    try:
        if version != "1.0":
            raise Http404
        try:
            shapefile = Shapefile.objects.get(id=shapefile_id, created_by=request.user)
        except Shapefile.DoesNotExist:
            raise Http404
        zoom = int(zoom)
        x = int(x)
        y = int(y)
        if zoom < 0 or zoom > MAX_ZOOM_LEVEL:
            raise Http404
        xExtent = _unitsPerPixel(zoom) * TILE_WIDTH
        yExtent = _unitsPerPixel(zoom) * TILE_HEIGHT
        minLong = x * xExtent - 200037508.34
        minLat = y * yExtent - 200037508.34
        maxLong = minLong + xExtent
        maxLat = minLat + yExtent

        #create de mapnik.map object
        map = mapnik.Map(TILE_WIDTH, TILE_HEIGHT, "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +units=m +k=1.0 +nadgrids=@null +no_defs")
        map.background = mapnik.Color("#7391ad")

        #defining the feature layer
        geometryField = utils.calcGeometryField(shapefile.geom_type)
        query = '(select ' + geometryField + ', id from "shapefile_feature" where' + ' shapefile_id=' + str(shapefile.id) + ') as geom'

        datasource = mapnik.PostGIS(user=dbSettings['USER'],
                        password=dbSettings['PASSWORD'],
                        dbname=dbSettings['NAME'],
                        table=query,
                        srid=900913,
                        geometry_field=geometryField,
                        geometry_table='"shapefile_feature"',
                        simplify_geometries=True)

        featureLayer = mapnik.Layer("featureLayer")
        featureLayer.srs = "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +no_defs +over"
        featureLayer.datasource = datasource
        featureLayer.styles.append("featureLayerStyle")

        #defining the feature layer styles
        rule = mapnik.Rule()
        if shapefile.geom_type in ["Point", "MultiPoint"]:
            rule.symbols.append(mapnik.PointSymbolizer())
        elif shapefile.geom_type in ["LineString", "MultiLineString"]:
            rule.symbols.append(mapnik.LineSymbolizer(mapnik.Color("#000000"), 0.5))
        elif shapefile.geom_type in ["Polygon", "MultiPolygon"]:
            rule.symbols.append(mapnik.PolygonSymbolizer(mapnik.Color("#f7edee")))
            rule.symbols.append(mapnik.LineSymbolizer(mapnik.Color("#000000"), 0.5))
        style = mapnik.Style()
        style.rules.append(rule)

        #add new feature to the map
        map.append_style("featureLayerStyle", style)
        map.layers.append(featureLayer)
        map.srs = "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +units=m +k=1.0 +nadgrids=@null +no_defs"
        #rendering the map tile
        map.zoom_to_box(mapnik.Box2d(minLong, minLat, maxLong, maxLat))
        image = mapnik.Image(TILE_WIDTH, TILE_HEIGHT)
        mapnik.render(map, image)
        imageData = image.tostring('png')
        return HttpResponse(imageData, image)

    except:
        traceback.print_exc()
        return HttpResponse("")

def _unitsPerPixel(zoomLevel):
    return 400075016.68 / math.pow(2, zoomLevel)

In openlayers.html

map = new OpenLayers.Map('map',
          {maxResolution: 156543.0399,
            projection: new OpenLayers.Projection("EPSG:900913"),
            units: 'm',
            maxExtent: new OpenLayers.Bounds(-20037508.34,-20037508.34,20037508.34,20037508.34),
            numZoomLevels: 21});
tiledLayer = new OpenLayers.Layer.TMS('TMS',
                "{{ tmsURL }}",
                {serviceVersion: "1.0",
                layername: "{{ shapefile.id }}",
                type: 'png',
                projection: map.projection});
map.addLayer(tiledLayer);

What am I missging? why the features are not rendered at their normal sizes and positions?

Answer 1

I can spot two things:

1. Change every instance of EPSG:900913 to EPSG:3857 -- EPSG:900913 (Or EPSG:GOOGLE if you squint hard enough) never existed in the EPSG database. It's a long story, but EPSG:3857 is the correct code.
2. It looks like you have one too many zeroes when you initialize your minLat and minLong variables. Compare it to your new OpenLayers.Bounds in the HTML.

Answer 2

Many different tools implement conversions between tile z/x/y and bounding boxes for a given projection. This problem is not specific to Mapnik at all. Why not use a library like TileStache, Invar, or MapProxy?