Oblicz odległość między 2 współrzędnymi GPS


Odpowiedzi:


406

Oblicz odległość między dwiema współrzędnymi według szerokości i długości geograficznej , w tym implementacji Javascript.

Zachód i południowe są ujemne. Pamiętaj, że minuty i sekundy są poza 60, więc S31 30 'wynosi -31,50 stopnia.

Nie zapomnij przekonwertować stopni na radiany . Wiele języków ma tę funkcję. Lub jego kalkulacja prosta 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;
}

Oto kilka przykładów użycia:

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

17
W przypadku, gdy nie jest to oczywiste, metoda toRad () jest dostosowywanie się do numeru prototypu, takich jak: Number.prototype.toRad = function() { return this * (Math.PI / 180); }; . Lub, jak wskazano poniżej, możesz zastąpić (Math.PI/2)0,0174532925199433 (... dowolną precyzję, którą uznasz za niezbędną) dla zwiększenia wydajności.
Vinney Kelly,

44
Jeśli ktoś, szczególnie ci z was, którzy nie szukają komentarzy na końcu linii, wpatrują się w tę formułę i szukają jednostki odległości, jednostka jest km. :)
Dylan Knowles

1
@ VinneyKelly Mała literówka, ale nie zastępuje (Math.PI / 180) nie (Math.PI / 2), dzięki za pomoc wszystkim
Patrick Murphy

1
@ChristianKRider Spójrz na pierwszą linię. Zastanów się, co Rzwykle oznacza matematyka, a następnie wyszukaj odpowiednie ilości związane z Ziemią, aby sprawdzić, czy liczby się zgadzają.
Pozew funduszu Moniki z

3
W przypadku jednostek imperialnych (mil) możesz zmienić earthRadiusKmna var earthRadiusMiles = 3959;fyi.
chapeljuice

59

Szukaj haversine w Google; oto moje rozwiązanie:

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

3
Możesz zastąpić (M_PI / 180.0) 0,0174532925199433, aby uzyskać lepszą wydajność.
Hlung

3
Pod względem wydajności: można obliczyć sin (dlat / 2.0) tylko raz, zapisać go w zmiennej a1, a zamiast pow (, 2) DUŻO lepiej użyć a1 * a1. To samo dotyczy drugiego pow (, 2).
pms

71
Tak, lub po prostu użyj kompilatora po latach 60.
prawej

17
Nie ma potrzeby „optymalizacji” (M_PI / 180.0) do stałej, której nikt nie rozumie bez kontekstu. Kompilator oblicza te stałe warunki dla Ciebie!
Patrick Cornelissen,

2
@ TõnuSamuel Dziękuję bardzo za komentarz. Bardzo to doceniam. Ma to sens, że kompilator z włączoną optymalizacją (-O) może wstępnie obliczyć operacje stałych, dzięki czemu ręczne zwijanie jest bezużyteczne. Przetestuję to, kiedy będę miał czas.
Hlung

44

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

Oto skrzypce .NET tego rozwiązania , dzięki czemu możesz go przetestować za pomocą własnych Lat / Longs.


1
Dodałem również sprawdzone skrzypce .NET, aby ludzie mogli łatwo to przetestować.
Pure.Krome

7
.Net Framework ma wbudowaną metodę GeoCoordinate.GetDistanceTo. Należy odwołać się do zestawu System.Device. Artykuł MSDN msdn.microsoft.com/en-us/library/…
fnx

27

Wersja Java algorytmu Haversine na podstawie odpowiedzi Romana Makarowa na ten wątek

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

}

@Radu upewnij się, że używasz go poprawnie i nie wymieniasz miejsc lat / log podczas przekazywania ich do dowolnej metody.
Paulo Miguel Almeida,

1
Otrzymałem dość bliską odpowiedź przy użyciu tej formuły. Dokładność oparłem na tej stronie: movable-type.co.uk/scripts/latlong.html co dało mi 0.07149km, podczas gdy twoja formuła dała mi 0.07156dokładność około 99%
Janac Meena

