Converting from Lat/Lon (WGS84) to Mercator (UTM)

Question

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.

Answer 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

Answer 2

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".

Answer 3

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)