Converting from Lat/Lon (WGS84) to Mercator (UTM)
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.
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
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”.
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)
- Database Administration Tutorials
- Programming Tutorials & IT News
- Linux & DevOps World
- Entertainment & General News
- Games & eSport