¿Diferencia entre Vincenty y los cálculos de distancia de gran círculo?

El paquete de geopy de Python presenta dos técnicas de medición de distancia: las fórmulas Great Circle y Vincenty .

>>> from geopy.distance import great_circle
>>> from geopy.distance import vincenty
>>> p1 = (31.8300167,35.0662833) # (lat, lon) - https://goo.gl/maps/TQwDd
>>> p2 = (31.8300000,35.0708167) # (lat, lon) - https://goo.gl/maps/lHrrg
>>> vincenty(p1, p2).meters
429.16765838976664
>>> great_circle(p3, p4).meters
428.4088367903001

¿Cuál es la diferencia? ¿Qué medida de distancia se prefiere?

Según Wikipedia, la fórmula de Vincenty es más lenta pero más precisa :

Las fórmulas de Vincenty son dos métodos iterativos relacionados utilizados en la geodesia para calcular la distancia entre dos puntos en la superficie de un esferoide, desarrollado por Thaddeus Vincenty (1975a) Se basan en el supuesto de que la figura de la Tierra es un esferoide achatado, y por lo tanto son más precisos que métodos como la distancia de gran círculo que supone una Tierra esférica.

La diferencia de precisión está ~0.17%en una distancia de 428 metros en Israel. He hecho una prueba de velocidad rápida y sucia:

<class 'geopy.distance.vincenty'>       : Total 0:00:04.125913, (0:00:00.000041 per calculation)
<class 'geopy.distance.great_circle'>   : Total 0:00:02.467479, (0:00:00.000024 per calculation)

Código:

import datetime
from geopy.distance import great_circle
from geopy.distance import vincenty
p1 = (31.8300167,35.0662833)
p2 = (31.83,35.0708167)

NUM_TESTS = 100000
for strategy in vincenty, great_circle:
    before = datetime.datetime.now()
    for i in range(NUM_TESTS):
        d=strategy(p1, p2).meters
    after = datetime.datetime.now()
    duration = after-before
    print "%-40s: Total %s, (%s per calculation)" % (strategy, duration, duration/NUM_TESTS)

Para concluir: la fórmula de Vincenty duplica el tiempo de cálculo en comparación con el gran círculo, y su ganancia de precisión en el punto probado es ~ 0.17%.

Dado que el tiempo de cálculo es insignificante, se prefiere la fórmula de Vincenty para cada necesidad práctica.

Actualización : Siguiendo los perspicaces comentarios de whuber y cffk y la respuesta de cffk , estoy de acuerdo en que la ganancia de precisión debe compararse con el error, no con la medición. Por lo tanto, la fórmula de Vincenty es unos pocos órdenes de magnitud más precisa, no ~ 0.17%.

Si está utilizando geopy, las distancias great_circle y vincenty son igualmente convenientes de obtener. En este caso, casi siempre debe usar el que le da el resultado más preciso, es decir, vincenty. Las dos consideraciones (como usted señala) son la velocidad y la precisión.

Vincenty es dos veces más lento. Pero probablemente en una aplicación real el aumento del tiempo de ejecución es insignificante. Incluso si su aplicación requería un millón de cálculos de distancia, solo estamos hablando de una diferencia en tiempos de un par de segundos.

Para los puntos que usa, el error en vincenty es de 6 μm y el error en la distancia del gran círculo es de 0.75 m. Entonces diría que vincenty es 120000 veces más preciso (en lugar de 0,17% más preciso). Para puntos generales, el error en la distancia del gran círculo puede ser de hasta 0.5%. Entonces, ¿puedes vivir con un error de 0.5% en las distancias? Para uso casual (¿cuál es la distancia de Ciudad del Cabo a El Cairo?), Probablemente puedas. Sin embargo, muchas aplicaciones SIG tienen requisitos de precisión mucho más estrictos. (0.5% son 5m en 1km. Eso realmente hace la diferencia).

Casi todo el trabajo de mapeo serio se lleva a cabo en el elipsoide de referencia y, por lo tanto, tiene sentido que las distancias también se deben medir en el elipsoide. Tal vez pueda escapar con distancias de gran círculo hoy. Pero para cada nueva aplicación, deberá verificar si aún es aceptable. Lo mejor es usar la distancia elipsoidal desde el principio. Dormirás mejor por la noche.