24

Jest to bardzo łatwe do zrobienia z typem geograficznym w 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 jest SRID dla elipsoidalnego modelu Ziemi WGS84


19

Oto funkcja Haversine w Pythonie, której używam:

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}


11

Tutaj jest w C # (lat i long w radianach):

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

Jeśli twoje łat i long są w stopniach, podziel przez 180 / PI, aby zamienić na radiany.


1
Jest to obliczenie „sferycznego prawa cosinusów”, które jest najmniej dokładną i podatną na błędy metodą obliczania dużej odległości koła.
John Machin

11

Musiałem obliczyć wiele odległości między punktami dla mojego projektu, więc poszedłem dalej i spróbowałem zoptymalizować kod, który znalazłem tutaj. Przeciętnie w różnych przeglądarkach moja nowa implementacja działa 2 razy szybciej niż najbardziej pozytywna odpowiedź.

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
}

Możesz grać z moim jsPerf i zobaczyć wyniki tutaj .

Ostatnio musiałem zrobić to samo w Pythonie, więc oto implementacja Pythona :

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

I dla kompletności: Haversine na wiki.


11

Wersja PHP:

(Usuń wszystkie, deg2rad()jeśli współrzędne są już w radianach.)

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

1
Zmień lat1 i lat2 na $ lat1 i $ lat2.
zamiast

7

Funkcja T-SQL, której używam do wybierania rekordów według odległości dla centrum

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

Jest to obliczenie „sferycznego prawa cosinusów”, które jest najmniej dokładną i podatną na błędy metodą obliczania dużej odległości koła.
John Machin

5

Jeśli potrzebujesz czegoś dokładniejszego, spójrz na to .

Wzory Vincenta to dwie powiązane metody iteracyjne stosowane w geodezji do obliczania odległości między dwoma punktami na powierzchni sferoidy, opracowane przez Thaddeusa Vincentego (1975a). Opierają się one na założeniu, że Ziemia jest spłaszczoną sferoidą, a zatem są dokładniejsze niż metody takie jak odległość od wielkiego koła, która zakłada sferyczną Ziemię.

Pierwsza (bezpośrednia) metoda oblicza położenie punktu, który jest daną odległością i azymutem (kierunkiem) od innego punktu. Druga (odwrotna) metoda oblicza odległość geograficzną i azymut między dwoma danymi punktami. Były one szeroko stosowane w geodezji, ponieważ mają dokładność do 0,5 mm (0,020 ″) na elipsoidzie Ziemi.


5

I. W odniesieniu do metody „bułki tartej”

  1. Promień Ziemi jest różny dla różnych Lat. Należy to wziąć pod uwagę w algorytmie Haversine.
  2. Rozważ zmianę łożyska, która zamienia proste linie w łuki (które są dłuższe)
  3. Biorąc pod uwagę zmianę Prędkości, łuki zamieniają się w spirale (dłuższe lub krótsze niż łuki)
  4. Zmiana wysokości zamieni płaskie spirale w 3D (które znów są dłuższe). Jest to bardzo ważne w obszarach pagórkowatych.

Poniżej zobacz funkcję w C, która uwzględnia nr 1 i nr 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. Jest prostszy sposób, który daje całkiem dobre wyniki.

Według średniej prędkości.

Trip_distance = Trip_average_speed * Trip_time

Ponieważ prędkość GPS jest wykrywana przez efekt Dopplera i nie jest bezpośrednio związana z [Lon, Lat], można ją przynajmniej uznać za wtórną (tworzenie kopii zapasowych lub korekcję), jeśli nie jako główną metodę obliczania odległości.


4

Jeśli używasz platformy .NET, nie należy ponownie zakładać koła. Zobacz System.Device.Location . Podziękowania dla fnx w komentarzach w innej odpowiedzi .

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

3

Ten kod Lua został zaadaptowany z materiałów znalezionych w Wikipedii oraz w narzędziu GPSbabel Roberta Lipe'a :

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

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

