Conversión de Lat/Longs en Coordenadas X/Y

Question

Tengo el valor Lat/Long de un área pequeña en Melbourne; -37.803134,145.132377 y también una imagen plana de eso, que exporté del mapa de openstreet (imagen de Osmarender). Ancho de imagen: 1018 y Altura: 916

Me gustaría poder convertir, usando C++, la coordenada Lat/Long en una X, Y donde el punto reflejaría la ubicación.

Utilicé varias fórmulas que encontré en la web como la que figura a continuación, pero nada ayuda.

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

Sería de gran ayuda si alguien me da una explicación clara de cómo hacer esto. Cualquier código sería muy apreciado.

Answer 1

Necesita más información que un solo par de latitud/longitud para poder hacer esto.

En esta etapa, la información que ha proporcionado no se encuentra dos cosas:

• tamaño del área cubre su imagen (en términos de latitud/longitud)? En función de lo que ha proporcionado, no sé si la imagen muestra un área de un metro de ancho o un kilómetro de ancho.
• ¿en qué lugar de su imagen se refieren sus coordenadas de referencia (-37.803134, 145.132377)? ¿Es una de las esquinas? En el medio en algún lugar?

También voy a suponer que su imagen está alineada norte/sur; por ejemplo, no tiene el norte apuntando hacia la esquina superior izquierda. Eso tendería a complicar las cosas.

El enfoque más fácil es determinar exactamente qué coordenadas lat/lon corresponden al (0, 0) píxel y al (1017, 915) píxel. A continuación, puede calcular el píxel correspondiente a una coordenada lat/lon determinada a través del interpolation.

Para delinear brevemente ese proceso, imagine que su (-37.803134, 145.132377) latitud/lon corresponde a su (0, 0) píxel, y que ha descubierto que su (1017, 915) píxel corresponde al lat/lon (-37.798917, 145.138535). Asumiendo la convención usual con un píxel (0, 0) en la esquina inferior izquierda, esto significa que el norte está arriba en la imagen.

Entonces, si usted está interesado en el objetivo de coordenadas (-37.801465, 145.134984), se puede averiguar el número de píxeles correspondiente a la imagen que es el siguiente:

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

Es decir, los correspondientes píxel es 362 píxeles desde la parte inferior de la imagen. A continuación, puede hacer exactamente lo mismo para la colocación horizontal de píxeles, pero utilizando longitudes y X píxeles en su lugar.

La pieza ((targetLat - minLat)/(maxLat - minLat)) determina la distancia que se encuentra entre las dos coordenadas de referencia, y da 0 para indicar que está en la primera, 1 para indicar la segunda y los números intermedios para indicar las ubicaciones intermedias. Entonces, por ejemplo, produciría 0.25 para indicar que estás al 25% del camino al norte entre las dos coordenadas de referencia. El último bit lo convierte en los píxeles equivalentes.

HTH!

EDIT Bueno, según su comentario puedo ser un poco más específico.Teniendo en cuenta que parece estar utilizando la esquina superior izquierda como punto de referencia principal, voy a utilizar las siguientes definiciones:

minLat = -37.803134 
maxLat = -37.806232 
MAP_HEIGHT = 916 

Entonces, si utilizamos un ejemplo coordenada (-37.804465, 145.134984), la Y de coordenadas del píxel correspondiente, relativa a la esquina superior izquierda, es decir:

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

por lo tanto, el píxel correspondiente es de 393 píxeles hacia abajo desde la parte superior. Te dejaré calcular el equivalente horizontal para ti - es básicamente lo mismo. NOTA El -1 con el MAP_HEIGHT se debe a que si se inicia en cero, el número máximo de píxeles es 915, no 916.

EDIT: Una cosa que me gustaría tener la oportunidad de señalar es que este es una aproximación. En realidad, no hay una relación lineal simple entre las coordenadas de latitud y longitud y otras formas de coordenadas cartesianas por varias razones, incluidas las proyecciones que se utilizan al hacer mapas, y el hecho de que la Tierra no es una esfera perfecta. En áreas pequeñas, esta aproximación es lo suficientemente cercana como para no hacer una diferencia significativa, pero a mayor escala las discrepancias pueden ser evidentes. Dependiendo de tus necesidades, YMMV. (Mi agradecimiento a uray, cuya respuesta a continuación me recordó que este es el caso).

Answer 2

Si está buscando la conversión precisa de geodésico (lote, lan) en su coordenada cartesiana definida (x, y metros del punto de referencia), puede hacer con mi fragmento de código aquí, esta función aceptará coordenada geodésica en radián y la salida el resultado en x, y

entrada:

• refLat, refLon: de coordenadas geodésico que ha definido como 0,0 en cartesiano de coordenadas (la unidad está en radianes)
• lat , lon: geodésico coordenadas que desea calcular la coordenada cartesiana (la unidad está en radianes)
• xOffset, yOffset: el resultado en de coordenadas cartesianas x, y (la unidad está en metros)

código:

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

Answer 3

No me preocuparía demasiado sobre la curvatura de la tierra. No he usado openstreetmap anteriormente, pero acabo de echarle un vistazo rápido y parece que usan una proyección de Mercator.

Lo que simplemente significa que han aplanado el planeta sobre un rectángulo, haciendo que X sea proporcional a la longitud, y que Y sea casi exactamente proporcional a Latitude.

Así que puede seguir adelante y usar las fórmulas simples de Mac y será muy preciso. Su Latitude se reducirá en mucho menos que el valor de un píxel para el pequeño mapa con el que está tratando. Incluso en un mapa del tamaño de Victoria solo obtendría un error del 2-3%.

diverscuba23 señaló que debe elegir un elipsoide ... openstreetmap usa WGS84, y también la mayoría de los mapas modernos. Sin embargo, tenga en cuenta que muchos mapas en Australia usan el AGD66 anterior, que puede diferir en 100-200 metros más o menos.?

Answer 4

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