ADENDA (mayo de 2017)

En respuesta a la respuesta dada por @ craig-hicks. El método vincenty () en geopy tiene un defecto potencialmente fatal: arroja un error para puntos casi antipodales. La documentación en el código sugiere aumentar el número de iteraciones. Pero esta no es una solución general porque el método iterativo utilizado por vincenty () es inestable para tales puntos (cada iteración lo lleva más lejos de la solución correcta).

¿Por qué caracterizo el problema como "potencialmente fatal"? Debido a que cualquier uso de la función de distancia dentro de otra biblioteca de software debe ser capaz de manejar la excepción. Manejarlo devolviendo un NaN o la distancia del gran círculo puede no ser satisfactorio, porque la función de distancia resultante no obedecerá a la desigualdad del triángulo que impide su uso, por ejemplo, en árboles de puntos estratégicos.

La situación no es completamente sombría. Mi paquete python GeographicLib calcula la distancia geodésica con precisión y sin ningún fallo. La solicitud de extracción de geopy # 144 cambia la función de distancia de geopy para usar el paquete Geographiclib si está disponible. Lamentablemente, esta solicitud de extracción ha estado en el limbo desde Augest 2016.

ADENDA (mayo de 2018)

geopy 1.13.0 ahora usa el paquete Geographiclib para calcular distancias. Aquí hay una muestra de llamada (basada en el ejemplo de la pregunta original):

>>> from geopy.distance import great_circle
>>> from geopy.distance import geodesic
>>> p1 = (31.8300167,35.0662833) # (lat, lon) - https://goo.gl/maps/TQwDd
>>> p2 = (31.8300000,35.0708167) # (lat, lon) - https://goo.gl/maps/lHrrg
>>> geodesic(p1, p2).meters
429.1676644986777
>>> great_circle(p1, p2).meters
428.28877358686776

Mis disculpas por publicar una segunda respuesta aquí, pero aprovecho la oportunidad para responder a la solicitud de @ craig-hicks para proporcionar comparaciones de precisión y sincronización para varios algoritmos para calcular la distancia geodésica. Esto parafrasea un comentario que hago a mi solicitud de extracción # 144 para geopy, que permite el uso de una de las dos implementaciones de mi algoritmo para geodésicas para ser utilizadas dentro de geopy, una es una implementación nativa de python, geodésica (geografía) , y la otra utiliza una implementación en C, geodésica (pyproj) .

Aquí hay algunos datos de tiempo. Los tiempos están en microsecs por llamada

method                          dist    dest
geopy great_circle              20.4    17.1
geopy vincenty                  40.3    30.4
geopy geodesic(pyproj)          37.1    31.1
geopy geodesic(geographiclib)  302.9   124.1

Aquí está la precisión de los cálculos geodésicos basados ​​en mi conjunto de pruebas geodésicas . Los errores se dan en unidades de micras (1e-6 m)

method                        distance destination
geopy vincenty                 205.629  141.945
geopy geodesic(pyproj)           0.007    0.013
geopy geodesic(geographiclib)    0.011    0.010

He incluido la solicitud de extracción de Hannosche # 194 que corrige un error grave en la función de destino. Sin esta solución, el error en el cálculo del destino para vincenty es de 8,98 metros.

El 19.2% de los casos de prueba fallaron con vincenty.distance (iteraciones = 20). Sin embargo, el conjunto de prueba está sesgado hacia los casos que provocarían este fallo.

Con puntos aleatorios en el elipsoide WGS84, se garantiza que el algoritmo Vincenty fallará 16.6 de cada 1000000 veces (la solución correcta es un punto fijo inestable del método Vincenty).

Con la implementación de geopy de Vincenty e iteraciones = 20, la tasa de falla es 82.8 por 1000000. Con iteraciones = 200, la tasa de falla es 21.2 por 1000000.

Aunque estas tasas son pequeñas, las fallas pueden ser bastante comunes. Por ejemplo, en un conjunto de datos de 1000 puntos aleatorios (piense en los aeropuertos del mundo, tal vez), calcular la matriz de distancia completa fallaría en promedio 16 veces (con iteraciones = 20).

