Bereken de afstand tussen twee breedte- en lengtegraadpunten? (Haversine-formule)

Hoe bereken ik de afstand tussen twee punten gespecificeerd door breedte- en lengtegraad?

Ter verduidelijking, ik wil graag de afstand in kilometers; de punten gebruiken het WGS84-systeem en ik zou graag de relatieve nauwkeurigheid van de beschikbare benaderingen willen begrijpen.


Antwoord 1, autoriteit 100%

Deze linkkan nuttig voor u zijn, omdat deze de gebruik van de Haversine-formuleom de afstand te berekenen.

Uittreksel:

Dit script [in Javascript] berekent grootcirkelafstanden tussen de twee punten
dat wil zeggen, de kortste afstand over het aardoppervlak met behulp van de
Haversine -formule.

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

Antwoord 2, autoriteit 34%

Ik moest veel afstanden tussen de punten berekenen voor mijn project, dus ik ging door en probeerde de code te optimaliseren, die ik hier heb gevonden. Mijn nieuwe implementatie werkt gemiddeld in verschillende browsers twee keer snellerdan het meest gewaardeerde antwoord.

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
}

Je kunt met mijn jsPerf spelen en de resultaten hierbekijken.

Onlangs moest ik hetzelfde doen in python, dus hier is een python-implementatie:

from math import cos, asin, sqrt, pi
def distance(lat1, lon1, lat2, lon2):
    p = 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...

En voor de volledigheid: Haversineop wiki.


Antwoord 3, autoriteit 6%

Hier is een C#-implementatie:

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

Antwoord 4, autoriteit 5%

Hier is een Java-implementatie van de Haversine-formule.

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

Houd er rekening mee dat we het antwoord hier afronden op de dichtstbijzijnde km.


Antwoord 5, autoriteit 3%

Heel erg bedankt voor dit alles. Ik heb de volgende code gebruikt in mijn Objective-C iPhone-app:

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

Breedtegraad en lengtegraad zijn in decimalen. Ik heb min() niet gebruikt voor de asin()-aanroep omdat de afstanden die ik gebruik zo klein zijn dat ze het niet nodig hebben.

Het gaf onjuiste antwoorden totdat ik de waarden in Radians doorgaf – nu is het vrijwel hetzelfde als de waarden verkregen uit de Apple Map-app 🙂

Extra update:

Als je iOS4 of hoger gebruikt, biedt Apple een aantal methoden om dit te doen, zodat dezelfde functionaliteit wordt bereikt met:

-(double)kilometresBetweenPlace1:(CLLocationCoordinate2D) place1 andPlace2:(CLLocationCoordinate2D) place2 {
    MKMapPoint  start, finish;
    start = MKMapPointForCoordinate(place1);
    finish = MKMapPointForCoordinate(place2);
    return MKMetersBetweenMapPoints(start, finish) / 1000;
}

Antwoord 6, autoriteit 3%

Dit is een eenvoudige PHP-functie die een zeer redelijke benadering geeft (minder dan +/-1% foutmarge).

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

Zoals eerder gezegd; de aarde is GEEN bol. Het is als een oude, oude honkbal waarmee Mark McGwire besloot te oefenen – het zit vol met deuken en stoten. De eenvoudigere berekeningen (zoals deze) behandelen het als een bol.

Verschillende methoden kunnen min of meer nauwkeurig zijn, afhankelijk van waar u zich op deze onregelmatige eivorm bevindt EN hoe ver uw punten uit elkaar liggen (hoe dichter ze zijn, hoe kleiner de absolute foutmarge). Hoe nauwkeuriger uw verwachting, hoe complexer de wiskunde.

Voor meer info: wikipedia geografische afstand


Antwoord 7, autoriteit 2%

Ik post hier mijn werkvoorbeeld.

Maak een lijst van alle punten in de tabel met een afstand tussen een aangewezen punt (we gebruiken een willekeurig punt – lat:45.20327, long:23.7806) van minder dan 50 KM, met latitude & longitude, in MySQL (de tabelvelden zijn coord_lat en coord_long):

