3

i am trying to render a map with Mapnik via Python. Therefore I first created a style file (as XML) and my code looks so far:

m = mapnik.Map(400, 400)
style = 'osm.xml'
mapnik.load_map(m, style)

box = mapnik.Box2d(347950, 5207105, 425404, 5316784)
m.zoom_to_box(box)

mapnik.render_to_file(m, 'map.png')

The problem is that the map shows the wrong part of the world. I think somewhere in the conversion is an error.

I tried the following values:

Lat/Lon: 47/7 ==> X/Y: 32T 347950mE 5207105mN

Lat/Lon: 48/8 ==> X/Y: 32U 425404mE 5316784mN

As I enter these lat/lon values on google maps it shows the mid east part of france southwest part of germany. Calculating the map with the calculated values i get the coast of france, especialle the city named "Port Vendres" in the lower left corner.

What am i doing wrong here?

The values in the database seem to bee in mercator form and the data was imported using osm2pgsql with no special options.

Update: I printed out m.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 +wktext +no_defs +over

and the scale_denominator(): -8.92857142857

Update 2: I read the documentation of osm2pgsql and it says that it imports the data as spherical mercator by default. So i assume that the values are stored in that format and the projection is correct.

  • 1
    I think you're confusing UTM zone 32 North (32632) with "web Mercator" (3857). If you want UTM, the Proj4 string is all wrong. – mkennedy Feb 13 '15 at 18:32
  • Ah ok. Thanks for clearing that. But is it possible to only calculate in Mercator or do I need the way over UTM? Sorry for the beginner question ... – Freddy Feb 13 '15 at 18:45
1

I guess most tools around the rendering stack rely heavily on WGS84 and so it might be not a good idea to alter the projection of the metadata. Instead I would recommend to reproject the rendered image itself?

Doing it on the data level seems to be very tricky: http://wiki.openstreetmap.org/wiki/Mapnik_GB_Projection

  • Thank you very much for this comment. It gave me the tip to look at the projection in detail. – Freddy Mar 2 '15 at 14:38
1

I finally solved the problem with the following code:

I left the code (almost) as it was before (especially the proj4) and as it is in the question. I then added the following:

lat_min = 47.975
lon_min = 7.8

lat_max = 48.0
lon_max = 7.875

x_min = 6378137.0 * num * lon_min

a = lat_min * 0.017453292519943295;
y_min =  3189068.5 * math.log((1.0 + math.sin(a)) / (1.0 - math.sin(a)))

print x_min
print y_min

x_max = 6378137.0 * num * lon_max
a = lat_max * 0.017453292519943295;
y_max =  3189068.5 * math.log((1.0 + math.sin(a)) / (1.0 - math.sin(a)))

# Freiburg im Breisgau
box = mapnik.Box2d(x_min, y_min, x_max, y_max)

This gives me the correct values and the correct box i want to have. Thanks to mkennedy with the tip of "Web Mercator".

1

I came up with this solution, found on http://bit.ly/1Qyd094

bounds = (-6.5, 49.5, 2.1, 59)

# spherical mercator (most common target map projection of osm data imported with osm2pgsql)
merc = mapnik.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')

# long/lat in degrees, aka ESPG:4326 and "WGS 84" 
longlat = mapnik.Projection('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')

bbox = mapnik.Envelope(bounds)

transform = mapnik.ProjTransform(longlat,merc)
merc_bbox = transform.forward(bbox)

m.zoom_to_box(merc_bbox)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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