将经纬度转换为X/Y坐标

Question

将经纬度转换为X/Y坐标

13

我拥有墨尔本一个小区域的纬度/经度值：-37.803134,145.132377，并且还有一张该区域的平面图像，是我从OpenStreetMap（Osmarender Image）导出的。图片宽度：1018，高度：916。

我想使用C++将纬度/经度转换为X、Y坐标，使得该点反映该位置。

我尝试了一些在网上找到的公式，如下所示，但没有任何帮助。

var y = ((-1 * lat) + 90) * (MAP_HEIGHT / 180);
var x = (lon + 180) * (MAP_WIDTH / 360);

如果有人能给我清晰的解释如何做到这一点，那将非常有帮助。任何代码都将不胜感激。

- Verve Innovation

4个回答

9

如果您需要将大地坐标（经度、纬度）精确转换为您定义的笛卡尔坐标系（距离参考点x、y米），您可以使用我的代码片段。此函数将接受以弧度表示的大地坐标，并输出x、y的结果。

输入：

- refLat、refLon：您在笛卡尔坐标系中定义为0,0的大地坐标（单位为弧度） - lat、lon：您要计算其笛卡尔坐标的大地坐标（单位为弧度） - xOffset、yOffset：笛卡尔坐标x、y的结果（单位为米）

代码：

#define GD_semiMajorAxis 6378137.000000
#define GD_TranMercB     6356752.314245
#define GD_geocentF      0.003352810664

