25 ответов:
вычислить расстояние между двумя координатами по широте и долготе, включая реализацию Javascript.
Запад и Южный расположение отрицательное. Помните, что минуты и секунды находятся вне 60, поэтому S31 30 ' -31.50 градусов.
Не забудьте преобразование градусов в радианы. Многие языки имеют эту функцию. Или его простой расчет:
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) / / расстояние между одинаковыми точками должно быть 0 Ноль distanceInKmBetweenCoordinates(51.5, 0, 38.8, -77.1) // из Лондона в Арлингтоне 5918.185064088764
ищите haversine с Google; Вот мое решение:
#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; }
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; }вот .NET Скрипка этого, так что вы можете проверить его с вашим собственным широта/долгота.
Java-версия алгоритма Haversine на основе ответа Романа Макарова на этот поток
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.
SELECT geography::Point(lat1, lon1, 4326).STDistance(geography::Point(lat2, lon2, 4326)) -- computes distance in meters using eliptical model, accurate to the mm4326-это SRID для модели элипсоидной Земли WGS84
Это зависит от того, насколько точно вам это нужно, если вам нужна точная точность, лучше всего посмотреть на алгоритм с использованием эллипсоида, а не сферы, например, алгоритм Винсенти, который точен до мм.http://en.wikipedia.org/wiki/Vincenty%27s_algorithm
вот функция Haversine в Python, которую я использую:
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}
вот он в C# (lat и long в радианах):
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)); }Если ваши lat и long находятся в градусах, то разделите на 180/PI, чтобы преобразовать в радианы.
версия 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, которую я использую для выбора записей по расстоянию для центра
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
мне нужно было рассчитать много расстояний между точками для моего проекта, поэтому я пошел вперед и попытался оптимизировать код, который я нашел здесь. В среднем в разных браузерах моя новая реализация работает в 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 на вики.
Если вам нужно что-то более точное, то есть смотреть на это.
формулы Винсенти являются двумя связанными итерационными методами, используемыми в геодезии чтобы вычислить расстояние между двумя точками на поверхности сфероид, разработанный Thaddeus Vincenty (1975a) они основаны на предположение, что фигура Земли представляет собой сплюснутый сфероид, и следовательно, более точны, чем такие методы, как расстояние большого круга которые предполагают сферическую Земля.
первый (прямой) метод вычисляет местоположение точки, которая является задается расстояние и Азимут (направление) от другой точки. Второй (обратный) метод вычисляет географическое расстояние и Азимут между двумя заданными точками. Они широко используются в геодезии потому что они с точностью до 0.5 мм (0.020") на Земле эллипсоид.
I. относительно метода "панировочных сухарей"
- радиус Земли отличается на разных лат. Это должно быть учтено в алгоритме Haversine.
- рассмотрим изменение подшипника, которое превращает прямые линии в арки (которые длиннее)
- принимая во внимание изменение скорости превратит арки в спирали (которые длиннее или короче, чем арки)
- изменение высоты превратит плоские спирали в 3D спирали (которые длиннее снова.) Это очень важно для холмистых районов.
ниже смотрите функцию в 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; }второй. Есть более простой способ, который дает довольно хорошие результаты.
По Средней Скорости.
Trip_distance = Trip_average_speed * Trip_time
поскольку скорость GPS определяется эффектом Доплера и не имеет прямого отношения к [Lon, Lat], ее можно по крайней мере рассматривать как вторичную (резервную или коррекция), если не в качестве основного метода расчета расстояния.
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); }
это версия от "Генриха Вилинского" адаптированная для MySQL и километров:
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;
Я думаю, вы хотите его вдоль кривизны Земли. Две точки и центр Земли находятся на плоскости. Центр Земли-это центр окружности на этой плоскости и две точки (примерно) по периметру этого круга. Из этого вы можете рассчитать расстояние, выяснив, какой угол от одной точки до другой.
Если точки не совпадают по высоте, или если вам нужно учитывать, что Земля не является идеальной сферой, она получает немного сложнее.
этот код Lua адаптирован из материалов, найденных в Википедии и в Роберте Липе GPSbabel:
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
недавно мне пришлось сделать то же самое. Я нашел этой сайт, чтобы быть очень полезным, объясняя сферической тригонометрии с примерами, которые легко следовать.
вы можете найти реализацию этого (с некоторым хорошим объяснением) в F# on fssnip
вот важные части:
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
Мне нужно реализовать это в PowerShell, надеюсь, что это может помочь кому-то еще. Некоторые заметки об этом методе
- не разделяйте ни одну из строк или расчет будет неправильным
- для расчета в КМ уберите * 1000 при расчете $ distance
- изменить $earthsRadius = 3963.19059 и удалить * 1000 в расчете $ distance чтобы вычислить расстояние в милях
Я использую Haversine, как и другие сообщения указанная формула Винсенти гораздо более точна
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 }
вот быстрая реализация из ответа
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 }
Я взял верхний ответ и использовал его в программе 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 }Я каррированные функции для того, чтобы иметь возможность легко производить функции, которые имеют одну из двух локаций основные и требуют лишь пары широта/долгота, чтобы произвести на расстоянии.
// может, опечатка ?
У нас есть неиспользуемая переменная dlon в GetDirection,
Я предполагаю, чтоdouble y = Math.Sin(dlon) * Math.Cos(lat2); // cannot use degrees in Cos ?должно быть
double y = Math.Sin(dlon) * Math.Cos(dlat);
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) ) }
Comments