Translate extent to use with gdal_translate

gdal_transform_extent_to_epsg3857.py

# coding: utf-8

import subprocess
import sys
import getopt
import math

# Formulas from http://wiki.openstreetmap.org/wiki/Mercator#Python_Implementation
def merc_x(lon):
    r_major = 6378137.000
    return r_major * math.radians(lon)


def merc_y(lat):
    if lat > 89.5:
        lat = 89.5
    if lat < -89.5:
        lat = -89.5
    r_major = 6378137.000
    r_minor = 6356752.3142
    temp = r_minor / r_major
    eccent = math.sqrt(1 - temp**2)
    phi = math.radians(lat)
    sinphi = math.sin(phi)
    con = eccent * sinphi
    com = eccent / 2
    con = ((1.0 - con) / (1.0 + con))**com
    ts = math.tan((math.pi / 2 - phi) / 2) / con
    y = 0 - r_major * math.log(ts)
    return y


def usage():
    print """
    You should run the example with:
    python gdal_transform_extent_to_epsg3857.py "5.11 45.21 5.12 45.20"
    """


def main(argv):
    try:
        opts, args = getopt.getopt(argv, "hg:d")
        if len(args) == 0:
            usage()
            sys.exit(2)
        extent = args[0]
        # extend = "5.11 45.21 5.12 45.20"
        coords = [float(i) for i in extent.split(" ")]
        x1t = merc_x(coords[0])
        y1t = merc_y(coords[1])
        x2t = merc_x(coords[2])
        y2t = merc_y(coords[3])
        cmd = "gdal_translate -of png -a_srs epsg:3857 -projwin \
            %s %s %s %s osm.xml OSM.png" % (
            str(x1t),
            str(y1t),
            str(x2t),
            str(y2t)
        )

        print cmd
        subprocess.Popen(cmd, stdout=subprocess.PIPE, shell=True)

    except getopt.GetoptError:
        usage()
        sys.exit(2)

if __name__ == "__main__":
    main(sys.argv[1:])