How can I get latitude, longitude from x, y on a Mercator map (JPEG)?

Question

I have a Mercator projection map as a JPEG and I would like to know how to relate a given x, y coordinate to its latitude and longitude. I've looked at the Gudermannian function but I honestly don't understand how to take that function and apply it. Namely, what input is it expecting? The implementation I found (JavaScript) seems to take a range between -PI and PI, but what's the correlation between my y-value in pixels and that range?

Also, I found this function which takes a latitude and returns the tile for Google Maps, which also uses Mercator. It would seem that if I knew how to inverse this function, I'd be pretty close to having my answer.

/*<summary>Get the vertical tile number from a latitude
using Mercator projection formula</summary>*/

    private int getMercatorLatitude(double lati)
    {
        double maxlat = Math.PI;

        double lat = lati;

        if (lat > 90) lat = lat - 180;
        if (lat < -90) lat = lat + 180;

        // conversion degre=>radians
        double phi = Math.PI * lat / 180;

        double res;
        //double temp = Math.Tan(Math.PI / 4 - phi / 2);
        //res = Math.Log(temp);
        res = 0.5 * Math.Log((1 + Math.Sin(phi)) / (1 - Math.Sin(phi)));
        double maxTileY = Math.Pow(2, zoom);
        int result = (int)(((1 - res / maxlat) / 2) * (maxTileY));

        return (result);
    }

Answer 1

Here is some code for you... Let me know if you need more explanation.

    /// <summary>
    /// Calculates the Y-value (inverse Gudermannian function) for a latitude. 
    /// <para><see cref="http://en.wikipedia.org/wiki/Gudermannian_function"/></para>
    /// </summary>
    /// <param name="latitude">The latitude in degrees to use for calculating the Y-value.</param>
    /// <returns>The Y-value for the given latitude.</returns>
    public static double GudermannianInv(double latitude)
    {
        double sign = Math.Sign(latitude);
        double sin = Math.Sin(latitude * RADIANS_PER_DEGREE * sign);
        return sign * (Math.Log((1.0 + sin) / (1.0 - sin)) / 2.0);
    }

    /// <summary>
    /// Returns the Latitude in degrees for a given Y.
    /// </summary>
    /// <param name="y">Y is in the range of +PI to -PI.</param>
    /// <returns>Latitude in degrees.</returns>
    public static double Gudermannian(double y)
    {
        return Math.Atan(Math.Sinh(y)) * DEGREES_PER_RADIAN;
    }

Answer 2

Erich Mirabal's answer was completely correct (if not completely complete).

I have just tested it using a 'theoretical 256x256 Mercator tile' (Google's single tile version of a world map).

Here's a little more code (JavaScript, but easy to follow) to elucidate.

I live in Australia, at a latitude of about -33°.

convertRange(
    GudermannianInv(-33), 
    [Math.PI, - Math.PI], 
    [0, 256]
);

152.88327883810192

If you count 152 pixels down from the top of the tile, you will find Australia. I have also verified this answer is correct by comparing the result to known-good functions.

To be sure, we can reverse that calculation:

Gudermannian(
    convertRange(
        152.88, 
        [0, 256], 
        [Math.PI, - Math.PI]
));

And we are returned -32.99613291758226.

The tricky part isn't in the Gudermannian function, but in the conversion between two scales.

Fortunately, being rather lazy, and hating these kind of scaling problems, I already had a little function to do that messy conversion for me.

    /**
     * convert number from _n_ of r1[0] .. r1[1] to _n_ of r2[0] .. r2[1]
     * @example `convertRange( 5, [0, 10], [0, 100] ) === 50`
     *
     * @param {number} value
     * @param {array<number>} r1 old range
     * @param {array<number>} r2 new range
     * @returns {number} value adjusted for new range
     */
    function convertRange( value, r1, r2 ) {
        return ( value - r1[0] )
             * ( r2[1] - r2[0] )
             / ( r1[1] - r1[0] )
             +   r2[0];
    }

And the JavaScript versions of the original functions are naturally:

function Gudermannian(y) {
    return Math.atan(Math.sinh(y)) * (180 / Math.PI)
}

function GudermannianInv(latitude)
{
    var sign = Math.sign(latitude);
    var sin  = Math.sin(
                          latitude 
                        * (Math.PI / 180) 
                        * sign
    );
    return sign * (
        Math.log(
            (1 + sin) / (1 - sin)
        ) / 2
    );
}

Answer 3

Google, etc., use "spherical Mercator", the Mercator projection using a spherical Earth model rather than the slower and more complex elliptical equations.

The transformations are available as part of the OpenLayers code:

http://docs.openlayers.org/library/spherical_mercator.html

Answer 4

I've done something similiar. Especially if you have an image from a part of the world. A cropped map or just not a complete world map: https://stackoverflow.com/a/10401734/730823

Answer 5

An important note when performing an inverse is that there is no such thing as "the mercator map" as is the case with most other map projections. Each mercator map in existence is different depending on the input phi value. According to wikipedia google uses 85.051129, and other map providers use 85.05113. Therefore the input values to Gudermannian must be scaled based on e.g. GudermannianInv(85.05113).