Lijst van alle met AFSTAND<50, in kilometers (beschouwd als straal van de aarde 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

Het bovenstaande voorbeeld is getest in MySQL 5.0.95 en 5.5.16 (Linux).


Antwoord 8, autoriteit 2%

In de andere antwoorden een implementatie in ontbreekt.

Het berekenen van de afstand tussen twee punten is vrij eenvoudig met de functie distmuit het pakket geosphere:

distm(p1, p2, fun = distHaversine)

waar:

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

Aangezien de aarde niet perfect bolvormig is, is de Vincenty-formule voor ellipsoïdenwaarschijnlijk de beste manier om afstanden te berekenen. Dus in het geospherepakket dat je dan gebruikt:

distm(p1, p2, fun = distVincentyEllipsoid)

Natuurlijk hoef je niet per se het geospherepakket te gebruiken, je kunt ook de afstand in basis Rberekenen met een functie:

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

Antwoord 9

De haversine is zeker een goede formule voor waarschijnlijk de meeste gevallen, andere antwoorden bevatten het al, dus ik ga de ruimte niet innemen. Maar het is belangrijk op te merken dat het niet uitmaakt welke formule wordt gebruikt (ja niet slechts één). Vanwege het enorme nauwkeurigheidsbereik en de benodigde rekentijd. De keuze van de formule vereist wat meer denkwerk dan een eenvoudig antwoord.

Dit bericht van een persoon bij nasa, is het beste dat ik heb gevonden bij het bespreken van de opties

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

Bijvoorbeeld als u rijen sorteert op afstand in een straal van 100 mijl. De platte aarde formule zal veel sneller zijn dan de haversine.

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;

Merk op dat er slechts één cosinus en één vierkantswortel is. Vs 9 van hen op de Haversine-formule.


Antwoord 10

Er zou een eenvoudigere en correctere oplossing kunnen zijn: de omtrek van de aarde is 40.000 km op de evenaar, ongeveer 37.000 op de cyclus van Greenwich (of welke lengtegraad dan ook). Dus:

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

Ik ben het ermee eens dat het moet worden afgesteld, omdat ik zelf zei dat het een ellipsoïde is, dus de straal die met de cosinus moet worden vermenigvuldigd, varieert. Maar het is iets nauwkeuriger. Vergeleken met Google Maps en het heeft de fout aanzienlijk verminderd.


Antwoord 11

pip install haversine

Python-implementatie

Oorsprong is het centrum van de aangrenzende Verenigde Staten.

from haversine import haversine, Unit
origin = (39.50, 98.35)
paris = (48.8567, 2.3508)
haversine(origin, paris, unit=Unit.MILES)

Om het antwoord in kilometers te krijgen, stelt u eenvoudig unit=Unit.KILOMETERSin (dat is de standaardinstelling).


Antwoord 12

Ik vind het niet leuk om nog een antwoord toe te voegen, maar de Google maps API v.3 heeft sferische geometrie (en meer). Na het converteren van uw WGS84 naar decimale graden kunt u dit doen:

<script src="https://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));

Geen woord over hoe nauwkeurig de berekeningen van Google zijn of zelfs welk model wordt gebruikt (hoewel er wel ‘sferisch’ staat in plaats van ‘geoïde’. Trouwens, de ‘rechte lijn’-afstand zal duidelijk verschillen van de afstand als die er is. reist over het aardoppervlak en dat is wat iedereen lijkt te veronderstellen.


Antwoord 13

Alle bovenstaande antwoorden gaan ervan uit dat de aarde een bol is. Een nauwkeurigere benadering zou echter die van een afgeplatte sferoïde zijn.

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

Antwoord 14

Je kunt de build in CLLocationDistance gebruiken om dit te berekenen:

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

In jouw geval, als je kilometers wilt, deel je gewoon door 1000.


Antwoord 15

Hier is een typescriptimplementatie van de Haversine-formule

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

Antwoord 16

Zoals gezegd, moet bij een nauwkeurige berekening rekening worden gehouden met het feit dat de aarde geen perfecte bol is. Hier zijn enkele vergelijkingen van de verschillende algoritmen die hier worden aangeboden:

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

Over kleine afstanden lijkt het algoritme van Keerthana samen te vallen met dat van Google Maps. Google Maps lijkt geen enkel eenvoudig algoritme te volgen, wat suggereert dat dit hier misschien de meest nauwkeurige methode is.

Hoe dan ook, hier is een Javascript-implementatie van het algoritme van 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));
}

Antwoord 17

Hier is de SQL-implementatie om de afstand in km te berekenen,

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;

Voor meer details over de implementatie door de taal te programmeren, kun je gewoon het php-script doorlopen dat hier


Antwoord 18

Dit script [in PHP] berekent de afstanden tussen de twee punten.

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

Antwoord 19

Java-implementatie volgens Haversine-formule

