#!/usr/bin/python2

import sys
import re
from math import *
from mapnik import *

webwidth, webheight = 1600, 1200
outwidth, outheight = 0, 0
outimage = 'renderosm.png' 
stylesheet = '/data/openstreetmap/mapnik-style-old-svn/osm.xml'

proj= Projection("+init=epsg:3857")
canvassize= 2 * proj.forward(Coord(180, 0)).x
tilesize= 256

def usageexit():
    print("""
Usage:
renderosm.py [options] <OpenStreetMap URL>
renderosm.py [options] <lat1,long1,lat2,long2>

Options:
-o <output.png>     Output file
-w <width>          Width of output image in pixels (default 3000)
-h <height>         Height of output image in pixels

The zoom level on which the map appearance depends will be computed by the
Mapnik library from the output image size.  The zoom level from the
OpenStreetMap URL is used only for determining the map region.  The map
rectangle represented by an OpenStreetMap URL is assumed to be the equivalent
of a 1600 by 1200 browser window or the same area with a different aspect ratio
if both -w and -h are present.
""")
    exit()

def latlong2tile(lat, long, zoom):
    tilecoord= proj.forward(Coord(long, lat))
    tilex= 2 ** (zoom - 1) + tilecoord.x * 2 ** zoom / canvassize
    tiley= 2 ** (zoom - 1) - 1 - tilecoord.y * 2 ** zoom / canvassize
    return tilex, tiley

def osmurl2box(lat, long, zoom, pixwidth, pixheight):
    # half window size in projection coordinates:
    halfwidth= canvassize / 2 ** zoom / tilesize * pixwidth / 2
    halfheight= canvassize / 2 ** zoom / tilesize * pixheight / 2
    centre= proj.forward(Coord(long, lat))
    return Box2d(centre.x - halfwidth, centre.y - halfheight, centre.x + halfwidth, centre.y + halfheight)

if len(sys.argv) == 1 or sys.argv[1] == "-h" or sys.argv[1] == "--help":
    usageexit()

while len(sys.argv) > 1 and sys.argv[1][0] == '-':
    if len(sys.argv) == 2:
        print("Missing mandatory argument of " + sys.argv[1] + " option.")
        usageexit()
    if sys.argv[1] == "-o":
        outimage= sys.argv[2]
    elif sys.argv[1] == "-w":
        outwidth= int(sys.argv[2])
    elif sys.argv[1] == "-h":
        outheight= int(sys.argv[2])
    sys.argv[1:]= sys.argv[3:]

if len(sys.argv) == 1:
    print("Missing map region argument.")
    usageexit()
elif len(sys.argv) > 2:
    print("Too many arguments.")
    usageexit()

osmurlre= re.compile(r"^https?://(?:www\.)?(?:openstreetmap|osm)\.org/.*#map=(\d+)/([0-9.]+)/([0-9.]+)$")
georectre= re.compile(r"^([0-9.]+),([0-9.]+),([0-9.]+),([0-9.]+)$")
urlmatch= osmurlre.match(sys.argv[1])
geomatch= georectre.match(sys.argv[1])

if urlmatch:
    zoom, lat, long = urlmatch.groups()
    zoom, lat, long = int(zoom), float(lat), float(long)
    if outwidth != 0 and outheight != 0:
        area= webwidth * webheight
        webwidth= int(sqrt(area * outwidth / outheight))
    box= osmurl2box(lat, long, zoom, webwidth, webheight)
    aspect= float(webwidth) / webheight

elif geomatch:
    lat1, long1, lat2, long2 = geomatch.groups()
    corner1= proj.forward(Coord(float(long1), float(lat1)))
    corner2= proj.forward(Coord(float(long2), float(lat2)))
    box= Box2d(corner1.x, corner1.y, corner2.x, corner2.y)
    aspect= float(abs(corner1.x - corner2.x)) / abs(corner1.y - corner2.y)
    if outwidth != 0 and outheight != 0 and outheight != int(outwidth / aspect):
        print("Warning: output aspect ratio does not match aspect of rectangle exactly.  Adapting height.")
        outheight= int(outwidth / aspect)

else:
    print("Cannot understand map rectangle argument.")
    usageexit()


if outwidth == 0 and outheight == 0:
    outwidth= 3000
    outheight= int(outwidth / aspect)
elif outwidth == 0:
    outwidth= int(outheight * aspect)
elif outheight == 0:
    outheight= int(outwidth / aspect)

map= Map(outwidth, outheight, "+init=epsg:3857")
load_map(map, stylesheet)
map.zoom_to_box(box)
render_to_file(map, outimage)