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

Эта ссылка может быть полезна вам, так как она подробно описывает использование формулы Хаверсина для расчета расстояния.

Выписка:

Этот скрипт [in Javascript] вычисляет расстояния большого круга между двумя точками – то есть кратчайшее расстояние по земной поверхности-с помощью ‘Гаверсинус формула.

function getDistanceFromLatLonInKm(lat1,lon1,lat2,lon2) {
  var R = 6371; // Radius of the earth in km
  var dLat = deg2rad(lat2-lat1);  // deg2rad below
  var dLon = deg2rad(lon2-lon1); 
  var a = 
    Math.sin(dLat/2) * Math.sin(dLat/2) +
    Math.cos(deg2rad(lat1)) * Math.cos(deg2rad(lat2)) * 
    Math.sin(dLon/2) * Math.sin(dLon/2)
    ; 
  var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); 
  var d = R * c; // Distance in km
  return d;
}

function deg2rad(deg) {
  return deg * (Math.PI/180)
}

Мне нужно было рассчитать большое количество расстояний между точками для моего проекта, поэтому я пошел дальше и попытался оптимизировать код, который я нашел здесь. В среднем в разных браузерах моя новая реализация работает в 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
}

Вы можете поиграть с my jsPerf и посмотреть результаты здесь .

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

from math import cos, asin, sqrt
def distance(lat1, lon1, lat2, lon2):
    p = 0.017453292519943295     #Pi/180
    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)) #2*R*asin...

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

Вот это реализация C# :

static class DistanceAlgorithm
{
    const double PIx = 3.141592653589793;
    const double RADIUS = 6378.16;

    /// <summary>
    /// Convert degrees to Radians
    /// </summary>
    /// <param name="x">Degrees</param>
    /// <returns>The equivalent in radians</returns>
    public static double Radians(double x)
    {
        return x * PIx / 180;
    }

    /// <summary>
    /// Calculate the distance between two places.
    /// </summary>
    /// <param name="lon1"></param>
    /// <param name="lat1"></param>
    /// <param name="lon2"></param>
    /// <param name="lat2"></param>
    /// <returns></returns>
    public static double DistanceBetweenPlaces(
        double lon1,
        double lat1,
        double lon2,
        double lat2)
    {
        double dlon = Radians(lon2 - lon1);
        double dlat = Radians(lat2 - lat1);

        double a = (Math.Sin(dlat / 2) * Math.Sin(dlat / 2)) + Math.Cos(Radians(lat1)) * Math.Cos(Radians(lat2)) * (Math.Sin(dlon / 2) * Math.Sin(dlon / 2));
        double angle = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
        return angle * RADIUS;
    }

}

Вот java реализация формулы Хаверсина.

public final static double AVERAGE_RADIUS_OF_EARTH_KM = 6371;
public int calculateDistanceInKilometer(double userLat, double userLng,
  double venueLat, double venueLng) {

    double latDistance = Math.toRadians(userLat - venueLat);
    double lngDistance = Math.toRadians(userLng - venueLng);

    double a = Math.sin(latDistance / 2) * Math.sin(latDistance / 2)
      + Math.cos(Math.toRadians(userLat)) * Math.cos(Math.toRadians(venueLat))
      * Math.sin(lngDistance / 2) * Math.sin(lngDistance / 2);

    double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));

    return (int) (Math.round(AVERAGE_RADIUS_OF_EARTH_KM * c));
}

Обратите внимание, что здесь мы округляем ответ до ближайшего километра.

Большое спасибо за все это. Я использовал следующий код в своем приложении Objective-C iPhone:

const double PIx = 3.141592653589793;
const double RADIO = 6371; // Mean radius of Earth in Km

double convertToRadians(double val) {

   return val * PIx / 180;
}

