如何从Mercator地图（JPEG）的x、y坐标获取纬度和经度？

Question

我有一张用墨卡托投影的JPEG地图，我想知道如何将给定的x、y坐标与其纬度和经度相关联。我查看了Gudermannian函数，但是我不知道如何使用它。具体来说，它需要哪些输入？我找到的实现（JavaScript）似乎接受-PI到PI之间的范围，但是我的像素y值和该范围之间有什么关系呢？

此外，我找到了一个函数，它接受纬度并返回Google Maps的瓦片，这也使用了墨卡托投影。如果我知道如何反转此函数，我就可以非常接近得到我的答案了。

/*<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);
    }

- Matthew Wensing

---

如果我没记错的话，Google 使用的是等距投影而不是墨卡托投影。 - Stefano Borini

Virtual Earth和Google都使用Mercator投影。 - Erich Mirabal

此外，在使用墨卡托投影时，最大有用纬度不是+-90度，而是大约+-85.05112878度。在极点处该值为无穷大，因此您必须将其限制并忽略极点。 - Erich Mirabal

为了真正解决这个问题，你还需要知道在处理瓦片时的缩放级别。 - Erich Mirabal

我有点不清楚。你想要做 (x, y) -> (lat, long) 还是 (lat, long) -> (x,y)？ - ntownsend

Answer 1

以下是一些代码... 如果需要更多解释，请告诉我。

    /// <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;
    }

- Erich Mirabal

---

那么，如果我有一张Mercator地图图片，高度为1588像素，并且我想知道当y = 677时的纬度，我应该如何计算？我需要将677转换为+PI到-PI之间的值，然后调用Gudermannian(y_in_terms_of_pi)函数吗？我意识到这是错误的，但你可以看出我在思维上卡住了... - Matthew Wensing

显然问题在于“用线性范围表示615”。那么我该怎么做呢？ - Matthew Wensing

你基本上需要做这样的事情：lat = Gudermannian(Ymax - ((y / Height) * (Ymax - Ymin))); 其中Ymax和Ymin是通过取+-85.05112878的反Gudermannian（或者您的图像最大和最小纬度边界）得到的，而Height是您的图像大小。如果您正在平铺，则只需将Ymax和Ymin替换为瓦片的边界，将Height替换为瓦片的大小即可。这有意义吗？ - 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) 返回的是4.42964421829。我期望它返回30.0. :-( - Matthew Wensing

在我的墨卡托世界地图瓦片上，重新缩放为1588高，像素615（从顶部）穿越西班牙的卡塔赫纳附近和西西里岛的卡塔尼亚附近，大约在北纬37度32分左右，因此我不认为30.0是正确的。 - Roger

Answer 2

Erich Mirabal的回答是完全正确的（如果不是完全完整的话）。

我刚刚使用了一个“理论上的256x256墨卡托瓦片”进行了测试（谷歌的单个瓦片版本的世界地图）。

以下是一些代码（JavaScript，但易于理解）以阐明。

我住在澳大利亚，纬度约为-33°。

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

152.88327883810192

如果你从瓷砖顶部向下数152个像素，你会找到澳大利亚。我还通过将结果与已知的好函数进行比较来验证此答案是正确的。

为了确保，我们可以反过来计算：

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

我们得到的结果是-32.99613291758226。

棘手的部分不在古德曼函数中，而在于两个尺度之间的转换。

幸运的是，由于我相当懒惰并且讨厌这些缩放问题，我已经有了一个小函数来帮我解决这个混乱的转换问题。

    /**
     * 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];
    }

原始函数的JavaScript版本自然是：

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

- Orwellophile

Answer 3

Google等使用“球形Mercator投影”，该投影使用球形地球模型而不是更慢且更复杂的椭圆方程式。这些转换可作为OpenLayers代码的一部分使用。

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

- Tim Sylvester

Answer 4

我做过类似的事情。特别是当你有来自世界某个地区的图像时。裁剪的地图或者不完整的世界地图：https://dev59.com/Z3I95IYBdhLWcg3w5SSh#10401734

- Raphael

Answer 5

进行反演时需要注意的一点是，与大多数其他地图投影不同，不存在所谓的“the mercator map”。每个存在的墨卡托地图都取决于输入的phi值。根据维基百科，Google使用85.051129，而其他地图提供商使用85.05113。因此，Gudermannian的输入值必须根据例如GudermannianInv(85.05113)进行缩放。

- J Tileson