计算2个GPS坐标之间的距离
我如何计算两个GPS坐标之间的距离（使用经度和纬度）？
根据纬度和经度计算两个坐标之间的距离 ，包括一个Javascript实现。
西部和南部地区是负面的。 请记住分钟和秒钟是在60以外，所以S31 30'是-31.50度。
不要忘记将度数转换为弧度 。 许多语言都有这个function。 或者简单的计算： radians = degrees * PI / 180 。
function degreesToRadians(degrees) { return degrees * Math.PI / 180; } function distanceInKmBetweenEarthCoordinates(lat1, lon1, lat2, lon2) { var earthRadiusKm = 6371; var dLat = degreesToRadians(lat2-lat1); var dLon = degreesToRadians(lon2-lon1); lat1 = degreesToRadians(lat1); lat2 = degreesToRadians(lat2); var a = Math.sin(dLat/2) * Math.sin(dLat/2) + Math.sin(dLon/2) * Math.sin(dLon/2) * Math.cos(lat1) * Math.cos(lat2); var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); return earthRadiusKm * c; }
这里有一些用法的例子：
> distanceInKmBetweenCoordinates(0,0,0,0) // Distance between same points should be 0 0 > distanceInKmBetweenCoordinates(51.5, 0, 38.8, -77.1) // From London to Arlington 5918.185064088764
用Googlesearch海峡; 这里是我的解决scheme：
#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; }
Haversine的C＃版本
double _eQuatorialEarthRadius = 6378.1370D; double _d2r = (Math.PI / 180D); private int HaversineInM(double lat1, double long1, double lat2, double long2) { return (int)(1000D * HaversineInKM(lat1, long1, lat2, long2)); } private double HaversineInKM(double lat1, double long1, double lat2, double long2) { double dlong = (long2 - long1) * _d2r; double dlat = (lat2 - lat1) * _d2r; double a = Math.Pow(Math.Sin(dlat / 2D), 2D) + Math.Cos(lat1 * _d2r) * Math.Cos(lat2 * _d2r) * Math.Pow(Math.Sin(dlong / 2D), 2D); double c = 2D * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1D - a)); double d = _eQuatorialEarthRadius * c; return d; }
这是一个.NET的小提琴 ，所以你可以用自己的纬度/经度testing它。
基于Roman Makarov的Haversinealgorithm的Java版本回复此主题
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; } }
在SQL Server 2008中使用地理types很容易。
SELECT geography::Point(lat1, lon1, 4326).STDistance(geography::Point(lat2, lon2, 4326)) -- computes distance in meters using eliptical model, accurate to the mm
4326是WGS84 elipsoidal地球模型的SRID
这取决于你需要的准确度，如果你需要准确的精确度，最好看一个使用椭球的algorithm，而不是像Vincentyalgorithm这样精确到毫米的球体。 http://en.wikipedia.org/wiki/Vincenty%27s_algorithm
这是在C＃中（经度和弧度）：
double CalculateGreatCircleDistance(double lat1, double long1, double lat2, double long2, double radius) { return radius * Math.Acos( Math.Sin(lat1) * Math.Sin(lat2) + Math.Cos(lat1) * Math.Cos(lat2) * Math.Cos(long2 - long1)); }
如果你的经度和纬度都是度数，那么用180 / PI除以弧度。
下面是我使用的Python中的Haversine函数：
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}
PHP版本：
（如果你的坐标已经是弧度，去掉所有的deg2rad() 。）
$R = 6371; // km $dLat = deg2rad($lat2-$lat1); $dLon = deg2rad($lon2-$lon1); $lat1 = deg2rad($lat1); $lat2 = deg2rad($lat2); $a = sin($dLat/2) * sin($dLat/2) + sin($dLon/2) * sin($dLon/2) * cos($lat1) * cos($lat2); $c = 2 * atan2(sqrt($a), sqrt(1-$a)); $d = $R * $c;
一个T-SQL函数，用于按距离为中心selectlogging
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
如果你需要更准确的东西，那么看看这个 。
Vincenty公式是大地测量中用来计算球体表面上两点间距的两种相关的迭代方法，由Thaddeus Vincenty（1975a）开发。他们基于假定地球的形状是扁球体，因此比诸如假定球形地球的大圆距离的方法更精确。
第一种（直接）方法计算距离另一点的给定距离和方位（方向）的点的位置。 第二种方法计算两个给定点之间的地理距离和方位angular。 它们已被广泛用于大地测量，因为它们精确到地球椭球体的0.5毫米（0.020“）以内。
一，关于“面包屑”的方法
- 不同纬度的地球半径是不同的。 Haversinealgorithm必须考虑到这一点。
- 考虑轴承的变化，把直线变成拱形（更长）
- 考虑到速度的变化会使拱门变成螺旋状（比拱门更长或更短）
- 高度变化将变成平坦的螺旋到三维螺旋（再次更长）。 这对丘陵地区非常重要。
下面看看C中考虑＃1和＃2的函数：
double calcDistanceByHaversine(double rLat1, double rLon1, double rHeading1, double rLat2, double rLon2, double rHeading2){ double rDLatRad = 0.0; double rDLonRad = 0.0; double rLat1Rad = 0.0; double rLat2Rad = 0.0; double a = 0.0; double c = 0.0; double rResult = 0.0; double rEarthRadius = 0.0; double rDHeading = 0.0; double rDHeadingRad = 0.0; if ((rLat1 < -90.0) || (rLat1 > 90.0) || (rLat2 < -90.0) || (rLat2 > 90.0) || (rLon1 < -180.0) || (rLon1 > 180.0) || (rLon2 < -180.0) || (rLon2 > 180.0)) { return -1; }; rDLatRad = (rLat2 - rLat1) * DEGREE_TO_RADIANS; rDLonRad = (rLon2 - rLon1) * DEGREE_TO_RADIANS; rLat1Rad = rLat1 * DEGREE_TO_RADIANS; rLat2Rad = rLat2 * DEGREE_TO_RADIANS; a = sin(rDLatRad / 2) * sin(rDLatRad / 2) + sin(rDLonRad / 2) * sin( rDLonRad / 2) * cos(rLat1Rad) * cos(rLat2Rad); if (a == 0.0) { return 0.0; } c = 2 * atan2(sqrt(a), sqrt(1 - a)); rEarthRadius = 6378.1370 - (21.3847 * 90.0 / ((fabs(rLat1) + fabs(rLat2)) / 2.0)); rResult = rEarthRadius * c; // Chord to Arc Correction based on Heading changes. Important for routes with many turns and U-turns if ((rHeading1 >= 0.0) && (rHeading1 < 360.0) && (rHeading2 >= 0.0) && (rHeading2 < 360.0)) { rDHeading = fabs(rHeading1 - rHeading2); if (rDHeading > 180.0) { rDHeading -= 180.0; } rDHeadingRad = rDHeading * DEGREE_TO_RADIANS; if (rDHeading > 5.0) { rResult = rResult * (rDHeadingRad / (2.0 * sin(rDHeadingRad / 2))); } else { rResult = rResult / cos(rDHeadingRad); } } return rResult; }
II。 有一个更简单的方法可以给出相当好的结果。
平均速度。
Trip_distance = Trip_average_speed * Trip_time
由于GPS速度是由多普勒效应检测到的，并且与[Lon，Lat]没有直接关系，所以如果不是作为主距离计算方法，它可以至less被认为是次要的（备份或校正）。
我需要为我的项目计算点之间的很多距离，所以我继续尝试优化代码，我已经在这里find了。 平均而言，在不同的浏览器中，我的新实现的运行速度比最上面的答案快2倍 。
function distance(lat1, lon1, lat2, lon2) { var p = 0.017453292519943295; // Math.PI / 180 var c = Math.cos; var a = 0.5 - c((lat2 - lat1) * p)/2 + c(lat1 * p) * c(lat2 * p) * (1 - c((lon2 - lon1) * p))/2; return 12742 * Math.asin(Math.sqrt(a)); // 2 * R; R = 6371 km }
你可以玩我的jsPerf并在这里看到结果 。
最近我需要在python中做同样的事情，所以这是一个python实现 ：
from math import cos, asin, sqrt def distance(lat1, lon1, lat2, lon2): p = 0.017453292519943295 a = 0.5 - cos((lat2 - lat1) * p)/2 + cos(lat1 * p) * cos(lat2 * p) * (1 - cos((lon2 - lon1) * p)) / 2 return 12742 * asin(sqrt(a))
为了完整起见：维基上的Haversine 。
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); }
这个Lua代码是从维基百科和Robert Lipe的GPSbabel工具中find的东西改编的：
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
我最近不得不做同样的事情。 我发现这个网站是很有帮助的解释球形触发与易于遵循的例子。
你可以在fssnip上find这个实现（有一些很好的解释）
这里是重要的部分：
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
这是适用于MySQL和公里的“Henry Vilinskiy”版本：
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;
我需要在PowerShell中实现这一点，希望它可以帮助别人。 关于这种方法的一些说明
- 不要分割任何一行或者计算错误
- 在KM中计算，在$ distance的计算中删除* 1000
- 更改$ earthsRadius = 3963.19059，并在$ distance的计算中移除* 1000来计算以英里为单位的距离
-
我正在使用Haversine，正如其他文章指出，Vincenty的公式更准确
Function MetresDistanceBetweenTwoGPSCoordinates($latitude1, $longitude1, $latitude2, $longitude2) { $Rad = ([math]::PI / 180); $earthsRadius = 6378.1370 # Earth's Radius in KM $dLat = ($latitude2 - $latitude1) * $Rad $dLon = ($longitude2 - $longitude1) * $Rad $latitude1 = $latitude1 * $Rad $latitude2 = $latitude2 * $Rad $a = [math]::Sin($dLat / 2) * [math]::Sin($dLat / 2) + [math]::Sin($dLon / 2) * [math]::Sin($dLon / 2) * [math]::Cos($latitude1) * [math]::Cos($latitude2) $c = 2 * [math]::ATan2([math]::Sqrt($a), [math]::Sqrt(1-$a)) $distance = [math]::Round($earthsRadius * $c * 1000, 0) #Multiple by 1000 to get metres Return $distance }
我想你想要沿着地球的曲率。 你的两点和地球的中心在飞机上。 地球的中心是该平面上一个圆的中心，两个点（大致）在该圆的周边。 从这里你可以通过找出从一点到另一点的angular度来计算距离。
如果分数不是相同的高度，或者如果你需要考虑到地球不是一个完美的球体，它会变得更加困难。
//可能是拼写错误？
我们在GetDirection中有一个未使用的variablesdlon，
我假设
double y = Math.Sin(dlon) * Math.Cos(lat2); // cannot use degrees in Cos ?
应该
double y = Math.Sin(dlon) * Math.Cos(dlat);
斯卡拉版本
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) ) }
如果您使用.NET，请不要转动方向盘。 请参阅System.Device.Location 。 在另一个答案中的评论中信用fnx。
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);
这里是答案的Swift实现
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 }
