Berechnen Sie die Entfernung zwischen 2 GPS-Koordinaten

Antworten:

406

Berechnen Sie den Abstand zwischen zwei Koordinaten nach Breite und Länge , einschließlich einer Javascript-Implementierung.

West und Süd - Standorte sind negativ. Denken Sie daran, dass Minuten und Sekunden von 60 abweichen, sodass S31 30 '-31,50 Grad beträgt.

Vergessen Sie nicht , Grad in Bogenmaß umzurechnen . Viele Sprachen haben diese Funktion. Oder es ist eine einfache Berechnung : 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;
}

Hier einige Anwendungsbeispiele:

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
quelle
17
Falls es nicht offensichtlich ist, ist die toRad () -Methode eine Anpassung des Number- Prototyps wie : Number.prototype.toRad = function() { return this * (Math.PI / 180); }; . Oder, wie unten angegeben, können Sie es durch (Math.PI/2)0.0174532925199433 (... unabhängig von der Genauigkeit, die Sie für erforderlich halten) ersetzen, um die Leistung zu steigern.
Vinney Kelly
44
Wenn jemand, insbesondere diejenigen unter Ihnen, die nicht nach Kommentaren am Zeilenende suchen, auf diese Formel starrt und nach einer Entfernungseinheit sucht, ist die Einheit km. :)
Dylan Knowles
1
@ VinneyKelly Kleiner Tippfehler, aber ersetzen (Math.PI / 180) nicht (Math.PI / 2), danke für jedermanns Hilfe
Patrick Murphy
1
@ChristianKRider Schau dir die erste Zeile an. Überlegen Sie, was Rnormalerweise in der Mathematik bedeutet, und suchen Sie dann nach relevanten erdbezogenen Größen, um festzustellen, ob die Zahlen übereinstimmen.
Fund Monica Klage
3
Für imperialen Einheiten (Meilen) Sie könnte sich ändern earthRadiusKmwerden var earthRadiusMiles = 3959;, zu ihrer Information.
Kapellensaft
59

Suchen Sie bei Google nach Haversine. Hier ist meine Lösung:

#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
quelle
3
Sie können (M_PI / 180.0) durch 0.0174532925199433 ersetzen, um eine bessere Leistung zu erzielen.
Hlung
3
In Bezug auf die Leistung: Man könnte sin (dlat / 2.0) nur einmal berechnen, in der Variablen a1 speichern und anstelle von pow (, 2) ist es VIEL besser, a1 * a1 zu verwenden. Gleiches gilt für die andere Leistung (, 2).
PMS
71
Ja, oder verwenden Sie einfach einen Compiler nach den 60ern.
Rechtsfalte
17
Es ist nicht erforderlich, (M_PI / 180.0) auf eine Konstante zu "optimieren", die niemand ohne Kontext versteht. Der Compiler berechnet diese festen Bedingungen für Sie!
Patrick Cornelissen
2
@ TõnuSamuel Vielen Dank für Ihren Kommentar. Ich weiß das wirklich zu schätzen. Es ist sinnvoll, dass der Compiler mit aktivierter Optimierung (-O) Operationen von Konstanten vorberechnen kann, wodurch das manuelle Reduzieren unbrauchbar wird. Ich werde es testen, wenn ich Zeit habe.
Hlung
44

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

Hier ist eine .NET-Geige davon , damit Sie sie mit Ihren eigenen Lat / Longs testen können.

Roman Makarov
quelle
1
Ich habe auch eine karierte .NET-Geige hinzugefügt, damit die Leute dies leicht testen können.
Pure.Krome
7
Das .Net Framework verfügt über eine integrierte Methode GeoCoordinate.GetDistanceTo. Auf die Assembly System.Device muss verwiesen werden. MSDN Artikel msdn.microsoft.com/en-us/library/…
fnx
27

Java-Version des Haversine-Algorithmus basierend auf der Antwort von Roman Makarov auf diesen Thread

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
quelle
@Radu stellen Sie sicher, dass Sie es richtig verwenden und keine Lat / Log-Stellen austauschen, wenn Sie sie an eine Methode übergeben.
Paulo Miguel Almeida
1
Mit dieser Formel habe ich eine ziemlich genaue Antwort erhalten. Ich habe die Genauigkeit anhand dieser Website ermittelt: movable-type.co.uk/scripts/latlong.html, die mir 0.07149km gab, während Ihre Formel mir 0.07156eine Genauigkeit von ungefähr 99% gab
Janac Meena
24

Dies ist mit dem Geografietyp in SQL Server 2008 sehr einfach.

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

4326 ist SRID für das elipsoidale Erdmodell WGS84

Marko Tintor
quelle
19

