Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

enter image description here

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?

share|improve this question
add comment

2 Answers

up vote 1 down vote accepted

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.
share|improve this answer
1  
Now that Mintx has pointed the coordinate value problem out (I can see it now!), the 2*20million value in def _unitsPerPixel(zoomLevel) also has one too many zeroes. –  mkennedy May 20 '13 at 21:47
    
Good catch! That was not an easy one to find. Thank you! –  Below the Radar May 20 '13 at 21:49
    
@mkennedy btw, great to meet you at APSG! :) –  Mintx May 22 '13 at 20:54
add comment

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?

share|improve this answer
    
I use Mapnik in GeoDjango, im not sure about adding another server like TileStache. I would like to have its features but dont know if its good to bind it inside geodjango. –  Below the Radar May 9 '13 at 19:51
    
Anyway im not using Mapnik xml and tilestache seems not to support mapnik.Map object directly –  Below the Radar May 19 '13 at 14:38
1  
creating mapnik xml from a python mapnik.Map object is easy: gist.github.com/springmeyer/5637596 –  Dane May 24 '13 at 17:24
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.