1
Myślę, że twoja funkcja GetDistance zwraca wartość w metrach
Przemek

Czy to jest poprawne? GetDirection () nie korzysta z „dlat”.
pyzaty

3

To jest wersja z „Henry Vilinskiy” dostosowana do MySQL i Kilometrów:

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;

MySQLpowiedziałSomething is wrong in your syntax near '' on line 8 // declare distance float;
Legionar

Jest to obliczenie „sferycznego prawa cosinusów”, które jest najmniej dokładną i podatną na błędy metodą obliczania odległości dużego koła
John Machin

3

oto implementacja Swift od odpowiedzi

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
}

3

wybrałem najlepszą odpowiedź i użyłem jej w programie 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
}

Curryję tę funkcję, aby móc łatwo tworzyć funkcje, które mają jedną z dwóch ustalonych lokalizacji i wymagają tylko pary lat / lon do wytworzenia odległości.


2

Chyba chcesz tego wzdłuż krzywizny ziemi. Twoje dwa punkty i środek ziemi znajdują się na płaszczyźnie. Środek ziemi jest środkiem koła na tej płaszczyźnie, a dwa punkty znajdują się (z grubsza) na obwodzie tego koła. Na tej podstawie możesz obliczyć odległość, sprawdzając, jaki jest kąt od jednego punktu do drugiego.

Jeśli punkty nie są tej samej wysokości lub jeśli musisz wziąć pod uwagę, że ziemia nie jest idealną kulą, staje się ona nieco trudniejsza.


2

Niedawno musiałem zrobić to samo. Okazało się, że ta strona internetowa jest bardzo pomocna w wyjaśnianiu sferycznych wyzwalaczy za pomocą przykładów, które były łatwe do naśladowania.


2

implementację tego (z dobrym wyjaśnieniem) znajdziesz w F # on fssnip

Oto ważne części:


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

2

Musiałem zaimplementować to w PowerShell, mam nadzieję, że może pomóc komuś innemu. Kilka uwag na temat tej metody

  1. Nie dziel żadnej linii, bo obliczenia będą niepoprawne
  2. Aby obliczyć w KM, usuń * 1000 z obliczenia odległości $
  3. Zmień $ earthsRadius = 3963.19059 i usuń * 1000 w obliczeniach $ odległości, aby obliczyć odległość w milach
  4. Używam Haversine, jak zauważyły ​​inne posty, formuły Vincenta są znacznie bardziej dokładne

    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  
    }
    

2

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

1

// Może błąd literowy?
Mamy nieużywaną zmienną dłoń w GetDirection,
Zakładam

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

Powinien być

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

1
To nie jest odpowiedź, to w najlepszym razie komentarz.
Kevin

1

Oto moja implementacja w 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

1

Wersja Dart

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

0

Myślę, że wciąż brakuje wersji algorytmu w języku 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
}

0

Oto wariant Kotlina:

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

Dlaczego zastosowałeś promień równikowy zamiast średniego promienia Ziemi?
user13044086

@ user13044086 Dobre pytanie. To dlatego, że wyprowadziłem to z wersji Java Paulo Miguela Almeidy. Wygląda na to, że wersja C # również korzysta z tej odległości. Inne wersje mają tutaj 6371, ale musisz zdać sobie sprawę, że wszystkie te algorytmy mogą nie radzić sobie doskonale z geoidą Ziemi. Zmodyfikuj to i użyj 6371. Jeśli powiesz mi, że prowadzi to do bardziej precyzyjnych wartości, zmienię moją odpowiedź.
Csaba Toth

1
6371.008 jest powszechnie używany, ponieważ minimalizuje błąd względny formuły, jak wyjaśniono w uwagach na stronie movable-type.co.uk/scripts/latlong.html#ellipsoid
użytkownik13044086

Słusznie! Zmienię swoją odpowiedź jutro
Csaba Toth

@ user13044086 Dzięki za link, edytowałem odpowiedź jakiś czas temu na podstawie tego
Csaba Toth
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.