double calculateDistance(double latPoint1, double lngPoint1, 
                         double latPoint2, double lngPoint2) {
    if(latPoint1 == latPoint2 && lngPoint1 == lngPoint2) {
        return 0d;
    }
    final double EARTH_RADIUS = 6371.0; //km value;
    //converting to radians
    latPoint1 = Math.toRadians(latPoint1);
    lngPoint1 = Math.toRadians(lngPoint1);
    latPoint2 = Math.toRadians(latPoint2);
    lngPoint2 = Math.toRadians(lngPoint2);
    double distance = Math.pow(Math.sin((latPoint2 - latPoint1) / 2.0), 2) 
            + Math.cos(latPoint1) * Math.cos(latPoint2)
            * Math.pow(Math.sin((lngPoint2 - lngPoint1) / 2.0), 2);
    distance = 2.0 * EARTH_RADIUS * Math.asin(Math.sqrt(distance));
    return distance; //km value
}

Antwoord 20

Om de afstand tussen twee punten op een bol te berekenen, moet je de Grote Cirkel-berekening.

Er zijn een aantal C/C++-bibliotheken om te helpen met kaartprojectie op MapToolsals dat nodig is projecteer uw afstanden opnieuw op een plat oppervlak. Hiervoor heeft u de projectiestring van de verschillende coördinatenstelsels nodig.

Misschien vindt u MapWindowook een handig hulpmiddel om de punten te visualiseren. Ook als open source is het een handige gids voor het gebruik van de proj.dll-bibliotheek, die de belangrijkste open source-projectiebibliotheek lijkt te zijn.


Antwoord 21

Hier is mijn Java-implementatie voor het berekenen van afstand via decimale graden na wat zoeken. Ik gebruikte de gemiddelde straal van de wereld (van wikipedia) in km. Als u resultaatmijlen wilt, gebruik dan de wereldradius in mijlen.

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

Antwoord 22

Hier is de geaccepteerde antwoordimplementatie geport naar Java voor het geval iemand het nodig heeft.

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

Antwoord 23

hier is een voorbeeld in postgressql (in km, voor miles-versie, vervang 1.609344 door 0.8684-versie)

CREATE OR REPLACE FUNCTION public.geodistance(alat float, alng float, blat  
float, blng  float)
  RETURNS float AS
$BODY$
DECLARE
    v_distance float;
BEGIN
    v_distance = asin( sqrt(
            sin(radians(blat-alat)/2)^2 
                + (
                    (sin(radians(blng-alng)/2)^2) *
                    cos(radians(alat)) *
                    cos(radians(blat))
                )
          )
        ) * cast('7926.3352' as float) * cast('1.609344' as float) ;
    RETURN v_distance;
END 
$BODY$
language plpgsql VOLATILE SECURITY DEFINER;
alter function geodistance(alat float, alng float, blat float, blng float)
owner to postgres;

Antwoord 24

Voor wie op zoek is naar een Excel-formule op basis van WGS-84 & GRS-80-normen:

=ACOS(COS(RADIANS(90-Lat1))*COS(RADIANS(90-Lat2))+SIN(RADIANS(90-Lat1))*SIN(RADIANS(90-Lat2))*COS(RADIANS(Long1-Long2)))*6371

Bron


Antwoord 25

Hier is de implementatie VB.NET, deze implementatie geeft je het resultaat in KM of Miles op basis van een Enum-waarde die je doorgeeft.

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

Antwoord 26

Ik heb de berekening gecomprimeerd door de formule te vereenvoudigen.

Hier is het in Ruby:

include Math
earth_radius_mi = 3959
radians = lambda { |deg| deg * PI / 180 }
coord_radians = lambda { |c| { :lat => radians[c[:lat]], :lng => radians[c[:lng]] } }
# from/to = { :lat => (latitude_in_degrees), :lng => (longitude_in_degrees) }
def haversine_distance(from, to)
  from, to = coord_radians[from], coord_radians[to]
  cosines_product = cos(to[:lat]) * cos(from[:lat]) * cos(from[:lng] - to[:lng])
  sines_product = sin(to[:lat]) * sin(from[:lat])
  return earth_radius_mi * acos(cosines_product + sines_product)
end

Antwoord 27

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

Chuck’s oplossing, ook geldig voor mijlen.


Antwoord 28

Gebruik in Mysql de volgende functie en geef de parameters door als 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;

Antwoord 29

Hier is nog een geconverteerde naar Rubycode:

include Math
#Note: from/to = [lat, long]
def get_distance_in_km(from, to)
  radians = lambda { |deg| deg * Math.PI / 180 }
  radius = 6371 # Radius of the earth in kilometer
  dLat = radians[to[0]-from[0]]
  dLon = radians[to[1]-from[1]]
  cosines_product = Math.sin(dLat/2) * Math.sin(dLat/2) + Math.cos(radians[from[0]]) * Math.cos(radians[to[1]]) * Math.sin(dLon/2) * Math.sin(dLon/2)
  c = 2 * Math.atan2(Math.sqrt(cosines_product), Math.sqrt(1-cosines_product)) 
  return radius * c # Distance in kilometer