Parece que el paquete geopy.distance ofrece una función "distancia ()" que por defecto es vincenty (). Recomendaría usar distance () por principio, ya que es la recomendación del paquete, en caso de que alguna vez se desvíe de vincenty () en el futuro (por poco probable que sea). Sigue leyendo:

Esta nota de documentación se incluye en el código fuente de la función vincenty () que especificó:

Nota: Esta implementación de la distancia de Vincenty no converge para algunos puntos válidos. En algunos casos, se puede obtener un resultado al aumentar el número de iteraciones ( iterationsargumento de palabra clave, dado en la clase __init__, con un valor predeterminado de 20). Puede ser preferible usar: class:, .great_circleque es marginalmente menos preciso, pero siempre produce un resultado.

El código fuente con el comentario / nota anterior se puede encontrar en https://github.com/geopy/geopy/blob/master/geopy/distance.py Desplácese hacia abajo hasta la definición de vincenty ()

Sin embargo, la función de distancia predeterminada utilizada por ese paquete al calcular la distancia () es la función vincenty (), lo que implica que la falta de convergencia no es catastrófica y se devuelve una respuesta razonable; lo más importante es que no se genera una excepción.

Actualización: como lo señala "cffk", la función vincenty () lanza explícitamente una excepción ValueError cuando el algoritmo no converge, aunque no está documentado en la descripción de la función. Por lo tanto, la documentación tiene errores.

Ya sea que use vincenty o haversine o la ley esférica de cosenos, es prudente tomar conciencia de cualquier problema potencial con el código que planea usar, cosas que debe vigilar y mitigar, y cómo se trata con problemas de vincenty vs haversine vs sloc diferirá a medida que uno se dé cuenta de los problemas de cada uno que están al acecho / casos extremos, que pueden o no ser conocidos popularmente. El experimentado programador lo sabe. Los novatos no pueden. Espero evitarles la frustración cuando un fragmento de un foro hace algo inesperado, en ciertos casos. Si uno va a usar seriamente alguna versión de cualquiera de estos, vincenty, haversine, sloc, entonces SE, SO, Reddit, Quora, etc., pueden haber proporcionado ayuda limitada en la codificación inicial de una solución, pero eso no significa que su solución o 'respuesta' aceptada está libre de problemas. Si un proyecto es lo suficientemente importante, merece una cantidad razonable de investigación adecuada. Lea el manual, lea los documentos, y si existe una revisión de código de ese código, léalo. Copiar y pegar un fragmento o una esencia que se votó cientos o más veces no significa que su seguridad sea exhaustiva y segura.

La intrigante respuesta publicada por cffk plantea el punto de estar al tanto de los casos al acecho, en soluciones empaquetadas, que pueden producir excepciones u otras dificultades. . Las afirmaciones específicas hechas en esa publicación están más allá de mi presupuesto de tiempo para perseguir en la actualidad, pero deduzco que de hecho hay problemas al acecho en ciertos paquetes, incluida al menos una implementación de Vincent, respecto de la cual al menos una persona ha propuesto mejorar de una forma u otra, para minimizar o eliminar el riesgo de encontrar esas dificultades. No añadiré más a ese tema relacionado con vincenty (siendo demasiado ignorante de él), sino que me volveré a haversine, al menos en parte sobre el tema con el OP.

La fórmula de Haversine publicada popularmente, ya sea en Python u otro idioma, porque probablemente usará la especificación de punto flotante IEEE 754 en la mayoría de todos los sistemas de inteligencia e inteligencia actuales, y procesadores ARM, powerPC, etc. También es susceptible a errores de excepción raros pero reales y repetibles muy cerca o a una distancia de arco de 180 grados, puntos antipodales, debido a aproximaciones de punto flotante y redondeo. Es posible que algunos novatos aún no hayan sido mordidos por esta situación. Debido a que esta especificación de fp se aproxima y redondea, esto no significa que cualquier código que llame a fp64 pueda causar errores de excepción, no. Pero algo de código, Es posible que algunas fórmulas no tengan casos tan obvios en los que las aproximaciones y redondeos de IEEE 754 fp64 puedan hacer que un valor se salga ligeramente del dominio de un método matemático que se espera que evalúe sin problemas dicho valor. Un ejemplo ... sqrt (). Si un valor negativo llega a un sqrt (), como sqrt (-0.00000000000000000122739), habrá un error de excepción. En la fórmula de Haversine, la forma en que progresa hacia una solución, hay dos métodos sqrt () en atan2 (). losuna que se calcula y luego se usa en el sqrt (), puede, en los puntos antipodales del globo, ligeramente desviada por debajo de 0.0 o por encima de 1.0, muy ligeramente debido a las aproximaciones y redondeos de fp64, rara vez, pero de manera repetible. La repetibilidad confiable y consistente, en este contexto, hace de este un riesgo de excepción, un caso límite para proteger, mitigar, en lugar de una casualidad aleatoria aislada. Aquí hay un ejemplo de un fragmento corto de python3 de haversine, sin la protección necesaria:

