Fórmula Haversine en Python (rumbo y distancia entre dos puntos GPS)

Problema

Me gustaría saber cómo obtener la distancia y el rumbo entre 2 puntos GPS . He investigado sobre la fórmula haversine. Alguien me dijo que también podría encontrar el rodamiento usando los mismos datos.

Editar

Todo funciona bien, pero el rodamiento no funciona bien todavía. El rodamiento produce resultados negativos, pero debe estar entre 0 y 360 grados. El conjunto de datos debe hacer el rodamiento horizontal 96.02166666666666 y es:

 Start point: 53.32055555555556 , -1.7297222222222221 Bearing: 96.02166666666666 Distance: 2 km Destination point: 53.31861111111111, -1.6997222222222223 Final bearing: 96.04555555555555 

Aquí está mi nuevo código:

 from math import * Aaltitude = 2000 Oppsite = 20000 lat1 = 53.32055555555556 lat2 = 53.31861111111111 lon1 = -1.7297222222222221 lon2 = -1.6997222222222223 lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2]) dlon = lon2 - lon1 dlat = lat2 - lat1 a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2 c = 2 * atan2(sqrt(a), sqrt(1-a)) Base = 6371 * c Bearing =atan2(cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(lon2-lon1), sin(lon2-lon1)*cos(lat2)) Bearing = degrees(Bearing) print "" print "" print "--------------------" print "Horizontal Distance:" print Base print "--------------------" print "Bearing:" print Bearing print "--------------------" Base2 = Base * 1000 distance = Base * 2 + Oppsite * 2 / 2 Caltitude = Oppsite - Aaltitude a = Oppsite/Base b = atan(a) c = degrees(b) distance = distance / 1000 print "The degree of vertical angle is:" print c print "--------------------" print "The distance between the Balloon GPS and the Antenna GPS is:" print distance print "--------------------" 

Aquí hay una versión de Python:

 from math import radians, cos, sin, asin, sqrt def haversine(lon1, lat1, lon2, lat2): """ Calculate the great circle distance between two points on the earth (specified in decimal degrees) """ # convert decimal degrees to radians lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2]) # haversine formula dlon = lon2 - lon1 dlat = lat2 - lat1 a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2 c = 2 * asin(sqrt(a)) r = 6371 # Radius of earth in kilometers. Use 3956 for miles return c * r 

La mayoría de estas respuestas son “redondeando” el radio de la tierra. Si los compara con otras calculadoras de distancia (como geopy), estas funciones estarán desactivadas.

Esto funciona bien:

 from math import radians, cos, sin, asin, sqrt def haversine(lat1, lon1, lat2, lon2): R = 3959.87433 # this is in miles. For Earth radius in kilometers use 6372.8 km dLat = radians(lat2 - lat1) dLon = radians(lon2 - lon1) lat1 = radians(lat1) lat2 = radians(lat2) a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2 c = 2*asin(sqrt(a)) return R * c # Usage lon1 = -103.548851 lat1 = 32.0004311 lon2 = -103.6041946 lat2 = 33.374939 print(haversine(lat1, lon1, lat2, lon2)) 

El cálculo del rumbo es incorrecto, necesita intercambiar las entradas a atan2.

  bearing = atan2(sin(long2-long1)*cos(lat2), cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(long2-long1)) bearing = degrees(bearing) bearing = (bearing + 360) % 360 

Esto le dará la orientación correcta.

Puedes probar lo siguiente:

 from haversine import haversine haversine((45.7597, 4.8422),(48.8567, 2.3508),miles = True) 243.71209416020253 

Puedes resolver el problema del rodamiento negativo añadiendo 360 °. Desafortunadamente, esto podría dar como resultado rodamientos de más de 360 ​​° para rodamientos positivos. Este es un buen candidato para el operador de módulo, así que en general debe agregar la línea.

 Bearing = (Bearing + 360) % 360 

Al final de tu método.