Hier ist eine Haversine-Funktion in Python, die ich verwende:

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
quelle
16

Es hängt davon ab, wie genau Sie es benötigen. Wenn Sie eine genaue Genauigkeit benötigen, ist es am besten, einen Algorithmus zu betrachten, bei dem ein Ellipsoid anstelle einer Kugel verwendet wird, z. B. der Algorithmus von Vincenty, der auf mm genau ist. http://en.wikipedia.org/wiki/Vincenty%27s_algorithm

Seanb
quelle
11

Hier ist es in C # (lat und long im Bogenmaß):

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

Wenn Lat und Long in Grad angegeben sind, dividieren Sie durch 180 / PI, um sie im Bogenmaß umzurechnen.

Mike Chamberlain
quelle
1
Dies ist die Berechnung des "sphärischen Kosinusgesetzes", die die am wenigsten genaue und fehleranfälligste Methode zur Berechnung eines Großkreisabstands ist.
John Machin
11

Ich musste viele Entfernungen zwischen den Punkten für mein Projekt berechnen, also habe ich versucht, den Code zu optimieren, den ich hier gefunden habe. Im Durchschnitt läuft meine neue Implementierung in verschiedenen Browsern zweimal schneller als die am besten bewertete Antwort.

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
}

Sie können mit meinem jsPerf spielen und die Ergebnisse hier sehen .

Vor kurzem musste ich dasselbe in Python tun, daher hier eine Python-Implementierung :

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

Und der Vollständigkeit halber: Haversine im Wiki.

Salvador Dali
quelle
11

PHP-Version:

(Entfernen Sie alle, deg2rad()wenn Ihre Koordinaten bereits im Bogenmaß liegen.)

$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
quelle
1
Bitte ändern Sie lat1 und lat2 in $ lat1 und $ lat2.
stattdessen
7

Eine T-SQL-Funktion, mit der ich Datensätze nach Entfernung für ein Zentrum auswähle

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
quelle
Dies ist die Berechnung des "sphärischen Kosinusgesetzes", die die am wenigsten genaue und fehleranfälligste Methode zur Berechnung eines Großkreisabstands ist.
John Machin
5

Wenn Sie etwas genaueres benötigen, schauen Sie sich das an .

Vincentys Formeln sind zwei verwandte iterative Methoden, die in der Geodäsie verwendet werden, um den Abstand zwischen zwei Punkten auf der Oberfläche eines Sphäroids zu berechnen, entwickelt von Thaddeus Vincenty (1975a). Sie basieren auf der Annahme, dass die Figur der Erde ein abgeflachter Sphäroid ist, und daher sind genauer als Methoden wie der Großkreisabstand, die eine sphärische Erde annehmen.

Die erste (direkte) Methode berechnet die Position eines Punktes, der eine bestimmte Entfernung und einen bestimmten Azimut (Richtung) von einem anderen Punkt ist. Die zweite (inverse) Methode berechnet die geografische Entfernung und den Azimut zwischen zwei gegebenen Punkten. Sie sind in der Geodäsie weit verbreitet, da sie auf dem Erdellipsoid auf 0,5 mm genau sind.

dsmelser
quelle
5

I. In Bezug auf die "Breadcrumbs" -Methode

  1. Der Erdradius ist bei verschiedenen Lat unterschiedlich. Dies muss im Haversine-Algorithmus berücksichtigt werden.
  2. Betrachten Sie einen Lagerwechsel, bei dem gerade Linien zu Bögen (die länger sind) werden.
  3. Unter Berücksichtigung der Geschwindigkeitsänderung werden Bögen zu Spiralen (die länger oder kürzer als Bögen sind).
  4. Durch Höhenänderung werden flache Spiralen zu 3D-Spiralen (die wieder länger sind). Dies ist sehr wichtig für hügelige Gebiete.

Unten sehen Sie die Funktion in C, die # 1 und # 2 berücksichtigt:

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. Es gibt einen einfacheren Weg, der ziemlich gute Ergebnisse liefert.

Nach Durchschnittsgeschwindigkeit.

Trip_distance = Trip_average_speed * Trip_time

Da die GPS-Geschwindigkeit durch den Doppler-Effekt erkannt wird und nicht direkt mit [Lon, Lat] zusammenhängt, kann sie zumindest als sekundär (Sicherung oder Korrektur) betrachtet werden, wenn nicht als Methode zur Berechnung der Hauptentfernung.

Tod Samay
quelle
4

Wenn Sie .NET verwenden, aktivieren Sie das Rad nicht erneut. Siehe System.Device.Location . Gutschrift an fnx in den Kommentaren in einer anderen Antwort .

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 Partridge
quelle
3