import math as m

a = m.sin(dlat / 2)**2 + m.cos(lat1) * m.cos(lat2) * m.sin(dlon / 2)**2
c = 2 * m.atan2(m.sqrt(a), m.sqrt(1 - a))
distance = Radius * c

Muy cerca o en puntos antipodales, un cálculo en la primera línea de la fórmula puede desviarse negativamente, raramente, pero repetidamente con esas mismas coordenadas de lat lon. Para proteger / corregir esos raros casos, simplemente se puede añadir, después de la de un cálculo, como se ve a continuación:

import math as m

note = ''

a = m.sin(dlat / 2)**2 + m.cos(lat1) * m.cos(lat2) * m.sin(dlon / 2)**2
if a < 0.0: a = 0.0 ; note = '*'
if a > 1.0: a = 1.0 ; note = '**'
c = 2 * m.atan2(m.sqrt(a), m.sqrt(1 - a))
distance = Radius * c

# note = '*'  # a went below 0.0 and was normalized back to 0.0
# note = '**' # a went above 1.0 and was normalized back to max of 1.0

Por supuesto, no mostré toda la función aquí, sino un fragmento corto como se publica con tanta frecuencia. Pero este muestra la protección para el sqrt (), probando el a , y normalizándolo si es necesario, también ahorrando la necesidad de probar todo excepto. La nota = '' arriba es para evitar que la etapa de bytecode proteste contra esa nota se está utilizando antes de que se le asigne un valor, si se devuelve con el resultado de la función.

Con este simple cambio, de la adición de los dos a las pruebas, el sqrt () funciones estarán contentos, y el código de ahora tiene un adicional de nota que se puede devolver al código de llamada, que alerta de que el resultado ha sido ligeramente normalizado, y por qué. A algunos les puede interesar, a otros no, pero está ahí, evitando un error de excepción, que 'de otra manera' puede ocurrir. Un intento excepto el bloque puede detectar la excepción, pero no solucionarlo, a menos que se escriba explícitamente para hacerlo. Parece más fácil de código de corrección de la línea (s) inmediatamente después de la una línea de cálculo. La entrada completamente depurada no debería requerir un intento, excepto el bloqueo aquí.

Resumen, si usa haversine, codificado explícitamente en lugar de usar un paquete o una biblioteca, sin importar su idioma de elección, sería una buena idea probar y normalizar una vuelta al rango necesario de 0.0 <= a <= 1.0 en orden para proteger la siguiente línea con sus cálculos c . Pero la mayoría de los fragmentos de código de Haversine no lo muestran y no mencionan el riesgo.

Experiencia: durante pruebas exhaustivas en todo el mundo, en incrementos de 0.001 grados, llené un disco duro con combinaciones de lat lon que causaron una excepción, una excepción repetible confiable y consistente, durante un mes de pruebas colaterales de la confiabilidad del enfriamiento de la CPU fan, y mi paciencia. Sí, desde entonces he eliminado la mayoría de esos registros, ya que su propósito era principalmente probar el punto (si el juego de palabras está permitido). Pero tengo algunos registros más cortos de 'valores de problemas de lat lon', guardados para fines de prueba.