-(double)kilometresBetweenPlace1:(CLLocationCoordinate2D) place1 andPlace2:(CLLocationCoordinate2D) place2 {

        double dlon = convertToRadians(place2.longitude - place1.longitude);
        double dlat = convertToRadians(place2.latitude - place1.latitude);

        double a = ( pow(sin(dlat / 2), 2) + cos(convertToRadians(place1.latitude))) * cos(convertToRadians(place2.latitude)) * pow(sin(dlon / 2), 2);
        double angle = 2 * asin(sqrt(a));

        return angle * RADIO;
}

Широта и долгота указаны в десятичной системе счисления. Я не использовал min() для вызова asin(), поскольку расстояния, которые я использую, настолько малы, что они не требуют этого.

Он давал неверные ответы, пока я не передал значения в радианах - теперь это почти то же самое, что значения, полученные из приложения Карты Apple :-)

Дополнительное обновление:

Если вы используете iOS4 или более позднюю версию, то Apple предоставляет некоторые методы для этого, чтобы та же функциональность была достигнута с помощью:

-(double)kilometresBetweenPlace1:(CLLocationCoordinate2D) place1 andPlace2:(CLLocationCoordinate2D) place2 {

    MKMapPoint  start, finish;


    start = MKMapPointForCoordinate(place1);
    finish = MKMapPointForCoordinate(place2);

    return MKMetersBetweenMapPoints(start, finish) / 1000;
}

Это простая функция PHP, которая даст очень разумное приближение (при допустимой погрешности +/-1%).

<?php
function distance($lat1, $lon1, $lat2, $lon2) {

    $pi80 = M_PI / 180;
    $lat1 *= $pi80;
    $lon1 *= $pi80;
    $lat2 *= $pi80;
    $lon2 *= $pi80;

    $r = 6372.797; // mean radius of Earth in km
    $dlat = $lat2 - $lat1;
    $dlon = $lon2 - $lon1;
    $a = sin($dlat / 2) * sin($dlat / 2) + cos($lat1) * cos($lat2) * sin($dlon / 2) * sin($dlon / 2);
    $c = 2 * atan2(sqrt($a), sqrt(1 - $a));
    $km = $r * $c;

    //echo '<br/>'.$km;
    return $km;
}
?>

Как уже было сказано, земля-это сфера. Это похоже на старый - престарый бейсбол, с которым Марк McGwire решил попрактиковаться-он полон вмятин и ударов. Более простые вычисления (как это) рассматривают его как сферу.

Различные методы могут быть более или менее точными в зависимости от того, где вы находитесь на этом неправильном яйце AND, как далеко друг от друга находятся ваши точки (чем ближе они находятся, тем меньше абсолютная погрешность). Чем точнее ваши ожидания, тем сложнее математика.

Для получения дополнительной информации: Википедия географическое расстояние

Я публикую здесь свой рабочий пример.

Перечислите все точки в таблице, имеющие расстояние между обозначенной точкой (мы используем случайную точку-lat:45.20327, long:23.7806) меньше 50 KM, с широтой & долготы, в MySQL (поля таблицы coord_lat и coord_long):

Перечислите все имеющие DISTANCE<50, в километрах (считается радиус Земли 6371 KM):

SELECT denumire, (6371 * acos( cos( radians(45.20327) ) * cos( radians( coord_lat ) ) * cos( radians( 23.7806 ) - radians(coord_long) ) + sin( radians(45.20327) ) * sin( radians(coord_lat) ) )) AS distanta 
FROM obiective 
WHERE coord_lat<>'' 
    AND coord_long<>'' 
HAVING distanta<50 
ORDER BY distanta desc

Приведенный выше пример был протестирован в MySQL 5.0.95 и 5.5.16 (Linux).

В других ответах отсутствует реализация в r .

Вычисление расстояния между двумя точками довольно просто с помощью функции distm из пакета geosphere :

distm(p1, p2, fun = distHaversine)

где:

