将Lat/Longs转换为X/Y坐标

Question

我在墨尔本有一个小区域的Lat/Long值; -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);

如果任何人给我清楚解释如何做到这一点将是非常有帮助的.任何代码将非常感激.

Answer 1

您需要的信息不仅仅是单个纬度/经度对,而是能够做到这一点.

在此阶段,您提供的信息缺少两件事:

• 你的图像覆盖的区域有多大(就纬度/经度而言)？根据您提供的内容,我不知道图像是显示宽度为1米还是宽度为1千米的区域.
• 您的图像上的哪个点是您的参考坐标(-37.803134,145.132377)？这是角落之一吗？在某个地方的中间？

我还假设你的图像是北/南对齐的 - 例如,它没有指向左上角的北方.这往往会使事情复杂化.

最简单的方法是确切地确定哪些lat/lon坐标对应于(0,0)像素和(1017,915)像素.然后你可以通过插值找出对应于给定纬度/经度坐标的像素.

要简要概述该过程,假设您的(-37.803134,145.132377)lat/lon对应于您的(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%.最后一位将其转换为等效像素.

HTH!

编辑好的,根据你的评论,我可以更具体一点.鉴于您似乎使用左上角作为主要参考点,我将使用以下定义:

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,他的回答提醒我,情况就是这样).

Answer 2

如果你正在寻找准确的大地测量(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;
}