Precisión: ¿Será una y el resultado completo haversine perder cierta precisión mediante la normalización de esa vuelta pequeño fragmento en el dominio? No mucho, tal vez no más que las aproximaciones y redondeos fp64 que ya se estaban introduciendo, lo que provocó que esa ligera deriva fuera del dominio. Si ya ha encontrado aceptable a Haversine sobre Vincenty: más simple, más rápido, más fácil de personalizar, solucionar problemas y mantener, entonces Haversine puede ser una buena solución para su proyecto.

He usado haversine en una skysphere proyectada por encima de la cabeza para medir distancias angulares entre objetos en el cielo, tal como se ve desde una posición en la tierra, mapeando azimut y alt a coordenadas de la equivalencia de skysphere lat lon, sin elipsoide para tener en cuenta en absoluto, ya que La skysphere teórica proyectada es una esfera perfecta, cuando se trata de medir la distancia angular, mire los ángulos entre dos objetos desde una posición en la superficie de la tierra. Se adapta perfectamente a mis necesidades. Entonces, el haversine sigue siendo muy útil y muy preciso en ciertas aplicaciones (dentro de mis propósitos) ... pero si lo usa, ya sea en la tierra para SIG o navegación, o en observaciones y mediciones de objetos del cielo, proteja en el caso de puntos antipodales o muy cerca de puntos antipodales, probando uny empujándolo nuevamente a su dominio necesario cuando sea necesario.

El haversine desprotegido está en todo Internet, y solo he visto una publicación anterior de Usenet que mostró cierta protección, creo que alguien de JPL, y que puede haber sido anterior a 1985, anterior a la especificación de punto flotante IEEE 754. Otras dos páginas mencionaron posibles problemas cerca de los puntos antipodales, pero no describieron esos problemas ni cómo se podrían mitigarlos. Por lo tanto, existe una preocupación por los novatos (como yo) que no siempre entienden las buenas prácticas lo suficientemente bien como para seguir investigando y probar casos límite, de algún código que han copiado y pegado en un proyecto de confianza. La intrigante publicación de cffk fue refrescante, ya que era pública con este tipo de problemas, que a menudo no se mencionan, rara vez se codifican públicamente para protección en fragmentos, y rara vez se discuten de esta manera, en comparación con la cantidad de versiones no protegidas y no discutidas que se publican.

A partir de 20190923, la página wiki para la fórmula de Haversine de hecho menciona el problema posible en los puntos antipodales, debido a problemas de coma flotante en los dispositivos informáticos ... alentador ...

https://en.wikipedia.org/wiki/Haversine_formula

(debido a que esa página wiki no tiene, en este momento, un ancla html para la sección a la que me vincularía directamente, por lo tanto, después de que se cargue la página, haga una búsqueda en esa página del navegador para 'Al usar estas fórmulas' y usted vea el problema de la haversina con los puntos antipodales mencionados, más oficialmente).

Y este otro sitio también tiene una breve mención:

https://www.movable-type.co.uk/scripts/latlong.html

Si uno hace una búsqueda en esa página para 'incluir protección contra errores de redondeo', existe esto ...

Si atan2 no está disponible, c podría calcularse a partir de 2 ⋅ asin (min (1, √a)) (incluida la protección contra errores de redondeo).

Ahora hay una rara instancia en la que se mencionan errores de redondeo y se muestra protección para la versión asin (), pero no se menciona ni se muestra para la versión atan2 (). Pero al menos se menciona el riesgo de errores de redondeo.

En mi opinión, cualquier aplicación 24/7/365 que utilice Haversine necesita esta protección cerca de los puntos antipodales como un detalle importante y simple.

No sé qué paquetes de Haversine incluyen o no incluyen esta protección, pero si usted es nuevo en todo esto y va a utilizar las versiones 'snippet' publicadas popularmente, ahora sabe que necesita protección, y esa protección es muy simple de implementar, es decir, si no está utilizando vincenty, y no está utilizando un haversine empaquetado sin acceso fácil para modificar el código del paquete.

IOW, ya sea usando vincenty o haversine o sloc, uno debe ser consciente de cualquier problema con el código, cosas a tener en cuenta y mitigar, y la forma en que se trata de los problemas de vincenty vs haversine vs sloc será diferente a medida que uno se dé cuenta de cada uno problemas al acecho / casos extremos, que pueden o no ser conocidos popularmente.