p1 = longitude/latitude for point(s)
p2 = longitude/latitude for point(s)
# type of distance calculation
fun = distCosine / distHaversine / distVincentySphere / distVincentyEllipsoid 

Поскольку земля не является идеально сферической, формула Винсенти для эллипсоидов , вероятно, является лучшим способом вычисления расстояний. Таким образом, в пакете geosphere вы используете затем:

distm(p1, p2, fun = distVincentyEllipsoid)

Конечно, вам не обязательно использовать пакет geosphere , вы также можете рассчитать расстояние в базе R с помощью функции:

hav.dist <- function(long1, lat1, long2, lat2) {
  R <- 6371
  diff.long <- (long2 - long1)
  diff.lat <- (lat2 - lat1)
  a <- sin(diff.lat/2)^2 + cos(lat1) * cos(lat2) * sin(diff.long/2)^2
  b <- 2 * asin(pmin(1, sqrt(a))) 
  d = R * b
  return(d)
}

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

Это сообщение от человека из nasa, это лучшее, что я нашел при обсуждении вариантов

http://www.cs.nyu.edu/visual/home/proj/tiger/gisfaq.html

Например, если вы просто сортируете строки по расстоянию в радиусе 100 миль. Формула плоской Земли будет гораздо быстрее, чем хаверсин.

HalfPi = 1.5707963;
R = 3956; /* the radius gives you the measurement unit*/

a = HalfPi - latoriginrad;
b = HalfPi - latdestrad;
u = a * a + b * b;
v = - 2 * a * b * cos(longdestrad - longoriginrad);
c = sqrt(abs(u + v));
return R * c;

Обратите внимание, что существует только один Косинус и один квадратный корень. Против 9 из них на формула Гаверсинуса.

Вы можете использовать сборку в CLLocationDistance для расчета этого:

CLLocation *location1 = [[CLLocation alloc] initWithLatitude:latitude1 longitude:longitude1];
CLLocation *location2 = [[CLLocation alloc] initWithLatitude:latitude2 longitude:longitude2];
[self distanceInMetersFromLocation:location1 toLocation:location2]

- (int)distanceInMetersFromLocation:(CLLocation*)location1 toLocation:(CLLocation*)location2 {
    CLLocationDistance distanceInMeters = [location1 distanceFromLocation:location2];
    return distanceInMeters;
}

В вашем случае, если вы хотите километры просто разделить на 1000.

Мне не нравится добавлять еще один ответ, но Google maps API v.3 имеет сферическую геометрию (и многое другое). После преобразования вашего WGS84 в десятичные Градусы вы можете сделать это:

<script src="http://maps.google.com/maps/api/js?sensor=false&libraries=geometry" type="text/javascript"></script>  

distance = google.maps.geometry.spherical.computeDistanceBetween(
    new google.maps.LatLng(fromLat, fromLng), 
    new google.maps.LatLng(toLat, toLng));

