Как вычислить расстояние зная широту и долготу?
// Способ первый
type
TDistance = record
Distance: Double; // расстояние в метрах
Bearing: Double; // смещение в градусах
end;
implementation
uses
{...,} Math;
/// Начальная широта – градусы и сотые доли
/// Начальная долгота – градусы и сотые доли
/// Конечная широта – градусы и сотые доли
/// Конечная долгота – градусы и сотые доли
function GetDistance( StartLat, StartLong, EndLat, EndLong: Double ): TDistance;
const
// Константы, используемые для вычисления смещения и расстояния
D2R: Double = 0.017453; // Константа для преобразования градусов в радианы
R2D: Double = 57.295781; // Константа для преобразования радиан в градусы
a: Double = 6378137.0; // Основные полуоси
b: Double = 6356752.314245; // Не основные полуоси
e2: Double = 0.006739496742337; // Квадрат эксцентричности эллипсоида
f: Double = 0.003352810664747; // Выравнивание эллипсоида
var
// Переменные, используемые для вычисления смещения и расстояния
fPhimean: Double; // Средняя широта
fdLambda: Double; // Разница между двумя значениями долготы
fdPhi: Double; // Разница между двумя значениями широты
fAlpha: Double; // Смещение
fRho: Double; // Меридианский радиус кривизны
fNu: Double; // Поперечный радиус кривизны
fR: Double; // Радиус сферы Земли
fz: Double; // Угловое расстояние от центра сфероида
fTemp: Double; // Временная переменная, используемая в вычислениях
Distance: Double; // Вычисленное расстояния в метрах
Bearing: Double; // Вычисленное от и до смещение
begin
// Вычисляем разницу между двумя долготами и широтами и получаем среднюю широту
fdLambda := ( StartLong - EndLong ) * D2R;
fdPhi := ( StartLat - EndLat ) * D2R;
fPhimean := ( StartLat + EndLat ) / 2.0 * D2R;
// Вычисляем меридианные и поперечные радиусы кривизны средней широты
fTemp := 1 - e2 * Power( Sin( fPhimean ), 2 );
fRho := a * ( 1 - e2 ) / Power( fTemp, 1.5 );
fNu := a / Sqrt( 1 - e2 * Sin( fPhimean ) * Sin( fPhimean ) );
// Вычисляем угловое расстояние
fz := Sqrt( Power( Sin( fdPhi/2.0 ), 2 ) + Cos( EndLat * D2R ) *
Cos( StartLat * D2R ) * Power( Sin( fdLambda/2.0 ), 2 ) );
fz := 2 * ArcSin( fz );
// Вычисляем смещение
fAlpha := ArcSin( Cos( EndLat * D2R ) * Sin( fdLambda ) / Sin( fz ) );
// Вычисляем радиус Земли
// fR := fRho * fNu / ( ( fRho * Power( Sin( fAlpha ), 2 ) ) +
// ( fNu * Power( Cos( fAlpha ), 2 ) ) );
fR := fRho * fNu / ( fRho * Power( Sin( fAlpha ), 2 ) + fNu *
Power( Cos( fAlpha ), 2 ) );
// Получаем смещение и расстояние
Result.Distance := fz * fR;
if ( ( StartLat < EndLat ) and ( StartLong < EndLong ) ) then
Result.Bearing := Abs( fAlpha * R2D )
else
if ( ( StartLat < EndLat ) and ( StartLong > EndLong ) ) then
Result.Bearing := 360 - Abs( fAlpha * R2D )
else
if ( ( StartLat > EndLat ) and ( StartLong > EndLong ) ) then
Result.Bearing := 180 + Abs( fAlpha * R2D )
else if ( ( StartLat > EndLat ) and ( StartLong < EndLong ) ) then
Result.Bearing := 180 - Abs( fAlpha * R2D );
end;
procedure TForm1.Button1Click(Sender: TObject);
var
Distance: TDistance;
begin
Distance := GetDistance( 55.756339, 37.623612, 55.758928, 37.626016 );
ShowMessage( FloatToStr( Distance.Distance ) );
end;
// Способ второй [XE8]
uses
{...,} System.Math, FMX.Maps;
function GetDistance( const aStart, aEnd: TMapCoordinate ): Real;
const
EarthRadius = 6372795;
var
CosLatStart, SinLatStart, CosLatEnd, SinLatEnd,
Delta, CosDelta, SinDelta, X, Y: Real;
begin
try
CosLatStart := Cos( DegToRad( aStart.Latitude ) );
CosLatEnd := Cos( DegToRad( aEnd.Latitude ) );
SinLatStart := Sin( DegToRad( aStart.Latitude ) );
SinLatEnd := Sin( DegToRad( aEnd.Latitude ) );
Delta := ( DegToRad( aEnd.Longitude ) ) - ( DegToRad( aStart.Longitude ) );
CosDelta := Cos( Delta );
SinDelta := Sin( Delta );
Y := Sqrt( ( ( CosLatEnd * SinDelta ) * ( CosLatEnd * SinDelta ) ) +
( ( CosLatStart * SinLatEnd - SinLatStart * CosLatEnd * CosDelta ) *
( CosLatStart * SinLatEnd - SinLatStart * CosLatEnd * CosDelta ) ) );
X := SinLatStart * SinLatEnd + CosLatStart * CosLatEnd * CosDelta;
Result := Round( ArcTan2( Y, X ) * EarthRadius );
except
Result := -1;
end;
end;
procedure TForm1.Button1Click(Sender: TObject);
var
aStart, aEnd: TMapCoordinate;
begin
aStart.Latitude := 55.756339;
aStart.Longitude := 37.623612;
aEnd.Latitude := 55.758928;
aEnd.Longitude := 37.626016;
ShowMessage( FloatToStr( GetDistance( aStart, aEnd ) ) );
end;
// Способ третий
// Для модели эллипсоида WGS84
uses
{...,} System.Math;
function GetDistance(Lat1, Lon1, Lat2, Lon2: Double): Double;
var
fdLambda: Double;
fdPhi: Double;
fPhimean: Double;
fTemp: Double;
fRho: Double;
fNu: Double;
fz: Double;
fAlpha: Double;
fR: Double;
const
e2 = 0.0067394967423;
a = 6378137;
begin
Result := 0;
if (Abs(Lon1-Lon2) <= 0.0000014) and (Abs(Lat1-Lat2) <= 0.0000008) then
Exit
else
begin
Lat1 := Lat1 * Pi / 180;
Lon1 := Lon1 * Pi / 180;
Lat2 := Lat2 * Pi / 180;
Lon2 := Lon2 * Pi / 180;
fdLambda := Lon1 - Lon2;
fdPhi := Lat1 - Lat2;
fPhimean := (Lat1 + Lat2) / 2;
fTemp := 1 - e2 * Power(Sin(fPhimean), 2);
fRho := a * (1 - e2) / Power(fTemp, 1.5);
fNu := a / Sqrt(1 - e2 * Sin(fPhimean) * Sin(fPhimean));
fz := 2 * ArcSin(Sqrt(Power(Sin(fdPhi / 2.0), 2) + Cos(Lat2) *
Cos(Lat1) * Power(Sin(fdLambda / 2.0), 2)));
fAlpha := ArcSin(Cos(Lat2) * Sin(fdLambda) / Sin(fz));
fR := fRho * fNu / (fRho * Power(Sin(fAlpha), 2) + fNu * Power(Cos(fAlpha), 2));
Result := fz * fR;
end;
end;
procedure TForm1.Button3Click(Sender: TObject);
begin
ShowMessage(FloatToStr(GetDistance(55.756339, 37.623612, 55.758928, 37.626016)));
end;
// Способ четвертый
uses
{...,} System.Math;
type
TPointF = record
X: Double;
Y: Double;
end;
function GetGeoDistance(P1, P2: TPointF): Double;
var
x1, x2, y1, y2: Double;
begin
try
x1 := P1.X * pi / 180;
y1 := P1.Y * pi / 180;
x2 := P2.X * pi / 180;
y2 := P2.Y * pi / 180;
Result := 6378137 * 2 * ArcSin(Sqrt(Sqr(Sin((y1 - y2) / 2)) +
Cos(y1) * Cos(y2) * Sqr(Sin((x1 - x2) / 2)))) ;
except
Result := 0;
end;
end;
procedure TForm1.Button4Click(Sender: TObject);
var
p1, p2: TPointF;
begin
p1.X := 37.623612;
p1.Y := 55.756339;
p2.X := 37.626016;
p2.Y := 55.758928;
ShowMessage(FloatToStr(GetGeoDistance(p1, p2)));
end;
|