void geodeticOffsetInv( double refLat, double refLon,
                        double lat,    double lon, 
                        double& xOffset, double& yOffset )
{
    double a = GD_semiMajorAxis;
    double b = GD_TranMercB;
    double f = GD_geocentF;

    double L     = lon-refLon
    double U1    = atan((1-f) * tan(refLat));
    double U2    = atan((1-f) * tan(lat));
    double sinU1 = sin(U1); 
    double cosU1 = cos(U1);
    double sinU2 = sin(U2);
    double cosU2 = cos(U2);

    double lambda = L;
    double lambdaP;
    double sinSigma;
    double sigma;
    double cosSigma;
    double cosSqAlpha;
    double cos2SigmaM;
    double sinLambda;
    double cosLambda;
    double sinAlpha;
    int iterLimit = 100;
    do {
        sinLambda = sin(lambda);
        cosLambda = cos(lambda);
        sinSigma = sqrt((cosU2*sinLambda) * (cosU2*sinLambda) + 
                        (cosU1*sinU2-sinU1*cosU2*cosLambda) * 
                        (cosU1*sinU2-sinU1*cosU2*cosLambda) );
        if (sinSigma==0)
        {
            xOffset = 0.0;
            yOffset = 0.0;
            return ;  // co-incident points
        }
        cosSigma    = sinU1*sinU2 + cosU1*cosU2*cosLambda;
        sigma       = atan2(sinSigma, cosSigma);
        sinAlpha    = cosU1 * cosU2 * sinLambda / sinSigma;
        cosSqAlpha  = 1 - sinAlpha*sinAlpha;
        cos2SigmaM  = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
        if (cos2SigmaM != cos2SigmaM) //isNaN
        {
            cos2SigmaM = 0;  // equatorial line: cosSqAlpha=0 (§6)
        }
        double C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
        lambdaP = lambda;
        lambda = L + (1-C) * f * sinAlpha *
            (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
    } while (fabs(lambda-lambdaP) > 1e-12 && --iterLimit>0);

    if (iterLimit==0)
    {
        xOffset = 0.0;
        yOffset = 0.0;
        return;  // formula failed to converge
    }

    double uSq  = cosSqAlpha * (a*a - b*b) / (b*b);
    double A    = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
    double B    = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
    double deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
        B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
    double s = b*A*(sigma-deltaSigma);

    double bearing = atan2(cosU2*sinLambda,  cosU1*sinU2-sinU1*cosU2*cosLambda);
    xOffset = sin(bearing)*s;
    yOffset = cos(bearing)*s;
}

- uray

嘿，感谢您的输入。但这是什么公式？它充满了数学，很难理解。如果您告诉我如何完成它或提供一些参考链接，那将更容易让我理解。 - Verve Innovation

@ITion：与其将纬度/经度坐标之间的差异转换为像素位置，它将其改变为米数差异。如果我理解正确，这更加复杂，因为它考虑到了行星不是完美球体的事实——我认为在您所处理的小区域内这可能并不重要。 - Mac

@Mac：没错，对于ITion正在处理的领域，你可以简化事情，忘记地球实际上是一个椭球体，但是一旦你开始处理更大的区域，你就不能再忽略这个事实了...然后你就得想出你要相信地球是哪个椭球体...这是另一个完全不同的讨论...但其中一个更流行的是WGS-84。 - diverscuba23

3

我不太担心地球曲率的问题。

我以前没有使用过OpenStreetMap，但我刚刚查看了一下，它似乎使用了墨卡托投影。

这意味着他们已经将地球压缩成一个矩形，使得X与经度成比例，而Y几乎完全与纬度成比例。

因此，您可以使用Mac的简单公式，非常准确。您的纬度误差将小于处理的小地图的一个像素值。即使在维多利亚州这样大小的地图上，您也只会出现2-3％的误差。

diverscuba23指出，您必须选择一个椭球体... OpenStreetMap使用WGS84，大多数现代映射也是如此。然而，请注意，澳大利亚的许多地图使用较旧的AGD66，可能相差100-200米左右。

- Jake R

2

double testClass::getX(double lon, int width)
{
    // width is map width
    double x = fmod((width*(180+lon)/360), (width +(width/2)));

    return x;
}

double testClass::getY(double lat, int height, int width)
{
    // height and width are map height and width
    double PI = 3.14159265359;
    double latRad = lat*PI/180;

    // get y value
    double mercN = log(tan((PI/4)+(latRad/2)));
    double y     = (height/2)-(width*mercN/(2*PI));
    return y;
}

- AngryDuck

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Mac · Accepted Answer

为了完成此任务，您需要比一个纬度/经度对更多的信息。

在这个阶段，您提供的信息缺少两个方面：

您的图像覆盖的区域有多大（以纬度/经度表示）？根据您提供的信息，我不知道图像是否显示了1米宽或1公里宽的区域。
您参考坐标（-37.803134，145.132377）指的是图像上的哪个位置？是其中一个角落吗？还是中间的某个地方？

我还假设您的图像是朝南/北对齐的 - 例如，北方不会指向左上角。那会使事情变得复杂。

最简单的方法是确定（0,0）像素和（1017,915）像素分别对应于哪些纬度/经度坐标。然后，您可以通过插值确定与给定纬度/经度坐标相对应的像素。

简要概述该过程，假设您的（-37.803134，145.132377）纬度/经度对应于（0,0）像素，并且您已经发现您的（1017,915）像素对应于纬度/经度为（-37.798917，145.138535）的位置。假设像素（0,0）在左下角，这意味着图像向上是北方。

然后，如果您对目标坐标（-37.801465，145.134984）感兴趣，则可以通过以下方式确定相应的像素数量：

pixelY = ((targetLat - minLat) / (maxLat - minLat)) * (maxYPixel - minYPixel)
       = ((-37.801465 - -37.803134) / (-37.798917 - -37.803134)) * (915 - 0)
       = 362.138

也就是说，对应的像素距离图像底部362个像素。然后你可以按照同样的方式进行水平像素定位，但使用经度和X像素。

部分((targetLat - minLat) / (maxLat - minLat))计算出您在两个参考坐标之间的距离，并给出0表示您在第一个坐标处，1表示第二个坐标处，以及介于两者之间的数字表示其它位置。例如，它将产生0.25，表示您在两个参考坐标之间向北移动了25%。最后一部分将其转换为相应的像素。

希望有所帮助！

编辑好的，根据您的评论，我可以更具体地说明。假设您将左上角作为主要参考点，我将使用以下定义：

minLat = -37.803134
maxLat = -37.806232
MAP_HEIGHT = 916

然后，如果我们使用坐标示例(-37.804465, 145.134984)，相应像素的Y坐标，相对于左上角，是:

pixelY = ((targetLat - minLat) / (maxLat - minLat)) * (MAP_HEIGHT - 1)
       = ((-37.804465 - -37.803134) / (-37.806232 - -37.803134)) * 915
       = 393.11

因此，相应的像素距离顶部下方393个像素。我会让你自己计算水平等价物-基本上是相同的。注意在MAP_HEIGHT中使用-1的原因是，如果从零开始，则最大像素数为915，而不是916。

编辑：我想借此机会指出的一件事是，这是一种近似方法。实际上，纬度和经度坐标与其他形式的笛卡尔坐标之间没有简单的线性关系，原因包括制作地图时使用的投影方式以及地球不是一个完美的球体。在小范围内，此近似足够接近，没有显着的差异，但在较大的尺度上，可能会出现差异。根据您的需求，结果可能有所不同。（感谢uray，他下面的答案提醒我这种情况）。