Ни слова о том, насколько точны расчеты Google или даже какая модель используется (хотя она говорит "spherical", а не "geoid". Кстати, расстояние "straight line", очевидно, будет отличаться от расстояния, если человек путешествует по поверхности Земли, что, по-видимому, все и предполагают.

Реализация Python Происхождение-это центр сопредельных Соединенных Штатов.

from haversine import haversine
origin = (39.50, 98.35)
paris = (48.8567, 2.3508)
haversine(origin, paris, miles=True)

Чтобы получить ответ в километрах, просто установите miles=false.

Могло бы быть и более простое решение, и более правильное: периметр Земли составляет 40 000 км на экваторе, около 37 000 на Гринвичском (или любом долготном) цикле. Таким образом:

pythagoras = function (lat1, lon1, lat2, lon2) {
   function sqr(x) {return x * x;}
   function cosDeg(x) {return Math.cos(x * Math.PI / 180.0);}

   var earthCyclePerimeter = 40000000.0 * cosDeg((lat1 + lat2) / 2.0);
   var dx = (lon1 - lon2) * earthCyclePerimeter / 360.0;
   var dy = 37000000.0 * (lat1 - lat2) / 360.0;

   return Math.sqrt(sqr(dx) + sqr(dy));
};

Я согласен, что он должен быть точно настроен, так как я сам сказал, что это эллипсоид, поэтому радиус, который нужно умножить на косинус, меняется. Но это немного более точно. По сравнению с Google Maps и это действительно значительно уменьшило ошибку.

Все вышеприведенные ответы предполагают, что Земля-это сфера. Однако более точным приближением был бы сплюснутый сфероид.

a= 6378.137#equitorial radius in km
b= 6356.752#polar radius in km

def Distance(lat1, lons1, lat2, lons2):
    lat1=math.radians(lat1)
    lons1=math.radians(lons1)
    R1=(((((a**2)*math.cos(lat1))**2)+(((b**2)*math.sin(lat1))**2))/((a*math.cos(lat1))**2+(b*math.sin(lat1))**2))**0.5 #radius of earth at lat1
    x1=R*math.cos(lat1)*math.cos(lons1)
    y1=R*math.cos(lat1)*math.sin(lons1)
    z1=R*math.sin(lat1)

    lat2=math.radians(lat2)
    lons2=math.radians(lons2)
    R1=(((((a**2)*math.cos(lat2))**2)+(((b**2)*math.sin(lat2))**2))/((a*math.cos(lat2))**2+(b*math.sin(lat2))**2))**0.5 #radius of earth at lat2
    x2=R*math.cos(lat2)*math.cos(lons2)
    y2=R*math.cos(lat2)*math.sin(lons2)
    z2=R*math.sin(lat2)

    return ((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)**0.5

Это скорее зависит от того, насколько точным вы хотите быть и на каком датуме определяются lat и long. Очень, очень приблизительно вы делаете небольшую сферическую тригонометрию, но поправка на то, что Земля не является сферой, делает формулы более сложными.

Вот это реализация typescript формула Гаверсинуса

static getDistanceFromLatLonInKm(lat1: number, lon1: number, lat2: number, lon2: number): number {
    var deg2Rad = deg => {
        return deg * Math.PI / 180;
    }

    var r = 6371; // Radius of the earth in km
    var dLat = deg2Rad(lat2 - lat1);   
    var dLon = deg2Rad(lon2 - lon1);
    var a =
        Math.sin(dLat / 2) * Math.sin(dLat / 2) +
        Math.cos(deg2Rad(lat1)) * Math.cos(deg2Rad(lat2)) *
        Math.sin(dLon / 2) * Math.sin(dLon / 2);
    var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
    var d = r * c; // Distance in km
    return d;
}

Этот скрипт [in PHP] вычисляет расстояния между двумя точками.

public static function getDistanceOfTwoPoints($source, $dest, $unit='K') {
        $lat1 = $source[0];
        $lon1 = $source[1];
        $lat2 = $dest[0];
        $lon2 = $dest[1];

        $theta = $lon1 - $lon2;
        $dist = sin(deg2rad($lat1)) * sin(deg2rad($lat2)) +  cos(deg2rad($lat1)) * cos(deg2rad($lat2)) * cos(deg2rad($theta));
        $dist = acos($dist);
        $dist = rad2deg($dist);
        $miles = $dist * 60 * 1.1515;
        $unit = strtoupper($unit);

        if ($unit == "K") {
            return ($miles * 1.609344);
        }
        else if ($unit == "M")
        {
            return ($miles * 1.609344 * 1000);
        }
        else if ($unit == "N") {
            return ($miles * 0.8684);
        } 
        else {
            return $miles;
        }
    }

Вот реализация SQL для расчета расстояния в км,

SELECT UserId, ( 3959 * acos( cos( radians( your latitude here ) ) * cos( radians(latitude) ) * 
cos( radians(longitude) - radians( your longitude here ) ) + sin( radians( your latitude here ) ) * 
sin( radians(latitude) ) ) ) AS distance FROM user HAVING
distance < 5  ORDER BY distance LIMIT 0 , 5;

Чтобы вычислить расстояние между двумя точками на сфере, вам нужно сделать расчет Большого Круга .

Существует ряд библиотек C/C++, которые помогут с картографической проекцией на MapTools , если вам нужно перепроецировать свои расстояния на плоскую поверхность. Для этого вам понадобится проекционная строка из различных систем координат.

Вы также можете найти MapWindow полезным инструментом для визуализации точек. Кроме того, как его открытый исходный код это полезное руководство по использованию библиотеки proj.dll, которая, по-видимому, является основной библиотекой проекции с открытым исходным кодом.

Вот общепринятая реализация ответа, перенесенная на Java, если кому-то это нужно.

package com.project529.garage.util;


/**
 * Mean radius.
 */
private static double EARTH_RADIUS = 6371;

/**
 * Returns the distance between two sets of latitudes and longitudes in meters.
 * <p/>
 * Based from the following JavaScript SO answer:
 * http://stackoverflow.com/questions/27928/calculate-distance-between-two-latitude-longitude-points-haversine-formula,
 * which is based on https://en.wikipedia.org/wiki/Haversine_formula (error rate: ~0.55%).
 */
public double getDistanceBetween(double lat1, double lon1, double lat2, double lon2) {
    double dLat = toRadians(lat2 - lat1);
    double dLon = toRadians(lon2 - lon1);

    double a = Math.sin(dLat / 2) * Math.sin(dLat / 2) +
            Math.cos(toRadians(lat1)) * Math.cos(toRadians(lat2)) *
                    Math.sin(dLon / 2) * Math.sin(dLon / 2);
    double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
    double d = EARTH_RADIUS * c;

    return d;
}

public double toRadians(double degrees) {
    return degrees * (Math.PI / 180);
}

Как уже отмечалось, при точном расчете следует учитывать, что Земля не является совершенной сферой. Вот некоторые сравнения различных алгоритмов, предлагаемых здесь:

geoDistance(50,5,58,3)
Haversine: 899 km
Maymenn: 833 km
Keerthana: 897 km
google.maps.geometry.spherical.computeDistanceBetween(): 900 km

geoDistance(50,5,-58,-3)
Haversine: 12030 km
Maymenn: 11135 km
Keerthana: 10310 km
google.maps.geometry.spherical.computeDistanceBetween(): 12044 km

geoDistance(.05,.005,.058,.003)
Haversine: 0.9169 km
Maymenn: 0.851723 km
Keerthana: 0.917964 km
google.maps.geometry.spherical.computeDistanceBetween(): 0.917964 km

geoDistance(.05,80,.058,80.3)
Haversine: 33.37 km
Maymenn: 33.34 km
Keerthana: 33.40767 km
google.maps.geometry.spherical.computeDistanceBetween(): 33.40770 km

На небольших расстояниях алгоритм Киртаны, похоже, совпадает с алгоритмом Google Maps. Google Maps, похоже, не следует какому-либо простому алгоритму, предполагая, что это может быть самый точный метод здесь.

Во всяком случае, здесь-это реализация Javascript алгоритма Keerthana по :

function geoDistance(lat1, lng1, lat2, lng2){
    const a = 6378.137; // equitorial radius in km
    const b = 6356.752; // polar radius in km

    var sq = x => (x*x);
    var sqr = x => Math.sqrt(x);
    var cos = x => Math.cos(x);
    var sin = x => Math.sin(x);
    var radius = lat => sqr((sq(a*a*cos(lat))+sq(b*b*sin(lat)))/(sq(a*cos(lat))+sq(b*sin(lat))));

    lat1 = lat1 * Math.PI / 180;
    lng1 = lng1 * Math.PI / 180;
    lat2 = lat2 * Math.PI / 180;
    lng2 = lng2 * Math.PI / 180;

    var R1 = radius(lat1);
    var x1 = R1*cos(lat1)*cos(lng1);
    var y1 = R1*cos(lat1)*sin(lng1);
    var z1 = R1*sin(lat1);

    var R2 = radius(lat2);
    var x2 = R2*cos(lat2)*cos(lng2);
    var y2 = R2*cos(lat2)*sin(lng2);
    var z2 = R2*sin(lat2);

    return sqr(sq(x1-x2)+sq(y1-y2)+sq(z1-z2));
}

Вот реализация VB.NET, эта реализация даст вам результат в KM или милях на основе значения перечисления, которое вы передадите.

Public Enum DistanceType
    Miles
    KiloMeters
End Enum

Public Structure Position
    Public Latitude As Double
    Public Longitude As Double
End Structure

Public Class Haversine

    Public Function Distance(Pos1 As Position,
                             Pos2 As Position,
                             DistType As DistanceType) As Double

        Dim R As Double = If((DistType = DistanceType.Miles), 3960, 6371)

        Dim dLat As Double = Me.toRadian(Pos2.Latitude - Pos1.Latitude)

        Dim dLon As Double = Me.toRadian(Pos2.Longitude - Pos1.Longitude)

        Dim a As Double = Math.Sin(dLat / 2) * Math.Sin(dLat / 2) + Math.Cos(Me.toRadian(Pos1.Latitude)) * Math.Cos(Me.toRadian(Pos2.Latitude)) * Math.Sin(dLon / 2) * Math.Sin(dLon / 2)

        Dim c As Double = 2 * Math.Asin(Math.Min(1, Math.Sqrt(a)))

        Dim result As Double = R * c

        Return result

    End Function

    Private Function toRadian(val As Double) As Double

        Return (Math.PI / 180) * val

    End Function

End Class

function getDistanceFromLatLonInKm(lat1,lon1,lat2,lon2,units) {
  var R = 6371; // Radius of the earth in km
  var dLat = deg2rad(lat2-lat1);  // deg2rad below
  var dLon = deg2rad(lon2-lon1); 
  var a = 
    Math.sin(dLat/2) * Math.sin(dLat/2) +
    Math.cos(deg2rad(lat1)) * Math.cos(deg2rad(lat2)) * 
    Math.sin(dLon/2) * Math.sin(dLon/2)
    ; 
  var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); 
  var d = R * c; 
  var miles = d / 1.609344; 

if ( units == 'km' ) {  
return d; 
 } else {
return miles;
}}

Решение Чака, действительное и для Майлза тоже.

Вот моя java реализация для вычисления расстояния через десятичные Градусы после некоторого поиска. Я использовал средний радиус мира (из Википедии) в километрах. Если вы хотите получить результат миль, то используйте мировой радиус в милях.

public static double distanceLatLong2(double lat1, double lng1, double lat2, double lng2) 
{
  double earthRadius = 6371.0d; // KM: use mile here if you want mile result

  double dLat = toRadian(lat2 - lat1);
  double dLng = toRadian(lng2 - lng1);

  double a = Math.pow(Math.sin(dLat/2), 2)  + 
          Math.cos(toRadian(lat1)) * Math.cos(toRadian(lat2)) * 
          Math.pow(Math.sin(dLng/2), 2);

  double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));

  return earthRadius * c; // returns result kilometers
}

public static double toRadian(double degrees) 
{
  return (degrees * Math.PI) / 180.0d;
}

В Mysql используйте следующую функцию передайте параметры как с помощью POINT(LONG,LAT)

CREATE FUNCTION `distance`(a POINT, b POINT)
 RETURNS double
    DETERMINISTIC
BEGIN

RETURN

GLength( LineString(( PointFromWKB(a)), (PointFromWKB(b)))) * 100000; -- To Make the distance in meters

END;