Конвертировать широты / долготы в координаты X / Y - PullRequest
13 голосов
/ 10 февраля 2011

У меня есть значение широты / долготы маленькой области в Мельбурне; -37.803134,145.132377, а также плоское изображение того, что я экспортировал из openstreet map (Osmarender Image). Ширина изображения: 1018 и высота: 916

Я хотел бы иметь возможность преобразовать, используя C ++, широту / долготу в координаты X, Y, где точка будет отражать местоположение.

Я использовал различные формулы, которые я нашел в сети, как показано ниже, но ничего не помогает.

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

Было бы очень полезно, если бы кто-нибудь дал мне четкое объяснение того, как это сделать. Любой код будет высоко ценится.

Ответы [ 4 ]

24 голосов
/ 10 февраля 2011

Вам нужно больше информации, чем просто одна широта / долгота, чтобы сделать это.

На данном этапе в предоставленной вами информации отсутствуют две вещи:

  • Какую площадь покрывает ваше изображение (в широте / долготе)? Исходя из того, что вы предоставили, я не знаю, показывает ли изображение область шириной метр или километр.
  • к какой точке вашего изображения относится ваша опорная координата (-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 пикселя. Я позволю вам разработать горизонтальный эквивалент для себя - это в основном то же самое. ПРИМЕЧАНИЕ * * -1 с MAP_HEIGHT потому, что если вы начинаете с нуля, максимальное число пикселей будет 915, а не 916.

РЕДАКТИРОВАТЬ: Я хотел бы воспользоваться возможностью, чтобы указать, что это приблизительное. В действительности не существует простой линейной связи между координатами широты и долготы и другими формами декартовых координат по ряду причин, включая проекции, которые используются при создании карт, и тот факт, что Земля не является идеальной сферой. В небольших областях это приближение достаточно близко, чтобы не иметь существенных различий, но в более крупных масштабах расхождения могут стать очевидными. В зависимости от ваших потребностей, YMMV. (Мое спасибо Урайу, чей ответ ниже напомнил мне, что это так).

8 голосов
/ 14 февраля 2011

если вы ищете для точного преобразования геодезических (много, LAN) в ваш определенном декартовые координаты (х, у м от контрольной точки), вы можете сделать с моей фрагмент кода здесь, эта функция будет принимать геодезический координат в радианахи вывести результат в виде x, y

input:

  • 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;
}
3 голосов
/ 19 февраля 2011

Я бы не слишком беспокоился о искривлении земли. Я раньше не использовал openstreetmap, но я только что посмотрел, и оказалось, что они используют проекцию Меркатора.

Это просто означает, что они выровняли планету для вас в прямоугольник, в результате чего X пропорционально долготе, а Y почти точно пропорционально широте.

Так что вы можете пойти дальше и использовать простые формулы Mac, и вы будете очень точны. Ваша широта будет меньше, чем значение пикселя для маленькой карты, с которой вы имеете дело. Даже на карте размером с Викторию вы получите ошибку только в 2-3%.

diverscuba23 указал, что вы должны выбрать эллипсоид ... openstreetmap использует WGS84, как и большинство современных карт. Однако следует помнить, что на многих картах в Австралии используется более старый AGD66, который может отличаться на 100-200 метров или около того.

2 голосов
/ 08 апреля 2016
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;
}
...