Distancia entre dos puntos de la superficie terrestre mediante la fórmula de Haversine con Python

En el artículo anterior se determinó la distancia entre dos puntos de la superficie terrestre usando los valores de longitud y latitud (coordenadas geográficas en grados decimales) empleando como base para el cálculo el elipsoide WGS84. El resultado se comparó con la distancia euclideana determinada para los mismos puntos pero proyectados en coordenadas UTM zona 19N con el datum WGS84. La diferencia, en términos porcentuales, fue de 0,016 % (8.73 m en poco más de 56 Km). En este artículo se va a determinar la distancia entre los mismos puntos pero en lugar de usar el elipsoide se va a considerar la aproximación esférica de la superficie terrestre mediante la bien conocida fórmula Harvesine.

La formula Harvesine para la aproximación esférica de la distancia (d) entre dos puntos de la superficie terrestre se encuentra a continuación:

donde Φ1, Φ2 y λ1, λ2 se refieren a la latitud y a la longitud, expresadas ambas en radianes, de los puntos 1 y 2 respectivamente y r corresponde al radio terrestre (Ecuatorial 6378.1 km, Polar 6356.8 km, Medio 6371.0 km). El código Python para implementarla, conjuntamente con la determinación empleando la aproximación elipsoidal mediante el módulo pyproj, se encuentra a continuación:

#!/usr/bin/env/python
# -*- coding: utf-8

import pyproj, sys
from os import system
from math import sin,cos,sqrt,asin,pi

system("clear")

print "Elija la opción del elipsoide"

print "WGS84 = 1 GRS80 = 2 clrk66 = 3 intl = 4"
	
opcion = int(raw_input("opción = ? "))

long1 = -67.9766140755
lat1 = 9.36642904389
long2 = -67.4801781185
lat2 = 9.48364826587

if opcion == 1:
	name_ellps = "WGS84"

elif opcion == 2:
	name_ellps = "GRS80"

elif opcion == 3:
	name_ellps = "clrk66"

elif opcion == 4:
	name_ellps = "intl"

else:
	print "No existe la opción"
	sys.exit()

long1,lat1 = (long1,lat1)
long2,lat2 = (long2,lat2)
geod = pyproj.Geod(ellps=name_ellps)
angle1,angle2,distance = geod.inv(long1, lat1, long2, lat2)

print "La distancia es %0.2f metros basada en el elipsoide de" % distance, name_ellps

r = 6371000 #radio terrestre medio, en metros

c = pi/180 #constante para transformar grados en radianes

#Fórmula de haversine
d = 2*r*asin(sqrt(sin(c*(lat2-lat1)/2)**2 + cos(c*lat1)*cos(c*lat2)*sin(c*(long2-long1)/2)**2))

print "La distancia es %0.2f metros basada en la formula de haversine" % d

print "La diferencia, en valor absoluto, es %0.2f" % abs(distance-d), "metros"

print "La diferencia, en valor porcentual, es %0.2f" % abs((distance-d)*100/distance), "%"

La ejecución del código, en cónsola de bash, produce el siguiente resultado para la opción 1 (WGS84):

Elija la opción del elipsoide
WGS84 = 1 GRS80 = 2 clrk66 = 3 intl = 4
opción = ? 1
La distancia es 56042.15 metros basada en el elipsoide de WGS84
La distancia es 55994.13 metros basada en la formula de haversine
La diferencia, en valor absoluto, es 48.03 metros
La diferencia, en valor porcentual, es 0.09 %

Corresponde a 0.09 % (casi 48 m en aproximadamente 56 km); relativamente elevada si se le compara con el 0.016 % (8.73 m en aproximadamente 56 km) del procedimiento que emplea la proyección UTM.

El resultado fue corroborado empleado una aplicación en Java y la imagen del resultado se expone a continuación:

La aplicación produce el valor aquí obtenido para la distancia (en km) y posiciona los puntos en la zona correspondiente. No obstante, si les maravilla la espectacularidad de la aplicación web en Java, recuerden que el resultado obtenido mediante el módulo pyproj de Python sería el convencionalmente correcto.

Esta entrada fue publicada en Código Python, SIG, Software Libre. Guarda el enlace permanente.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s