Calcular la distancia entre 2 coordenadas GPS

Respuestas:

406

Calcule la distancia entre dos coordenadas por latitud y longitud , incluida una implementación de Javascript.

Las ubicaciones oeste y sur son negativas. Recuerde que los minutos y segundos están fuera de 60, entonces S31 30 'es -31.50 grados.

No olvides convertir grados a radianes . Muchos idiomas tienen esta función. O su cálculo simple: radians = degrees * PI / 180.

function degreesToRadians(degrees) {
  return degrees * Math.PI / 180;
}

function distanceInKmBetweenEarthCoordinates(lat1, lon1, lat2, lon2) {
  var earthRadiusKm = 6371;

  var dLat = degreesToRadians(lat2-lat1);
  var dLon = degreesToRadians(lon2-lon1);

  lat1 = degreesToRadians(lat1);
  lat2 = degreesToRadians(lat2);

  var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
          Math.sin(dLon/2) * Math.sin(dLon/2) * Math.cos(lat1) * Math.cos(lat2); 
  var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); 
  return earthRadiusKm * c;
}

Aquí hay algunos ejemplos de uso:

distanceInKmBetweenEarthCoordinates(0,0,0,0)  // Distance between same 
                                              // points should be 0
0

distanceInKmBetweenEarthCoordinates(51.5, 0, 38.8, -77.1) // From London
                                                          // to Arlington
5918.185064088764
cletus
fuente
17
En caso de que no es obvio, el método toRad () es una personalización para el número de prototipos, tales como: Number.prototype.toRad = function() { return this * (Math.PI / 180); }; . O, como se indica a continuación, puede reemplazar (Math.PI/2)con 0.0174532925199433 (... cualquier precisión que considere necesaria) para un mayor rendimiento.
Vinney Kelly
44
Si alguien, específicamente aquellos de ustedes que no buscan comentarios de fin de línea, está mirando esta fórmula y buscando una unidad de distancia, la unidad es km. :)
Dylan Knowles
1
@VinneyKelly Error tipográfico pequeño pero reemplaza (Math.PI / 180) no (Math.PI / 2), gracias por la ayuda de todos
Patrick Murphy
1
@ChristianKRider Mira la primera línea. Piensa en lo que Rnormalmente significa en matemáticas, luego busca cantidades relevantes relacionadas con la Tierra para ver si los números coinciden.
Financia la demanda de Mónica el
3
Para las unidades imperiales (millas) podrías cambiar earthRadiusKmpara ser var earthRadiusMiles = 3959;, para tu información.
chapeljuice
59

Busque haversine con Google; Aquí está mi solución:

#include <math.h>
#include "haversine.h"

#define d2r (M_PI / 180.0)

//calculate haversine distance for linear distance
double haversine_km(double lat1, double long1, double lat2, double long2)
{
    double dlong = (long2 - long1) * d2r;
    double dlat = (lat2 - lat1) * d2r;
    double a = pow(sin(dlat/2.0), 2) + cos(lat1*d2r) * cos(lat2*d2r) * pow(sin(dlong/2.0), 2);
    double c = 2 * atan2(sqrt(a), sqrt(1-a));
    double d = 6367 * c;

    return d;
}

double haversine_mi(double lat1, double long1, double lat2, double long2)
{
    double dlong = (long2 - long1) * d2r;
    double dlat = (lat2 - lat1) * d2r;
    double a = pow(sin(dlat/2.0), 2) + cos(lat1*d2r) * cos(lat2*d2r) * pow(sin(dlong/2.0), 2);
    double c = 2 * atan2(sqrt(a), sqrt(1-a));
    double d = 3956 * c; 

    return d;
}
Peter Greis
fuente
3
Puede reemplazar (M_PI / 180.0) con 0.0174532925199433 para un mejor rendimiento.
Hlung
3
En términos de rendimiento: uno podría calcular sin (dlat / 2.0) solo una vez, almacenarlo en la variable a1, y en lugar de pow (, 2) es MUCHO mejor usar a1 * a1. Lo mismo para el otro pow (, 2).
pms
71
Sí, o simplemente use un compilador posterior a los años 60.
derecha el
17
No hay necesidad de "optimizar" (M_PI / 180.0) a una constante que nadie entiende sin contexto. ¡El compilador calcula estos términos fijos para usted!
Patrick Cornelissen
2
@ TõnuSamuel Muchas gracias por tu comentario. Realmente lo aprecio. Tiene sentido que el compilador con optimización habilitada (-O) pueda calcular previamente las operaciones de las constantes, haciendo que el colapso manual sea inútil. Lo probaré cuando tenga tiempo.
Hlung
44

Versión C # de Haversine

double _eQuatorialEarthRadius = 6378.1370D;
double _d2r = (Math.PI / 180D);

private int HaversineInM(double lat1, double long1, double lat2, double long2)
{
    return (int)(1000D * HaversineInKM(lat1, long1, lat2, long2));
}

private double HaversineInKM(double lat1, double long1, double lat2, double long2)
{
    double dlong = (long2 - long1) * _d2r;
    double dlat = (lat2 - lat1) * _d2r;
    double a = Math.Pow(Math.Sin(dlat / 2D), 2D) + Math.Cos(lat1 * _d2r) * Math.Cos(lat2 * _d2r) * Math.Pow(Math.Sin(dlong / 2D), 2D);
    double c = 2D * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1D - a));
    double d = _eQuatorialEarthRadius * c;

    return d;
}

Aquí hay un violín de .NET de esto , para que pueda probarlo con sus propios Lat / Longs.

Roman Makarov
fuente
1
También he agregado un violín de .NET para que las personas puedan probarlo fácilmente.
Pure.Krome
77
.Net Framework tiene un método incorporado en GeoCoordinate.GetDistanceTo. Se debe hacer referencia al ensamblado System.Device. Artículo de MSDN msdn.microsoft.com/en-us/library/…
fnx
27

Versión Java del algoritmo de Haversine basado en la respuesta de Roman Makarov a este hilo

public class HaversineAlgorithm {

    static final double _eQuatorialEarthRadius = 6378.1370D;
    static final double _d2r = (Math.PI / 180D);

    public static int HaversineInM(double lat1, double long1, double lat2, double long2) {
        return (int) (1000D * HaversineInKM(lat1, long1, lat2, long2));
    }

    public static double HaversineInKM(double lat1, double long1, double lat2, double long2) {
        double dlong = (long2 - long1) * _d2r;
        double dlat = (lat2 - lat1) * _d2r;
        double a = Math.pow(Math.sin(dlat / 2D), 2D) + Math.cos(lat1 * _d2r) * Math.cos(lat2 * _d2r)
                * Math.pow(Math.sin(dlong / 2D), 2D);
        double c = 2D * Math.atan2(Math.sqrt(a), Math.sqrt(1D - a));
        double d = _eQuatorialEarthRadius * c;

        return d;
    }

}
Paulo Miguel Almeida
fuente
@Radu asegúrate de usarlo correctamente y no intercambiar lugares lat / log al pasarlos a cualquier método.
Paulo Miguel Almeida
1
Obtuve una respuesta razonablemente cercana usando esta fórmula. Basé la precisión usando este sitio web: movable-type.co.uk/scripts/latlong.html que me dio 0.07149km, mientras que su fórmula me dio 0.07156una precisión de aproximadamente 99%
Janac Meena
24

Esto es muy fácil de hacer con el tipo de geografía en SQL Server 2008.

SELECT geography::Point(lat1, lon1, 4326).STDistance(geography::Point(lat2, lon2, 4326))
-- computes distance in meters using eliptical model, accurate to the mm

4326 es SRID para el modelo de tierra elipsoidal WGS84

Marko Tintor
fuente
19

Aquí hay una función de Haversine en Python que uso:

from math import pi,sqrt,sin,cos,atan2

def haversine(pos1, pos2):
    lat1 = float(pos1['lat'])
    long1 = float(pos1['long'])
    lat2 = float(pos2['lat'])
    long2 = float(pos2['long'])

    degree_to_rad = float(pi / 180.0)

    d_lat = (lat2 - lat1) * degree_to_rad
    d_long = (long2 - long1) * degree_to_rad

    a = pow(sin(d_lat / 2), 2) + cos(lat1 * degree_to_rad) * cos(lat2 * degree_to_rad) * pow(sin(d_long / 2), 2)
    c = 2 * atan2(sqrt(a), sqrt(1 - a))
    km = 6367 * c
    mi = 3956 * c

    return {"km":km, "miles":mi}
PaulMcG
fuente
11

Aquí está en C # (lat y long en radianes):

double CalculateGreatCircleDistance(double lat1, double long1, double lat2, double long2, double radius)
{
    return radius * Math.Acos(
        Math.Sin(lat1) * Math.Sin(lat2)
        + Math.Cos(lat1) * Math.Cos(lat2) * Math.Cos(long2 - long1));
}

Si su lat y long están en grados, divida por 180 / PI para convertir a radianes.

Mike Chamberlain
fuente
1
Este es el cálculo de la "ley esférica del coseno", que es el método de cálculo menos preciso y más propenso a errores de una gran distancia circular.
John Machin
11

Necesitaba calcular muchas distancias entre los puntos para mi proyecto, así que seguí adelante e intenté optimizar el código que encontré aquí. En promedio, en diferentes navegadores, mi nueva implementación se ejecuta 2 veces más rápido que la respuesta más votada.

function distance(lat1, lon1, lat2, lon2) {
  var p = 0.017453292519943295;    // Math.PI / 180
  var c = Math.cos;
  var a = 0.5 - c((lat2 - lat1) * p)/2 + 
          c(lat1 * p) * c(lat2 * p) * 
          (1 - c((lon2 - lon1) * p))/2;

  return 12742 * Math.asin(Math.sqrt(a)); // 2 * R; R = 6371 km
}

Puedes jugar con mi jsPerf y ver los resultados aquí. .

Recientemente necesitaba hacer lo mismo en Python, así que aquí hay una implementación de Python :

from math import cos, asin, sqrt
def distance(lat1, lon1, lat2, lon2):
    p = 0.017453292519943295
    a = 0.5 - cos((lat2 - lat1) * p)/2 + cos(lat1 * p) * cos(lat2 * p) * (1 - cos((lon2 - lon1) * p)) / 2
    return 12742 * asin(sqrt(a))

Y en aras de la integridad: Haversine en wiki.

Salvador Dalí
fuente
11

Versión PHP:

(Elimine todo deg2rad()si sus coordenadas ya están en radianes).

$R = 6371; // km
$dLat = deg2rad($lat2-$lat1);
$dLon = deg2rad($lon2-$lon1);
$lat1 = deg2rad($lat1);
$lat2 = deg2rad($lat2);

$a = sin($dLat/2) * sin($dLat/2) +
     sin($dLon/2) * sin($dLon/2) * cos($lat1) * cos($lat2); 

$c = 2 * atan2(sqrt($a), sqrt(1-$a)); 
$d = $R * $c;
quape
fuente
1
Cambie lat1 y lat2 a $ lat1 nad $ lat2.
lugar, el
7

Una función T-SQL, que uso para seleccionar registros por distancia para un centro

Create Function  [dbo].[DistanceInMiles] 
 (  @fromLatitude float ,
    @fromLongitude float ,
    @toLatitude float, 
    @toLongitude float
  )
   returns float
AS 
BEGIN
declare @distance float

select @distance = cast((3963 * ACOS(round(COS(RADIANS(90-@fromLatitude))*COS(RADIANS(90-@toLatitude))+ 
SIN(RADIANS(90-@fromLatitude))*SIN(RADIANS(90-@toLatitude))*COS(RADIANS(@fromLongitude-@toLongitude)),15)) 
)as float) 
  return  round(@distance,1)
END
Henry Vilinskiy
fuente
Este es el cálculo de la "ley esférica del coseno", que es el método de cálculo menos preciso y más propenso a errores de una gran distancia circular.
John Machin
5

Si necesita algo más preciso, eche un vistazo a esto .

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.

El primer método (directo) calcula la ubicación de un punto que es una distancia y azimut (dirección) dados desde otro punto. El segundo método (inverso) calcula la distancia geográfica y el acimut entre dos puntos dados. Han sido ampliamente utilizados en geodesia porque son precisos dentro de 0.5 mm (0.020 ″) en el elipsoide de la Tierra.

dsmelser
fuente
5

I. Sobre el método de "migas de pan"

  1. El radio de la Tierra es diferente en diferentes Lat. Esto debe tenerse en cuenta en el algoritmo de Haversine.
  2. Considere el cambio de rumbo, que convierte las líneas rectas en arcos (que son más largos)
  3. Tener en cuenta el cambio de Velocidad convertirá los arcos en espirales (que son más largos o más cortos que los arcos)
  4. El cambio de altitud convertirá las espirales planas en espirales 3D (que son más largas nuevamente). Esto es muy importante para las zonas montañosas.

A continuación, vea la función en C que tiene en cuenta el n. ° 1 y n. ° 2:

double   calcDistanceByHaversine(double rLat1, double rLon1, double rHeading1,
       double rLat2, double rLon2, double rHeading2){
  double rDLatRad = 0.0;
  double rDLonRad = 0.0;
  double rLat1Rad = 0.0;
  double rLat2Rad = 0.0;
  double a = 0.0;
  double c = 0.0;
  double rResult = 0.0;
  double rEarthRadius = 0.0;
  double rDHeading = 0.0;
  double rDHeadingRad = 0.0;

  if ((rLat1 < -90.0) || (rLat1 > 90.0) || (rLat2 < -90.0) || (rLat2 > 90.0)
              || (rLon1 < -180.0) || (rLon1 > 180.0) || (rLon2 < -180.0)
              || (rLon2 > 180.0)) {
        return -1;
  };

  rDLatRad = (rLat2 - rLat1) * DEGREE_TO_RADIANS;
  rDLonRad = (rLon2 - rLon1) * DEGREE_TO_RADIANS;
  rLat1Rad = rLat1 * DEGREE_TO_RADIANS;
  rLat2Rad = rLat2 * DEGREE_TO_RADIANS;

  a = sin(rDLatRad / 2) * sin(rDLatRad / 2) + sin(rDLonRad / 2) * sin(
              rDLonRad / 2) * cos(rLat1Rad) * cos(rLat2Rad);

  if (a == 0.0) {
        return 0.0;
  }

  c = 2 * atan2(sqrt(a), sqrt(1 - a));
  rEarthRadius = 6378.1370 - (21.3847 * 90.0 / ((fabs(rLat1) + fabs(rLat2))
              / 2.0));
  rResult = rEarthRadius * c;

  // Chord to Arc Correction based on Heading changes. Important for routes with many turns and U-turns

  if ((rHeading1 >= 0.0) && (rHeading1 < 360.0) && (rHeading2 >= 0.0)
              && (rHeading2 < 360.0)) {
        rDHeading = fabs(rHeading1 - rHeading2);
        if (rDHeading > 180.0) {
              rDHeading -= 180.0;
        }
        rDHeadingRad = rDHeading * DEGREE_TO_RADIANS;
        if (rDHeading > 5.0) {
              rResult = rResult * (rDHeadingRad / (2.0 * sin(rDHeadingRad / 2)));
        } else {
              rResult = rResult / cos(rDHeadingRad);
        }
  }
  return rResult;
}

II Hay una manera más fácil que da muy buenos resultados.

Por velocidad media.

Trip_distance = Trip_average_speed * Trip_time

Dado que la velocidad del GPS es detectada por el efecto Doppler y no está directamente relacionada con [Lon, Lat], al menos se puede considerar como secundaria (respaldo o corrección) si no como método de cálculo de distancia principal.

Tod Samay
fuente
4

Si está utilizando .NET, no reviva la rueda. Ver System.Device.Location . Crédito a fnx en los comentarios en otra respuesta .

using System.Device.Location;

double lat1 = 45.421527862548828D;
double long1 = -75.697189331054688D;
double lat2 = 53.64135D;
double long2 = -113.59273D;

GeoCoordinate geo1 = new GeoCoordinate(lat1, long1);
GeoCoordinate geo2 = new GeoCoordinate(lat2, long2);

double distance = geo1.GetDistanceTo(geo2);
Tim Perdiz
fuente
3

Este código Lua está adaptado de lo que se encuentra en Wikipedia y en la herramienta GPSbabel de Robert Lipe :

local EARTH_RAD = 6378137.0 
  -- earth's radius in meters (official geoid datum, not 20,000km / pi)

local radmiles = EARTH_RAD*100.0/2.54/12.0/5280.0;
  -- earth's radius in miles

local multipliers = {
  radians = 1, miles = radmiles, mi = radmiles, feet = radmiles * 5280,
  meters = EARTH_RAD, m = EARTH_RAD, km = EARTH_RAD / 1000, 
  degrees = 360 / (2 * math.pi), min = 60 * 360 / (2 * math.pi)
}

function gcdist(pt1, pt2, units) -- return distance in radians or given units
  --- this formula works best for points close together or antipodal
  --- rounding error strikes when distance is one-quarter Earth's circumference
  --- (ref: wikipedia Great-circle distance)
  if not pt1.radians then pt1 = rad(pt1) end
  if not pt2.radians then pt2 = rad(pt2) end
  local sdlat = sin((pt1.lat - pt2.lat) / 2.0);
  local sdlon = sin((pt1.lon - pt2.lon) / 2.0);
  local res = sqrt(sdlat * sdlat + cos(pt1.lat) * cos(pt2.lat) * sdlon * sdlon);
  res = res > 1 and 1 or res < -1 and -1 or res
  res = 2 * asin(res);
  if units then return res * assert(multipliers[units])
  else return res
  end
end
Norman Ramsey
fuente
3
    private double deg2rad(double deg)
    {
        return (deg * Math.PI / 180.0);
    }

    private double rad2deg(double rad)
    {
        return (rad / Math.PI * 180.0);
    }

    private double GetDistance(double lat1, double lon1, double lat2, double lon2)
    {
        //code for Distance in Kilo Meter
        double theta = lon1 - lon2;
        double dist = Math.Sin(deg2rad(lat1)) * Math.Sin(deg2rad(lat2)) + Math.Cos(deg2rad(lat1)) * Math.Cos(deg2rad(lat2)) * Math.Cos(deg2rad(theta));
        dist = Math.Abs(Math.Round(rad2deg(Math.Acos(dist)) * 60 * 1.1515 * 1.609344 * 1000, 0));
        return (dist);
    }

    private double GetDirection(double lat1, double lon1, double lat2, double lon2)
    {
        //code for Direction in Degrees
        double dlat = deg2rad(lat1) - deg2rad(lat2);
        double dlon = deg2rad(lon1) - deg2rad(lon2);
        double y = Math.Sin(dlon) * Math.Cos(lat2);
        double x = Math.Cos(deg2rad(lat1)) * Math.Sin(deg2rad(lat2)) - Math.Sin(deg2rad(lat1)) * Math.Cos(deg2rad(lat2)) * Math.Cos(dlon);
        double direct = Math.Round(rad2deg(Math.Atan2(y, x)), 0);
        if (direct < 0)
            direct = direct + 360;
        return (direct);
    }

    private double GetSpeed(double lat1, double lon1, double lat2, double lon2, DateTime CurTime, DateTime PrevTime)
    {
        //code for speed in Kilo Meter/Hour
        TimeSpan TimeDifference = CurTime.Subtract(PrevTime);
        double TimeDifferenceInSeconds = Math.Round(TimeDifference.TotalSeconds, 0);
        double theta = lon1 - lon2;
        double dist = Math.Sin(deg2rad(lat1)) * Math.Sin(deg2rad(lat2)) + Math.Cos(deg2rad(lat1)) * Math.Cos(deg2rad(lat2)) * Math.Cos(deg2rad(theta));
        dist = rad2deg(Math.Acos(dist)) * 60 * 1.1515 * 1.609344;
        double Speed = Math.Abs(Math.Round((dist / Math.Abs(TimeDifferenceInSeconds)) * 60 * 60, 0));
        return (Speed);
    }

    private double GetDuration(DateTime CurTime, DateTime PrevTime)
    {
        //code for speed in Kilo Meter/Hour
        TimeSpan TimeDifference = CurTime.Subtract(PrevTime);
        double TimeDifferenceInSeconds = Math.Abs(Math.Round(TimeDifference.TotalSeconds, 0));
        return (TimeDifferenceInSeconds);
    }
Elanchezhian Babu P
fuente
1
Creo que su función GetDistance devuelve valor en metros
Przemek
¿Es esto correcto? GetDirection () no utiliza 'dlat'.
gordito
3

Esta es la versión de "Henry Vilinskiy" adaptada para MySQL y Kilómetros:

CREATE FUNCTION `CalculateDistanceInKm`(
  fromLatitude float,
  fromLongitude float,
  toLatitude float, 
  toLongitude float
) RETURNS float
BEGIN
  declare distance float;

  select 
    6367 * ACOS(
            round(
              COS(RADIANS(90-fromLatitude)) *
                COS(RADIANS(90-toLatitude)) +
                SIN(RADIANS(90-fromLatitude)) *
                SIN(RADIANS(90-toLatitude)) *
                COS(RADIANS(fromLongitude-toLongitude))
              ,15)
            )
    into distance;

  return  round(distance,3);
END;
Maxs
fuente
MySQLdichoSomething is wrong in your syntax near '' on line 8 // declare distance float;
Legionar
Este es el cálculo de la "ley esférica del coseno", que es el método de cálculo menos preciso y más propenso a errores de una gran distancia circular
John Machin
3

aquí está la implementación Swift de la respuesta

func degreesToRadians(degrees: Double) -> Double {
    return degrees * Double.pi / 180
}

func distanceInKmBetweenEarthCoordinates(lat1: Double, lon1: Double, lat2: Double, lon2: Double) -> Double {

    let earthRadiusKm: Double = 6371

    let dLat = degreesToRadians(degrees: lat2 - lat1)
    let dLon = degreesToRadians(degrees: lon2 - lon1)

    let lat1 = degreesToRadians(degrees: lat1)
    let lat2 = degreesToRadians(degrees: lat2)

    let a = sin(dLat/2) * sin(dLat/2) +
    sin(dLon/2) * sin(dLon/2) * cos(lat1) * cos(lat2)
    let c = 2 * atan2(sqrt(a), sqrt(1 - a))
    return earthRadiusKm * c
}
Sai Li
fuente
3

tomé la respuesta principal y la usé en un programa Scala

import java.lang.Math.{atan2, cos, sin, sqrt}

def latLonDistance(lat1: Double, lon1: Double)(lat2: Double, lon2: Double): Double = {
    val earthRadiusKm = 6371
    val dLat = (lat2 - lat1).toRadians
    val dLon = (lon2 - lon1).toRadians
    val latRad1 = lat1.toRadians
    val latRad2 = lat2.toRadians

    val a = sin(dLat / 2) * sin(dLat / 2) + sin(dLon / 2) * sin(dLon / 2) * cos(latRad1) * cos(latRad2)
    val c = 2 * atan2(sqrt(a), sqrt(1 - a))
    earthRadiusKm * c
}

Curri la función para poder producir fácilmente funciones que tengan una de las dos ubicaciones fijas y requieran solo un par de lat / lon para producir la distancia.

Peter Perháč
fuente
2

Supongo que lo quieres a lo largo de la curvatura de la tierra. Tus dos puntos y el centro de la tierra están en un plano. El centro de la tierra es el centro de un círculo en ese plano y los dos puntos están (aproximadamente) en el perímetro de ese círculo. A partir de eso, puede calcular la distancia descubriendo cuál es el ángulo de un punto a otro.

Si los puntos no son de la misma altura, o si necesita tener en cuenta que la tierra no es una esfera perfecta, se vuelve un poco más difícil.

Guge
fuente
2

Recientemente tuve que hacer lo mismo. Encontré que este sitio web es muy útil al explicar trigonometría esférica con ejemplos que fueron fáciles de seguir.

Twotymz
fuente
2

puede encontrar una implementación de esto (con una buena explicación) en F # en fssnip

Aquí están las partes importantes:


let GreatCircleDistance<[<Measure>] 'u> (R : float<'u>) (p1 : Location) (p2 : Location) =
    let degToRad (x : float<deg>) = System.Math.PI * x / 180.0<deg/rad>

    let sq x = x * x
    // take the sin of the half and square the result
    let sinSqHf (a : float<rad>) = (System.Math.Sin >> sq) (a / 2.0<rad>)
    let cos (a : float<deg>) = System.Math.Cos (degToRad a / 1.0<rad>)

    let dLat = (p2.Latitude - p1.Latitude) |> degToRad
    let dLon = (p2.Longitude - p1.Longitude) |> degToRad

    let a = sinSqHf dLat + cos p1.Latitude * cos p2.Latitude * sinSqHf dLon
    let c = 2.0 * System.Math.Atan2(System.Math.Sqrt(a), System.Math.Sqrt(1.0-a))

    R * c
Carsten
fuente
2

Necesitaba implementar esto en PowerShell, espero que pueda ayudar a alguien más. Algunas notas sobre este método

  1. No divida ninguna de las líneas o el cálculo será incorrecto.
  2. Para calcular en KM, elimine el * 1000 en el cálculo de $ distancia
  3. Cambie $ earthsRadius = 3963.19059 y elimine * 1000 en el cálculo de $ distancia para calcular la distancia en millas
  4. Estoy usando Haversine, ya que otras publicaciones han señalado que las fórmulas de Vincenty son mucho más precisas

    Function MetresDistanceBetweenTwoGPSCoordinates($latitude1, $longitude1, $latitude2, $longitude2)  
    {  
      $Rad = ([math]::PI / 180);  
    
      $earthsRadius = 6378.1370 # Earth's Radius in KM  
      $dLat = ($latitude2 - $latitude1) * $Rad  
      $dLon = ($longitude2 - $longitude1) * $Rad  
      $latitude1 = $latitude1 * $Rad  
      $latitude2 = $latitude2 * $Rad  
    
      $a = [math]::Sin($dLat / 2) * [math]::Sin($dLat / 2) + [math]::Sin($dLon / 2) * [math]::Sin($dLon / 2) * [math]::Cos($latitude1) * [math]::Cos($latitude2)  
      $c = 2 * [math]::ATan2([math]::Sqrt($a), [math]::Sqrt(1-$a))  
    
      $distance = [math]::Round($earthsRadius * $c * 1000, 0) #Multiple by 1000 to get metres  
    
      Return $distance  
    }
    
TheLukeMcCarthy
fuente
2

Versión Scala

  def deg2rad(deg: Double) = deg * Math.PI / 180.0

  def rad2deg(rad: Double) = rad / Math.PI * 180.0

  def getDistanceMeters(lat1: Double, lon1: Double, lat2: Double, lon2: Double) = {
    val theta = lon1 - lon2
    val dist = Math.sin(deg2rad(lat1)) * Math.sin(deg2rad(lat2)) + Math.cos(deg2rad(lat1)) *
      Math.cos(deg2rad(lat2)) * Math.cos(deg2rad(theta))
    Math.abs(
      Math.round(
        rad2deg(Math.acos(dist)) * 60 * 1.1515 * 1.609344 * 1000)
    )
  }
Przemek
fuente
1

// ¿Quizás un error tipográfico?
Tenemos un dlon variable no utilizado en GetDirection,
supongo

double y = Math.Sin(dlon) * Math.Cos(lat2);
// cannot use degrees in Cos ?

debiera ser

double y = Math.Sin(dlon) * Math.Cos(dlat);
Gerrie de Jager
fuente
1
Esta no es una respuesta, es, en el mejor de los casos, un comentario.
Kevin
1

Aquí está mi implementación en Elixir

defmodule Geo do
  @earth_radius_km 6371
  @earth_radius_sm 3958.748
  @earth_radius_nm 3440.065
  @feet_per_sm 5280

  @d2r :math.pi / 180

  def deg_to_rad(deg), do: deg * @d2r

  def great_circle_distance(p1, p2, :km), do: haversine(p1, p2) * @earth_radius_km
  def great_circle_distance(p1, p2, :sm), do: haversine(p1, p2) * @earth_radius_sm
  def great_circle_distance(p1, p2, :nm), do: haversine(p1, p2) * @earth_radius_nm
  def great_circle_distance(p1, p2, :m), do: great_circle_distance(p1, p2, :km) * 1000
  def great_circle_distance(p1, p2, :ft), do: great_circle_distance(p1, p2, :sm) * @feet_per_sm

  @doc """
  Calculate the [Haversine](https://en.wikipedia.org/wiki/Haversine_formula)
  distance between two coordinates. Result is in radians. This result can be
  multiplied by the sphere's radius in any unit to get the distance in that unit.
  For example, multiple the result of this function by the Earth's radius in
  kilometres and you get the distance between the two given points in kilometres.
  """
  def haversine({lat1, lon1}, {lat2, lon2}) do
    dlat = deg_to_rad(lat2 - lat1)
    dlon = deg_to_rad(lon2 - lon1)

    radlat1 = deg_to_rad(lat1)
    radlat2 = deg_to_rad(lat2)

    a = :math.pow(:math.sin(dlat / 2), 2) +
        :math.pow(:math.sin(dlon / 2), 2) *
        :math.cos(radlat1) * :math.cos(radlat2)

    2 * :math.atan2(:math.sqrt(a), :math.sqrt(1 - a))
  end