También hay una implementación vectorizada , que permite utilizar 4 matrices numpy en lugar de valores escalares para las coordenadas:

 def distance(s_lat, s_lng, e_lat, e_lng): # approximate radius of earth in km R = 6373.0 s_lat = s_lat*np.pi/180.0 s_lng = np.deg2rad(s_lng) e_lat = np.deg2rad(e_lat) e_lng = np.deg2rad(e_lng) d = np.sin((e_lat - s_lat)/2)**2 + np.cos(s_lat)*np.cos(e_lat) * np.sin((e_lng - s_lng)/2)**2 return 2 * R * np.arcsin(np.sqrt(d)) 

La Y en atan2 es, por defecto, el primer parámetro. Aquí está la documentación . Tendrá que cambiar sus entradas para obtener el ángulo de orientación correcto.

 bearing = atan2(sin(lon2-lon1)*cos(lat2), cos(lat1)*sin(lat2)in(lat1)*cos(lat2)*cos(lon2-lon1)) bearing = degrees(bearing) bearing = (bearing + 360) % 360 

Consulte este enlace: https://gis.stackexchange.com/questions/84885/whats-the-difference-between-vincenty-and-great-circle-distance-calculations

Esto en realidad da dos formas de distanciarse. Ellos son Haversine y Vincentys. De mi investigación llegué a saber que Vincentys es relativamente preciso. También use la statement de importación para hacer la implementación.

Aquí hay dos funciones para calcular la distancia y el rumbo, que se basan en el código de los mensajes anteriores y https://gist.github.com/jeromer/2005586 (tipo de tupla agregada para puntos geográficos en formato lat, lon para ambas funciones para mayor claridad) ). Probé ambas funciones y parecen funcionar bien.

 #coding:UTF-8 from math import radians, cos, sin, asin, sqrt, atan2, degrees def haversine(pointA, pointB): if (type(pointA) != tuple) or (type(pointB) != tuple): raise TypeError("Only tuples are supported as arguments") lat1 = pointA[0] lon1 = pointA[1] lat2 = pointB[0] lon2 = pointB[1] # convert decimal degrees to radians lat1, lon1, lat2, lon2 = map(radians, [lat1, lon1, lat2, lon2]) # haversine formula dlon = lon2 - lon1 dlat = lat2 - lat1 a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2 c = 2 * asin(sqrt(a)) r = 6371 # Radius of earth in kilometers. Use 3956 for miles return c * r def initial_bearing(pointA, pointB): if (type(pointA) != tuple) or (type(pointB) != tuple): raise TypeError("Only tuples are supported as arguments") lat1 = radians(pointA[0]) lat2 = radians(pointB[0]) diffLong = radians(pointB[1] - pointA[1]) x = sin(diffLong) * cos(lat2) y = cos(lat1) * sin(lat2) - (sin(lat1) * cos(lat2) * cos(diffLong)) initial_bearing = atan2(x, y) # Now we have the initial bearing but math.atan2 return values # from -180° to + 180° which is not what we want for a compass bearing # The solution is to normalize the initial bearing as shown below initial_bearing = degrees(initial_bearing) compass_bearing = (initial_bearing + 360) % 360 return compass_bearing pA = (46.2038,6.1530) pB = (46.449, 30.690) print haversine(pA, pB) print initial_bearing(pA, pB) 

Aquí hay una implementación numérica vectorizada de la Fórmula Haversine dada por @Michael Dunn, que ofrece una mejora de 10 a 50 veces mayor que en vectores grandes.

 from numpy import radians, cos, sin, arcsin, sqrt def haversine(lon1, lat1, lon2, lat2): """ Calculate the great circle distance between two points on the earth (specified in decimal degrees) """ #Convert decimal degrees to Radians: lon1 = np.radians(lon1.values) lat1 = np.radians(lat1.values) lon2 = np.radians(lon2.values) lat2 = np.radians(lat2.values) #Implementing Haversine Formula: dlon = np.subtract(lon2, lon1) dlat = np.subtract(lat2, lat1) a = np.add(np.power(np.sin(np.divide(dlat, 2)), 2), np.multiply(np.cos(lat1), np.multiply(np.cos(lat2), np.power(np.sin(np.divide(dlon, 2)), 2)))) c = np.multiply(2, np.arcsin(np.sqrt(a))) r = 6371 return c*r