views:

3066

answers:

3

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);
    }
A: 

pi would map to 0, and -pi would map to your image height. I believe the mapping is linear.

Totally not true. See Erich's answer.

Jimmy
-1. Not linear at all. The Gudermannian function is what determines the y/latitude of a point.
Erich Mirabal
sorry for being unclear, i meant the mapping between [-pi, pi] used as the y-input to the gudermannian function, and the [0, height] range in pixels :>
Jimmy
I undid the downvote, but I still don't think your comment is quite right. The mapping is not linear - it is hyperbolic. It looks linear'ish for the majority of the function, but grows quickly near the poles.
Erich Mirabal
then don't undo the downvote! :>
Jimmy
Well. I was hoping you would fix your statement :)
Erich Mirabal
well, now I realize that my answer is totally off mark, so it's better off as a 'known wrong answer' than a rehash of your answer :P
Jimmy
+3  A: 

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"/&gt;&lt;/para&gt;
    /// </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;
    }
Erich Mirabal
Thank you! Trying this out right now.
Matthew Wensing
So, if I have a mercator map image 1588 pixels tall, and I want to know the latitude where y = 677, I would figure out 677 in terms of +PI to -PI and call Gudermannian(y_in_terms_of_pi)? I realize that's wrong, but you can see where I am stuck mentally here ...
Matthew Wensing
For example, on a 1588 pixel tall mercator map, 30.0N is 615 pixels from the top. But if I express 615 in terms of a linear range from PI (0) to -PI (1588), I get 615 -> 0.70824318. And calling the above Gudermannian(0.70824318) yields 37.5587, not 30.0.
Matthew Wensing
Obviously the problem is the 'expressing 615 in terms of a linear range'. So how do I do this?
Matthew Wensing
You basically have to do something like this: lat = Gudermannian(Ymax - ((y / Height) * (Ymax - Ymin))); where Ymax and Ymin are given by taking the inverse Gudermannian of +-85.05112878 (or whatever the bounds max and min latitudes of your image are) and Height is the size of your image. If you are tiling, this will also work as long as you replace the Ymax and Ymin with the tile's bounds and Height with the tile's size. Does that make sense?
Erich Mirabal
def convert_y_to_lat(y): Ymax = gudermannian_inv(gudermannian(math.pi)); Ymin = gudermannian_inv(gudermannian(math.pi)) * -1; Height = 1588; lat = gudermannian(Ymax - ((y / Height) * (Ymax - Ymin))); return lat;convert_y_to_lat(615) is returning 4.42964421829. I expected it to return 30.0. :-(
Matthew Wensing
You can't use PI. That is only the limit which you will never reach. Instead, use the value of 85.05112878 degrees (converted to radians).
Erich Mirabal
I get 37.5587 when y = 615. For 30 degrees, y would have to be 655.17. What makes you so sure that 615 == 30 degrees?
Erich Mirabal
A: 

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

Tim Sylvester