Calculando la ruta más corta sobre un elipsoide con Python

Una geodésica (o línea geodésica) es el camino más corto posible entre dos puntos ubicados sobre una superficie dada. En el mundo de la cartografía y la geodesia, hallar el camino más corto entre dos ubicaciones sobre un elipsoide de revolución supone calcular la línea geodésica que los une (ver más en Wikipedia).

Considerar la Tierra como una esfera simplifica los cálculos, siendo el camino más corto el arco de círculo máximo entre ambos puntos (Wikipedia).

Desde hace algún tiempo, una de las principales librerías geo del ecosistema open source, Proj.4, ha venido introduciendo profundos cambios en su algoritmia de cálculos geodésicos. Concretamente, a partir de la versión 4.9 se ha portado a C dentro de Proj.4 parte de la librería C++ GeographicLib escrita por Charles Karney (más sobre los algoritmos de C. Karney).

Hasta la versión 4.8 esta librería venía utilizando los  algoritmos de Paul D. Thomas para ello. El intenso uso que hago diariamente de esta librería me ha conducido a realizar algunos testeos sobre cálculos de distancias geodésicas.

En este enlace pueden acceder a un benchmark que realizó Cayetano Benavent, con varias librerías para calcular la distancia geodésica entre diversas ubicaciones, conociendo las coordenadas de dichas ubicaciones (lo cual conlleva resolver el problema geodésico inverso): https://github.com/cayetanobv/GeodeticMusings

Además desarrolló una pequeña aplicación que, resolviendo el problema geodésico inverso con base en los algoritmos mencionados, calcula las líneas geodésicas entre los puntos de inicio y fin en formato GIS directamente (Shapefile y/o GeoJSON): GeodesicLinesToGIS.

Es importante resaltar que el programa solventa de forma eficiente el problema de líneas geodésicas que cruzan el antimeridiano a la hora de generar las geometrías en formato GIS. Este programa toma como base las excelentes librerías Pyproj, Fiona y Shapely.

En los ejemplos que mostramos a continuación se puede ver la diferencia entre la línea geodésica calculada (camino más corto; línea verde), la línea de rumbo (o loxodroma; línea roja) y una línea recta entre ambos puntos (línea discontinua negra).

Las diferencias máximas, como es lógico esperar, ocurren entre la proyección Mercator (ver Fig 1; la loxodroma es una línea recta) y la proyección Gnomónica (ver Fig 2; la geodésica es una línea recta). El problema de líneas geodésicas que cruzan la línea de 180º se observa solventado en las Fig 3 y 4.

Fig 1. Proyección Mercator. Líneas Geodésicas.

Fig 2. Proyección Gnomónica (centrada en 50W y 60N). Líneas Geodésicas

Fig 3. Líneas geodésicas que cruzan 180º, Proyección Mercator.

Fig 4. Líneas geodésicas que cruzan 180º, Proyección Mercator (Centrada en 150E).

Esta librería puede instalarse directamente desde el repositorio de aplicaciones de Python: https://pypi.python.org/pypi/GeodesicLinesToGIS

Todo el código fuente está en GiHub: https://github.com/GeographicaGS/GeodesicLinesToGIS

¡Hasta pronto!