Рассчитать расстояние между двумя координатами GPS

Question

Как рассчитать расстояние между двумя координатами GPS (с использованием широты и долготы)?

Answer 1

Calculate the distance between two coordinates by latitude and longitude, включая реализацию Javascript.

Запад и Южные местоположения отрицательные. Помните минуты и секунды из 60, поэтому S31 30 'составляет -31,50 градусов.

Не забывайте, что конвертирует градусы в радианы. У многих языков есть эта функция. Или его простой расчет: 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; 
} 

Вот некоторые примеры использования:

> distanceInKmBetweenCoordinates(0,0,0,0) // Distance between same points should be 0 
0 
> distanceInKmBetweenCoordinates(51.5, 0, 38.8, -77.1) // From London to Arlington 
5918.185064088764 

Answer 2

Я думаю, вы хотите его вдоль искривления земли. Ваши две точки и центр земли находятся на плоскости. Центр Земли - это центр круга на этой плоскости, и две точки (примерно) расположены по периметру этого круга. Из этого вы можете рассчитать расстояние, узнав, какой угол от одной точки к другой.

Если точки не совпадают с высотой, или если вам нужно учитывать, что земля не идеальная сфера, она становится немного сложнее.

Answer 3

Этот код Lua заимствован из материала найденного в Википедии и в GPSbabel инструмент Роберта Липе в:

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 

Answer 4

Недавно я должен был сделать тоже самое. Я нашел this сайт, чтобы быть очень полезным объяснением сферического триггера примерами, с которыми легко было следовать.

Answer 5

Это зависит от того, насколько точно вам это нужно, если вам нужна точная точность, лучше всего взглянуть на алгоритм с использованием эллипсоида, а не сферы, такой как алгоритм Винцентия, который точно соответствует миллиметру. http://en.wikipedia.org/wiki/Vincenty%27s_algorithm

Answer 6

Это очень легко сделать с типом географии в 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 является SRID для WGS84 elipsoidal модели Земли

Answer 7

Look для гаверсинуса с Google; вот мое решение:

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

Answer 8

Здесь находится в C# (лат и долго в радианах):

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

Если ваш лат и долго в градусах затем разделить на 180/PI для преобразования в радианы ,

Answer 9

Функция T-SQL, который я использую для выбора записей по расстоянию для центра

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 

Answer 10

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

Answer 11

// Может быть, ошибка опечатка?
У нас есть неиспользуемая переменная dlon в GetDirection,
Я полагаю

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

должен быть

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

Answer 12

Если вам нужно что-то более точным, то есть look at this.

формулы Vincenty являются две связанные итерационные методы, используемые в геодезии для вычисления расстояния между двумя точками на поверхности сфероида, разработанные Таддеус Винсенти (1975a) Они основаны на предположении, что фигура Земля является сплюснутым сфероидом, и поэтому являются более точными, чем методы, такие как расстояние большого круга , которые принимают сферическую Землю.

Первый (прямой) метод вычисляет местоположение точки, которая представляет собой заданное расстояние и азимут (направление) от другой точки. Второй метод (обратный) вычисляет географическое расстояние и азимут между двумя заданными точками. Они широко используются в геодезии , потому что они точны с точностью до 0,5 мм (0,020 ") на Земле. эллипсоид.

Answer 13

вы можете найти реализацию этого (с некоторым хорошим объяснением) в F # на fssnip

здесь важные части:

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

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

    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 

Answer 14

C# Версия гаверсинуса

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

Here's a .NET Fiddle of this, так что вы можете протестируйте его со своим собственным Lat/Longs.

Answer 15

Это версия от «Генри Vilinskiy», адаптированной для MySQL и Километры:

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; 

Answer 16

мне нужно реализовать это в PowerShell, надеюсь, что это может помочь кому-то еще. Некоторые замечания об этом методе

1. Не расколоть любой из линий или вычисление будет неправильно
2. Для расчета в КМ удалить * 1000 в расчете $ расстояние
3. Изменение $ earthsRadius = 3963.19059 и удалить * 1000 в расчете $ расстояния, на котором на calulate расстояния в милях

4. Я использую гаверсинус, как и другие сообщения указали формулы Vincenty является гораздо более точным

   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 
   } 
   

Answer 17

Java версия Haversine алгоритма на основе римского Makarov`s ответа на эту тему

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

} 

Answer 18

Вот функция Haversine в Python, который я использую:

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} 

Answer 19

I. Что касается "сухарей" метода

1. радиуса Земли отличается на разных Lat. Это необходимо учитывать в алгоритме Хаверсина.
2. Рассмотрит Bearing изменения, которое превращает прямые линии в арках (которые больше)
3. принимая изменений скорости во внимание превратят арку к спиралям (которые длиннее или короче, чем арки)
4. изменения высоты превратят плоские спирали 3D-спирали (которые длиннее). Это очень важно для холмистой местности.

Ниже смотрите функцию в C, который принимает # 1 и # 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. Существует более простой способ, который дает неплохие результаты.

Средняя скорость.

Trip_distance = Trip_average_speed * Trip_time

Так как GPS-Скорость определяется по эффекту Доплера и не имеет прямого отношения к [Lon, Lat] он может быть, по крайней мере рассматривать как вторичное (резервное копирование или исправление), если не в качестве основного расстояния метод расчета.

Answer 20

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

Answer 21

PHP версия:

(. Удалить все deg2rad() если ваши координаты уже в радианах)

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

Answer 22

мне нужно вычислить много расстояний между баллов за мой проект, поэтому я пошел вперед и попытался оптимизировать код, который я нашел здесь. В среднем в разных браузерах моя новая реализация работает в 2 раза быстрее, чем самый лучший ответ.

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 
} 

Вы можете играть с моим JSPerf и увидеть results here.

Недавно мне нужно сделать то же самое в питоне, так вот реализация по питона:

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

И для полноты картины: Haversine на вики.

Answer 23

Если вы используете .NET, не отвлекайте колесо. См. System.Device.Location. Отметьте fnx в комментариях в another answer.

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

Answer 24

здесь Свифт реализация от ответа

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 
} 

Answer 25

я взял верхний ответ и использовал его в программе 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 
} 

я кэрри функцию для того, чтобы быть в состоянии чтобы легко создавать функции, которые имеют одно из двух мест , фиксированное и для получения расстояния требуется только пара lat/lon.