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

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.

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

En esta etapa, la información que ha proporcionado le falta dos cosas:

  • ¿Cuán grande es el área que cubre la imagen (en términos de lat / lon)? 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é punto 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 averiguar exactamente qué coordenadas lat / lon corresponden al (0, 0) píxel y al (1017, 915) píxel. Luego puede calcular el píxel correspondiente a una determinada coordenada de latitud / longitud mediante la interpolación .

Para delinear brevemente ese proceso, imagina que tu (-37.803134, 145.132377) latitud / lon corresponde a tu (0, 0) píxel, y que has descubierto que tu (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.

Luego, si está interesado en la coordenada del objective (-37.801465, 145.134984), puede determinar el número correspondiente de píxeles en la imagen de la siguiente manera:

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

Es decir, el píxel correspondiente 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 parte ((targetLat - minLat) / (maxLat - minLat)) determina cuán lejos 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 en entre para indicar 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!

EDITAR Ok, basándome en tu comentario, puedo ser un poco más específico. Dado que parece que está utilizando la esquina superior izquierda como su punto de referencia principal, usaré las siguientes definiciones:

 minLat = -37.803134 maxLat = -37.806232 MAP_HEIGHT = 916 

Luego, si usamos una coordenada de ejemplo de (-37.804465, 145.134984), la coordenada Y del píxel correspondiente, relativa a la esquina superior izquierda, es:

 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 está 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 MAP_HEIGHT es porque si comienza en cero, el número máximo de píxeles es 915, no 916.

EDITAR: Una cosa que me gustaría aprovechar para señalar es que esta 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).

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

entrada:

  • refLat, refLon: coordenada geodésica que ha definido como 0,0 en coordenadas cartesianas (la unidad está en radianes)
  • lat, lon: coordenada geodésica para la que desea calcular su coordenada cartesiana (la unidad está en radianes)
  • xOffset, yOffset: el resultado en 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; } 

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 hacia un rectángulo, haciendo que X sea proporcional a la Longitud, y que Y sea casi exactamente proporcional a Latitude.

Así que puedes seguir adelante y usar las fórmulas simples de Mac y serás 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.

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