Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Render a static map using Mapnik given centre, zoom and image width/height.
#!/usr/bin/env python
# This script is like generate_image.py from the OSM Mapnik code,
# but it renders based on a given centre point, zoom and final image pixel size
# Author: Andrew Harvey <andrew.harvey4@gmail.com>
# License: CC0 http://creativecommons.org/publicdomain/zero/1.0/
#
# To the extent possible under law, the person who associated CC0
# with this work has waived all copyright and related or neighboring
# rights to this work.
# http://creativecommons.org/publicdomain/zero/1.0/
# For testing the tile 17/120590/78656 should match the output of
# ./render_map_preview.py -x 151.21170043945312 -y -33.871555787081476 -z 17 -w 256 -h 256
# that tile has centre (-33.871555787081476, 151.21170043945312) -> 16832809.4949049 -4011568.11846262
# and bounds (-33.87041555094183, 151.2103271484375),(-33.87269600798948, 151.21307373046875)
# which are projected as (gdaltransform -s_srs '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs' -t_srs '+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')
# xmin=16832656.6208484
# ymin=-4011415.24440605
# xmax=16832962.3689615
# ymax=-4011720.99251919
# dx=306 ie. so 306 in projected units for this one tile
try:
import mapnik2 as mapnik
except:
import mapnik
import getopt, sys, os
# target projection
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')
# WGS lat/long source projection of centre
longlat = mapnik.Projection('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
# ensure minimum Mapnik version
if not hasattr(mapnik,'mapnik_version') and not mapnik.mapnik_version() >= 600:
raise SystemExit('This script requires Mapnik >=0.6.0)')
if __name__ == "__main__":
# get arguments
try:
opts, args = getopt.getopt(sys.argv[1:], "m:x:y:z:w:h:o:", ["mapfile=", "centrex=", "centrey=", "zoom=", "width=", "height=", "output="])
except getopt.GetoptError, err:
print str(err)
#usage()
sys.exit(2)
for o, a in opts:
if o in ("-m", "--mapfile"):
mapfile = a
elif o in ("-x", "--centrex"):
centrex = float(a)
elif o in ("-y", "--centrey"):
centrey = float(a)
elif o in ("-z", "--zoom"):
zoom = int(a)
elif o in ("-w", "--width"):
width = int(a)
elif o in ("-h", "--height"):
height = int(a)
elif o in ("-o", "--output"):
output = a
# make a new Map object for the given mapfile
m = mapnik.Map(width, height)
mapnik.load_map(m, mapfile)
# ensure the target map projection is mercator
m.srs = merc.params()
# transform the centre point into the target coord sys
centre = mapnik.Coord(centrex, centrey)
transform = mapnik.ProjTransform(longlat, merc)
merc_centre = transform.forward(centre)
# 360/(2**zoom) degrees = 256 px
# so in merc 1px = (20037508.34*2) / (256 * 2**zoom)
# hence to find the bounds of our rectangle in projected coordinates + and - half the image width worth of projected coord units
dx = ((20037508.34*2*(width/2)))/(256*(2 ** (zoom)))
minx = merc_centre.x - dx
maxx = merc_centre.x + dx
# grow the height bbox, as we only accurately set the width bbox
m.aspect_fix_mode = mapnik.aspect_fix_mode.ADJUST_BBOX_HEIGHT
bounds = mapnik.Box2d(minx, merc_centre.y-10, maxx, merc_centre.y+10) # the y bounds will be fixed by mapnik due to ADJUST_BBOX_HEIGHT
m.zoom_to_box(bounds)
# render the map image to a file
mapnik.render_to_file(m, output)
sys.stdout.write('output image to %s!\n' % output)
@anthonyprior

This comment has been minimized.

Copy link

@anthonyprior anthonyprior commented Dec 8, 2015

A thousand thankyous this is exactly what I needed and was really struggling with!

@thisisbinod

This comment has been minimized.

Copy link

@thisisbinod thisisbinod commented Jan 29, 2016

Thanks a lot!

@cbdavide

This comment has been minimized.

Copy link

@cbdavide cbdavide commented Jul 7, 2017

Sorry for bother, but I'm trying to find a way to translate from the pixel coordinates of the map to lon/lat coordinates of the map, but I'm having a rough time finding some useful information, do you know if there is a way to accomplish this ?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment