Calculer la distance entre 2 coordonnées GPS

si vous avez besoin de quelque chose de plus précis alors avoir un regardez ce .

les formules de Vincenty sont deux méthodes itératives connexes utilisées en géodésie pour calculer la distance entre deux points sur la surface d'un sphéroïde, développé par Thaddeus Vincenty (1975a) ils sont basés sur le l'hypothèse que la figure de la Terre est un sphéroïde aplati, et par conséquent sont plus précis que des méthodes telles que la distance de grand-cercle qui supposent une Terre sphérique.

la première méthode (directe) calcule l'emplacement d'un point qui est Distance et azimut donnés (direction) d'un autre point. Deuxième (l'inverse) méthode calcule la distance géographique et l'azimut entre deux points donnés. Ils ont été largement utilisés en géodésie parce qu'ils sont précis à moins de 0.5 mm( 0.020") sur la Terre ellipsoïde.

I. concernant" Breadcrumbs "method

1. le rayon de la Terre est différent selon la Lat. Ceci doit être pris en considération dans Haversine algorithme.
2. Envisager un changement de Roulement, qui tourne, les lignes droites des arcs (qui sont plus)
3. si l'on tient compte du changement de vitesse, Les Arcs se transforment en spirales (qui sont plus longues ou plus courtes que les arcs)
4. changement D'Altitude tournera les spirales plates en spirales 3D (qui sont plus longues encore). C'est très important pour les régions montagneuses.

voir ci-dessous la fonction en C qui prend les numéros 1 et 2 en compte:

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. Il y a une façon plus facile qui donne de très bons résultats.

Par Vitesse Moyenne.

Trip_distance = Trip_average_speed * Trip_time

depuis la vitesse GPS est détecté par effet Doppler et n'est pas directement lié à [Lon,Lat] il peut être au moins considéré comme secondaire (sauvegarde ou correction) si ce n'est pas comme méthode de calcul de la distance principale.

    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);
    }

C'est la version de" Henry Vilinskiy "adaptée pour MySQL et kilomètres:

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;

si vous utilisez .NET ne ré-enfoncez pas la roue. Voir Système .Appareil.Lieu . Crédit à la fnx dans les commentaires dans une autre réponse .

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);

je suppose que vous voulez le long de la courbure de la terre. Vos deux points et le centre de la terre sont dans un avion. Le centre de la terre est le centre d'un cercle sur l'avion et les deux points sont (à peu près) sur le périmètre de ce cercle. De ce que vous pouvez calculer la distance par trouver ce que l'angle d'un point à l'autre.

si les points ne sont pas les mêmes hauteurs, ou si vous devez tenir compte que la terre n'est pas une parfaite sphère, c'est un peu plus difficile.

ce code Lua est adapté d'éléments trouvés sur Wikipedia et dans GPSBabel "de Robert Lipe outil:

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

j'ai récemment dû faire la même chose. J'ai trouvé ce site web pour être très utile expliquant trig sphérique avec des exemples qui étaient faciles à suivre avec.

vous pouvez trouver une mise en œuvre de ce (avec une bonne explication) dans F# sur fssnip

voici les parties importantes:

let GreatCircleDistance&lt[&ltMeasure&gt] 'u&gt (R : float&lt'u&gt) (p1 : Location) (p2 : Location) =
    let degToRad (x : float&ltdeg&gt) = System.Math.PI * x / 180.0&ltdeg/rad&gt

    let sq x = x * x
    // take the sin of the half and square the result
    let sinSqHf (a : float&ltrad&gt) = (System.Math.Sin &gt&gt sq) (a / 2.0&ltrad&gt)
    let cos (a : float&ltdeg&gt) = System.Math.Cos (degToRad a / 1.0&ltrad&gt)

    let dLat = (p2.Latitude - p1.Latitude) |&gt degToRad
    let dLon = (p2.Longitude - p1.Longitude) |&gt 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

j'avais besoin de mettre en œuvre ceci dans PowerShell, j'espère que cela pourra aider quelqu'un d'autre. Quelques notes sur cette méthode

1. ne divisez aucune des lignes ou le calcul sera erroné
2. pour calculer en KM supprimer le * 1000 dans le calcul de la distance $
3. Variation $earthsRadius = 3963.19059 et supprimer * 1000 dans le calcul de $distance la de calculez la distance en miles

4. J'utilise Haversine, comme d'autres Billets ont fait remarquer que les formules de Vincenty sont beaucoup plus précises

   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  
   }
   

voici L'implémentation Swift de la réponse

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
}

j'ai pris la réponse supérieure et l'ai utilisé dans un programme 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
}

je curry la fonction, afin d'être en mesure de produire facilement des fonctions qui ont l'un des deux emplacements fixe et ne nécessitent qu'une paire de lat/lon pour produire de la distance.

/ / peut-être une erreur typographique ?

Nous avons une variable non utilisée dlon dans GetDirection,

Je suppose

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

devrait être

double y = Math.Sin(dlon) * Math.Cos(dlat);

version 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)
    )
  }