Dieser Lua-Code wurde aus Wikipedia und dem GPSbabel- Tool von Robert Lipe übernommen :

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
quelle
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.
quelle
1
Ich denke, Ihre Funktion GetDistance gibt Wert in Metern zurück
Przemek
Ist das richtig? GetDirection () verwendet 'dlat' nicht.
Gubby
3

Dies ist eine Version von "Henry Vilinskiy", angepasst für MySQL und Kilometer:

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
quelle
MySQLsagteSomething is wrong in your syntax near '' on line 8 // declare distance float;
Legionär
Dies ist die Berechnung des "sphärischen Gesetzes der Kosinus", die die am wenigsten genaue und fehleranfälligste Methode zur Berechnung einer Großkreisentfernung ist
John Machin,
3

Hier ist die Swift-Implementierung aus der Antwort

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
quelle
3

Ich nahm die beste Antwort und verwendete sie in einem Scala-Programm

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
}

Ich habe die Funktion aktiviert, um auf einfache Weise Funktionen erstellen zu können, bei denen eine der beiden Positionen festgelegt ist und für die die Entfernung nur ein Paar Lat / Lon benötigt.

Peter Perháč
quelle
2

Ich denke du willst es entlang der Krümmung der Erde. Ihre beiden Punkte und der Erdmittelpunkt liegen auf einer Ebene. Der Erdmittelpunkt ist der Mittelpunkt eines Kreises auf dieser Ebene, und die beiden Punkte befinden sich (ungefähr) am Umfang dieses Kreises. Daraus können Sie die Entfernung berechnen, indem Sie den Winkel von einem Punkt zum anderen ermitteln.

Wenn die Punkte nicht die gleichen Höhen haben oder wenn Sie berücksichtigen müssen, dass die Erde keine perfekte Kugel ist, wird es etwas schwieriger.

Guge
quelle
2

Ich musste kürzlich das Gleiche tun. Ich fand diese Website sehr hilfreich, um sphärische Trigger anhand von Beispielen zu erklären, denen man leicht folgen konnte.

Twotymz
quelle
2

Eine Implementierung davon (mit einigen guten Erklärungen) finden Sie in F # auf fssnip

Hier sind die wichtigen Teile:


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
quelle
2

Ich musste dies in PowerShell implementieren, hoffe, es kann jemand anderem helfen. Einige Hinweise zu dieser Methode

  1. Teilen Sie keine der Zeilen, sonst ist die Berechnung falsch
  2. Um in KM zu berechnen, entfernen Sie die * 1000 in der Berechnung der $ Entfernung
  3. Ändern Sie $ earthsRadius = 3963.19059 und entfernen Sie * 1000 bei der Berechnung von $ distance the, um die Entfernung in Meilen zu berechnen
  4. Ich benutze Haversine, da andere Beiträge darauf hingewiesen haben, dass Vincentys Formeln viel genauer sind

    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
quelle
2

Scala-Version

  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
quelle
1

// Vielleicht ein Tippfehler? Ich nehme an,
wir haben eine nicht verwendete Variable dlon in GetDirection

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

sollte sein

double y = Math.Sin(dlon) * Math.Cos(dlat);
Gerrie de Jager
quelle
1
Dies ist keine Antwort, es ist bestenfalls ein Kommentar.
Kevin
1

Hier ist meine Implementierung in 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
quelle
1

Dart Version

Haversine-Algorithmus.

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
quelle
0

Ich denke, eine Version des Algorithmus in R fehlt noch:

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
quelle
0

Hier ist eine Kotlin-Variante:

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
quelle
Warum haben Sie den äquatorialen Radius anstelle des mittleren Erdradius verwendet?
user13044086
@ user13044086 Gute Frage. Das liegt daran, dass ich dies aus der Java-Version von Paulo Miguel Almeida abgeleitet habe. Anscheinend verwendet die C # -Version auch diesen Abstand. Andere Versionen hier haben 6371, aber dann müssen Sie erkennen, dass all diese Algorithmen möglicherweise nicht perfekt mit der Geoidform der Erde umgehen. Fühlen Sie sich frei, dies zu ändern und 6371 zu verwenden. Wenn Sie mir sagen, dass dies zu genaueren Werten führt, werde ich meine Antwort ändern.
Csaba Toth
1
6371.008 wird häufig verwendet, da es den relativen Fehler der Formel minimiert, wie in den Hinweisen auf der Seite movable-type.co.uk/scripts/latlong.html#ellipsoid
user13044086
Meinetwegen! Ich werde meine Antwort morgen bearbeiten
Csaba Toth
@ user13044086 Danke für den Link, ich habe meine Antwort vor einiger Zeit basierend darauf bearbeitet
Csaba Toth