13

我有墨尔本一个小区域的纬度/经度值；-37.803134,145.132377 以及我从openstreet 地图（Osmarender Image）导出的平面图像。图像宽度：1018 和高度：916

我希望能够使用 C++ 将 Lat/Long 转换为 X、Y 坐标，该点将反映该位置。

我使用了我在网上找到的各种公式，如下所示，但没有任何帮助。

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

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

4

4 回答 4

23

您需要更多信息，而不仅仅是一个纬度/经度对才能做到这一点。

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

  • 您的图像覆盖了多大的区域（以纬度/经度计）？根据您提供的信息，我不知道图像显示的区域是一米宽还是一公里宽。
  • 您的参考坐标（-37.803134、145.132377）指的是图像上的哪个点？它是角落之一吗？中间某处？

我还将假设您的图像是南北对齐的 - 例如，它没有指向左上角的北。这往往会使事情复杂化。

最简单的方法是准确计算出 (0, 0) 像素和 (1017, 915) 像素对应的纬度/经度坐标。然后您可以通过插值计算出与给定纬度/经度坐标对应的像素。

为了简要概述该过程，假设您的 (-37.803134, 145.132377) lat/lon 对应于您的 (0, 0) 像素，并且您发现您的 (1017, 915) 像素对应于 lat/lon (- 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 个像素。我会让你自己计算出水平等价物——基本上是一样的。注意-1带是因为如果MAP_HEIGHT从零开始，最大像素数是 915，而不是 916。

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

于 2011-02-10T04:52:53.377 回答
9

如果您正在寻找将大地测量 (lot,lan) 精确转换为定义的笛卡尔坐标（距离参考点的 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;
}
于 2011-02-14T05:36:23.993 回答
3

我不会太担心地球曲率。我以前没有使用过openstreetmap，但我只是快速浏览了一下，似乎他们使用了墨卡托投影。

这仅仅意味着他们已经为您将地球扁平化为一个矩形，使 X 与经度成正比，而 Y 几乎与纬度成正比。

所以你可以继续使用 Mac 的简单公式，你会非常准确。对于您正在处理的小地图，您的纬度将比一个像素的价值少得多。即使在维多利亚大小的地图上，您也只会得到 2-3% 的误差。

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

于 2011-02-19T06:21:40.737 回答
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;
}
于 2016-04-08T15:28:16.990 回答