end

Antwoord 30

Omdat dit de meest populaire discussie over het onderwerp is, zal ik hier mijn ervaring van eind 2019-begin 2020 toevoegen. Om aan de bestaande antwoorden toe te voegen: mijn focus was het vinden van een nauwkeurige EN snelle (d.w.z. gevectoriseerde) oplossing.

Laten we beginnen met wat hier het meest wordt gebruikt door antwoorden: de Haversine-aanpak. Het is triviaal om te vectoriseren, zie het voorbeeld in python hieronder:

def haversine(lat1, lon1, lat2, lon2):
    """
    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees)
    All args must be of equal length.
    Distances are in meters.
    Ref:
    https://stackoverflow.com/questions/29545704/fast-haversine-approximation-python-pandas
    https://ipython.readthedocs.io/en/stable/interactive/magics.html
    """
    Radius = 6.371e6
    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2
    c = 2 * np.arcsin(np.sqrt(a))
    s12 = Radius * c
    # initial azimuth in degrees
    y = np.sin(lon2-lon1) * np.cos(lat2)
    x = np.cos(lat1)*np.sin(lat2) - np.sin(lat1)*np.cos(lat2)*np.cos(dlon)
    azi1 = np.arctan2(y, x)*180./math.pi
    return {'s12':s12, 'azi1': azi1}

Nauwkeurigheid is het minst nauwkeurig. Wikipedia vermeldt gemiddeld 0,5% relatieve afwijking zonder bronnen. Mijn experimenten laten een kleinere afwijking zien. Hieronder is de vergelijking uitgevoerd op 100.000 willekeurige punten versus mijn bibliotheek, die tot op de millimeter nauwkeurig zou moeten zijn:

np.random.seed(42)
lats1 = np.random.uniform(-90,90,100000)
lons1 = np.random.uniform(-180,180,100000)
lats2 = np.random.uniform(-90,90,100000)
lons2 = np.random.uniform(-180,180,100000)
r1 = inverse(lats1, lons1, lats2, lons2)
r2 = haversine(lats1, lons1, lats2, lons2)
print("Max absolute error: {:4.2f}m".format(np.max(r1['s12']-r2['s12'])))
print("Mean absolute error: {:4.2f}m".format(np.mean(r1['s12']-r2['s12'])))
print("Max relative error: {:4.2f}%".format(np.max((r2['s12']/r1['s12']-1)*100)))
print("Mean relative error: {:4.2f}%".format(np.mean((r2['s12']/r1['s12']-1)*100)))

Uitvoer:

Max absolute error: 26671.47m
Mean absolute error: -2499.84m
Max relative error: 0.55%
Mean relative error: -0.02%

Dus een gemiddelde afwijking van 2,5 km op 100.000 willekeurige coördinatenparen, wat in de meeste gevallen goed kan zijn.

De volgende optie zijn de formules van Vincenty die tot op de millimeter nauwkeurig is, afhankelijk van de convergentiecriteria en die ook kunnen worden gevectoriseerd. Het heeft wel het probleem met convergentie in de buurt van antipodale punten. U kunt het op die punten laten convergeren door de convergentiecriteria te versoepelen, maar de nauwkeurigheid daalt tot 0,25% en meer. Buiten antipodale punten zal Vincenty resultaten leveren die dicht bij Geographiclib liggen met een relatieve fout van gemiddeld minder dan 1.e-6.

Geographiclib, hier genoemd, is echt de huidige gouden standaard. Het heeft verschillende implementaties en is redelijk snel, vooral als je de C++-versie gebruikt.

Als je van plan bent om Python te gebruiken voor iets boven de 10.000 punten, raad ik je aan om mijn gevectoriseerde implementatie te overwegen. Ik heb een geovectorslib-bibliotheek gemaakt met gevectoriseerde Vincenty-routine voor mijn eigen behoeften, die Geographiclib gebruikt als fallback voor bijna antipodale punten. Hieronder vindt u de vergelijking met Geographiclib voor 100.000 punten. Zoals u kunt zien, biedt het tot 20x verbetering voor inverse en 100x voor directemethoden voor 100.000 punten en de kloof zal met het aantal punten toenemen. Nauwkeurigheidsgewijs zal het binnen 1.e-5 rtol van Geographiclib liggen.

Direct method for 100,000 points
94.9 ms ± 25 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
9.79 s ± 1.4 s per loop (mean ± std. dev. of 7 runs, 1 loop each)
Inverse method for 100,000 points
1.5 s ± 504 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
24.2 s ± 3.91 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

LEAVE A REPLY

Please enter your comment!
Please enter your name here

eight − 1 =

Other episodes