El cálculo del punto medio lat / largo de Python da un resultado incorrecto cuando la longitud> 90

Tengo un problema con una función corta para calcular el punto medio de una línea cuando se dan la latitud y la longitud de los puntos en cada extremo. En pocas palabras, funciona correctamente cuando la longitud es mayor que -90 grados o menor que 90 grados. Para la otra mitad del planeta, proporciona un resultado algo aleatorio.

El código es una conversión python de javascript que se proporciona en http://www.movable-type.co.uk/scripts/latlong.html , y parece ajustarse a las versiones corregidas aquí y aquí . Al comparar con las dos versiones de stackoverflow, admito que no codifico en C # o Java, pero no puedo detectar dónde está mi error.

El código es el siguiente:

#!/usr/bin/python import math def midpoint(p1, p2): lat1, lat2 = math.radians(p1[0]), math.radians(p2[0]) lon1, lon2 = math.radians(p1[1]), math.radians(p2[1]) dlon = lon2 - lon1 dx = math.cos(lat2) * math.cos(dlon) dy = math.cos(lat2) * math.sin(dlon) lat3 = math.atan2(math.sin(lat1) + math.sin(lat2), math.sqrt((math.cos(lat1) + dx) * (math.cos(lat1) + dx) + dy * dy)) lon3 = lon1 + math.atan2(dy, math.cos(lat1) + dx) return(math.degrees(lat3), math.degrees(lon3)) p1 = (6.4, 45) p2 = (7.3, 43.5) print "Correct:", midpoint(p1, p2) p1 = (95.5,41.4) p2 = (96.3,41.8) print "Wrong:", midpoint(p1, p2) 

¿Alguna sugerencia?

Reemplace su código de configuración arg por:

 lat1, lon1 = p1 lat2, lon2 = p2 assert -90 <= lat1 <= 90 assert -90 <= lat2 <= 90 assert -180 <= lon1 <= 180 assert -180 <= lon2 <= 180 lat1, lon1, lat2, lon2 = map(math.radians, (lat1, lon1, lat2, lon2)) 

y ejecute su código de nuevo.

Actualice algunas sugerencias generales esperanzadamente útiles sobre cálculos con latitud / longitud:

  1. ¿Entrada lat / lon en grados o radianes?
  2. Compruebe la entrada lat / lon para el rango válido
  3. Marque OUTPUT lat / lon para el rango válido. La longitud tiene una discontinuidad en la línea de datos internacional.

La última parte de la rutina del punto medio podría cambiarse de manera útil para evitar un problema potencial con el uso a larga distancia:

 lon3 = lon1 + math.atan2(dy, math.cos(lat1) + dx) # replacement code follows: lon3d = math.degrees(lon3) if lon3d < -180: print "oops1", lon3d lon3d += 360 elif lon3d > 180: print "oops2", lon3d lon3d -= 360 return(math.degrees(lat3), lon3d) 

Por ejemplo, encontrar un punto medio entre Auckland, Nueva Zelanda (-36.9, 174.8) y Papeete, Tahití (-17.5, -149.5) produce oops2 194.270430902 en el camino a una respuesta válida (-28.355951246746923, -165.72956909809082)

En primer lugar, disculpas por dejarte otra respuesta. Tengo una solución al problema mencionado en esa respuesta que implica cómo encontrar el punto medio de dos puntos a cada lado de la línea de datos. Me encantaría simplemente agregar un comentario a la respuesta existente, pero no tengo la reputación de hacerlo.

La solución se encontró mirando los archivos de Javascript que activan la herramienta en http://www.movable-type.co.uk/scripts/latlong.html . Estoy usando shapely.geometry.Point. Si no desea instalar este paquete, usar tuplas en su lugar funcionaría igual.

 def midpoint(pointA, pointB): lonA = math.radians(pointA.x) lonB = math.radians(pointB.x) latA = math.radians(pointA.y) latB = math.radians(pointB.y) dLon = lonB - lonA Bx = math.cos(latB) * math.cos(dLon) By = math.cos(latB) * math.sin(dLon) latC = math.atan2(math.sin(latA) + math.sin(latB), math.sqrt((math.cos(latA) + Bx) * (math.cos(latA) + Bx) + By * By)) lonC = lonA + math.atan2(By, math.cos(latA) + Bx) lonC = (lonC + 3 * math.pi) % (2 * math.pi) - math.pi return Point(math.degrees(lonC), math.degrees(latC)) 

Espero que esto sea útil y que no se considere inapropiado, ya que es una respuesta a la pregunta planteada en la respuesta anterior.

No es una respuesta directa a la pregunta> 90, pero me ayudó en mi caso, así que la dejo aquí en caso de que pueda ayudar a otros.

Solo necesitaba colocar el punto medio en un mapa con lat, largo en grados. No utilicé conversión a radianes y se representa correctamente en el mapa utilizando una distancia euclidiana simple. Función de ejemplo:

 def midpoint_euclidean(x1,y1,x2,y2): dist_x = abs(x1-x2) / 2. dist_y = abs(y1-y2) / 2. res_x = x1 - dist_x if x1 > x2 else x2 - dist_x res_y = y1 - dist_y if y1 > y2 else y2 - dist_y return res_x, res_y