end
mroach
fuente
1

Dart Version

Algoritmo de Haversine.

import 'dart:math';

class GeoUtils {

  static double _degreesToRadians(degrees) {
    return degrees * pi / 180;
  }

  static double distanceInKmBetweenEarthCoordinates(lat1, lon1, lat2, lon2) {
    var earthRadiusKm = 6371;

    var dLat = _degreesToRadians(lat2-lat1);
    var dLon = _degreesToRadians(lon2-lon1);

    lat1 = _degreesToRadians(lat1);
    lat2 = _degreesToRadians(lat2);

    var a = sin(dLat/2) * sin(dLat/2) +
        sin(dLon/2) * sin(dLon/2) * cos(lat1) * cos(lat2);
    var c = 2 * atan2(sqrt(a), sqrt(1-a));
    return earthRadiusKm * c;
  }
}
abd3llatif
fuente
0

Creo que todavía falta una versión del algoritmo en R :

gpsdistance<-function(lat1,lon1,lat2,lon2){

# internal function to change deg to rad

degreesToRadians<- function (degrees) {
return (degrees * pi / 180)
}

R<-6371e3  #radius of Earth in meters

phi1<-degreesToRadians(lat1) # latitude 1
phi2<-degreesToRadians(lat2) # latitude 2
lambda1<-degreesToRadians(lon1) # longitude 1
lambda2<-degreesToRadians(lon2) # longitude 2

delta_phi<-phi1-phi2 # latitude-distance
delta_lambda<-lambda1-lambda2 # longitude-distance

a<-sin(delta_phi/2)*sin(delta_phi/2)+
cos(phi1)*cos(phi2)*sin(delta_lambda/2)*
sin(delta_lambda/2)

cc<-2*atan2(sqrt(a),sqrt(1-a))

distance<- R * cc

return(distance)  # in meters
}
shghm
fuente
0

Aquí hay una variación de Kotlin:

import kotlin.math.*

class HaversineAlgorithm {

    companion object {
        private const val MEAN_EARTH_RADIUS = 6371.0
        private const val D2R = Math.PI / 180.0
    }

    private fun haversineInKm(lat1: Double, lon1: Double, lat2: Double, lon2: Double): Double {
        val lonDiff = (lon2 - lon1) * D2R
        val latDiff = (lat2 - lat1) * D2R
        val latSin = sin(latDiff / 2.0)
        val lonSin = sin(lonDiff / 2.0)
        val a = latSin * latSin + (cos(lat1 * D2R) * cos(lat2 * D2R) * lonSin * lonSin)
        val c = 2.0 * atan2(sqrt(a), sqrt(1.0 - a))
        return EQATORIAL_EARTH_RADIUS * c
    }
}
Csaba Toth
fuente
¿Por qué usaste el radio ecuatorial en lugar del radio medio de la Tierra?
user13044086
@ user13044086 Buena pregunta. Es porque deduje esto de la versión Java de Paulo Miguel Almeida. Parece que la versión C # también está usando esa distancia. Otras versiones aquí tienen 6371, pero luego debes darte cuenta de que todos estos algoritmos pueden no manejar perfectamente la forma geoide de la Tierra. Siéntase libre de modificar esto y usar 6371. Si me dice que eso conduce a valores más precisos, cambiaré mi respuesta.
Csaba Toth
1
6371.008 se usa comúnmente porque minimiza el error relativo de la fórmula como se explica en las notas en la página movable-type.co.uk/scripts/latlong.html#ellipsoid
user13044086
¡Lo suficientemente justo! Editaré mi respuesta mañana
Csaba Toth
@ user13044086 Gracias por el enlace, edité mi respuesta hace un tiempo basado en eso
Csaba Toth