球面距离
Author: MrBaseball34
Calcuate distance between two global points
Note: Click Title to view in Edit Box for easier copying.
I found a VB unit to calculate two points on
a globe and thought I'd translate it to Delphi
for everyone. I also found another function
on Experts-Exchange, that was posted by
gary_williams, called GreatCircleDistance.
Here they are:
uses ...,Math;
// *****************************************************************
// this function get the arccos function from arctan function
// *****************************************************************
function ACos(ARad: Extended): Extended;
begin
Result := 0.0;
if Abs(ARad) <> 1 then
Result := PI/2 - ArcTan(ARad / Sqrt(1 - ARad * ARad))
else
if ARad = -1 then
Result := PI;
end;
// *****************************************************************
// this function converts decimal degrees to radians
// *****************************************************************
function Deg2Rad(ADeg: Extended): Extended;
begin
Result := (ADeg * PI) / 180;
end;
// *****************************************************************
// this function converts radians to decimal degrees
// *****************************************************************
function Rad2Deg(ARad: Extended): Extended;
begin
Result := (ARad * 180) / PI;
end;
// *****************************************************************
// this function returns the distance between two points
// ALat1, ALat2: Points of Lattitude
// ALon1, ALon2: Points of Longitude
// AUnit: Unit of measure to convert distance to
// m = Statue miles
// n = nautical miles
// k = kilomoters
// *****************************************************************
function Distance(ALat1, ALon1, ALat2, ALon2: Extended; AUnit: Char):Extended;
var
LTheta : Extended;
LDist : Extended;
begin
LTheta := ALon1 - ALon2;
LDist := Sin(Deg2Rad(ALat1)) * Sin(Deg2Rad(ALat2)) + Cos(Deg2Rad(ALat1)) *
Cos(Deg2Rad(ALat2)) * Cos(Deg2Rad(LTheta));
LDist := ACos(LDist);
LDist := Rad2Deg(LDist);
Result := LDist * 60 * 1.1515;
case AUnit of
'k':
Result := Result * 1.609344;
'n':
Result := Result * 0.8684;
end; {case}
end;
{
This function returns the approximate distance between two points on the
Earth's surface. Specify Lat1, Long1, Lat2, and Long2 in degrees.
By convention, south latitudes and west longitudes should be negative.
Radius should be set to 3963.1676 if you want miles. Note that this function
assumes a perfectly spherical earth.
Author: Gary Williams
}
function GreatCircleDistance(Lat1Deg,
Lon1Deg,
Lat2Deg,
Lon2Deg,
Radius: Double): Double;
var
Lat1Rad: Double;
Lon1Rad: Double;
Lat2Rad: Double;
Lon2Rad: Double;
begin
Lat1Rad := DegToRad(Lat1Deg);
Lon1Rad := DegToRad(Lon1Deg);
Lat2Rad := DegToRad(Lat2Deg);
Lon2Rad := DegToRad(Lon2Deg);
Result := Radius * ArcCos(Sin(Lat1Rad) *
Sin(Lat2Rad) +
Cos(Lat1Rad) *
Cos(Lat2Rad) *
Cos(Lon2Rad - Lon1Rad));
end;
procedure TForm1.Button1Click(Sender: TObject);
var
d: Double;
begin
// Austin, TX = 30.5144, -97.6555
// Round Rock, TX = 30.4257, -97.7142
d := GreatCircleDistance(30.5144, -97.6555, 30.4257, -97.7142, 3963.1676);
Label1.Caption := FloatToStr(RoundTo(d, -4));
Label2.Caption := FloatToStr(RoundTo(Distance(30.5144, -97.6555, 30.4257, -97.7142, 'm'), -4)) + ' Statue miles' + #13 +
FloatToStr(RoundTo(Distance(30.5144, -97.6555, 30.4257, -97.7142, 'k'), -4)) + ' Kilomoeters' + #13 +
FloatToStr(RoundTo(Distance(30.5144, -97.6555, 30.4257, -97.7142, 'n'), -4)) + ' Nautical miles';
end;
---------------------------------------
// -----------------------------------------------------------------------------
// GeoDegToRad
// PURPOSE : convert 2 Lat/Lon pairs from DEG to RAD
// -----------------------------------------------------------------------------
procedure GeoDegToRad(Latitude1, Longitude1, Latitude2, Longitude2: double;
var Lat1Rad, Lon1Rad, Lat2Rad, Lon2Rad: double);
begin
Lat1Rad := DegToRad(Latitude1);
Lon1Rad := DegToRad(Longitude1);
Lat2Rad := DegToRad(Latitude2);
Lon2Rad := DegToRad(Longitude2);
end;
// -----------------------------------------------------------------------------
// DistanceRad
// PURPOSE : calculate the distance between two points (given as geographic coordinates in Latitude and Longitude)
// -----------------------------------------------------------------------------
function DistanceRad(Latitude1, Longitude1, Latitude2, Longitude2: double): double;
var
d, dLa1, dLa2, dLo1, dLo2: double;
begin
if (Latitude1 = Latitude2) and (Longitude1 = Longitude2) then
Result := NOVALIDRESULT
else
begin
GeoDegToRad(Latitude1, Longitude1, Latitude2, Longitude2, dLa1, dLo1, dLa2, dLo2);
d := (Sin(dLa1) * Sin(dLa2) + (Cos(dLa1) * Cos(dLa2) * Cos(dLo2 - dLo1)));
Result := arccos(d);
end;
end;
// -----------------------------------------------------------------------------
// GetBearingBetweenCoordsRad
// PURPOSE : Get the bearing between two Lat/Lon's
// -----------------------------------------------------------------------------
function GetBearingBetweenCoordsRad(Latitude1, Longitude1, Latitude2,
Longitude2: double): double;
var
d: double;
dLa1, dLa2, dLo1, dLo2: double;
begin
try
d := DistanceRad(Latitude1, Longitude1, Latitude2, Longitude2);
if d = 0 then
begin
Result := 0;
Exit;
end;
GeoDegToRad(Latitude1, Longitude1, Latitude2, Longitude2, dLa1, dLo1, dLa2, dLo2);
if Sin(dlo2 - dlo1) < 0 then
begin
Result := 2 * pi - ArcCos((sin(dla2) - sin(dla1) * cos(d)) / (sin(d) * cos(dla1)));
end
else
begin
Result := ArcCos((sin(dla2) - sin(dla1) * cos(d)) / (sin(d) * cos(dla1)));
end;
except
Result := 0;
end;
end;
function GetBearingBetweenCoordsDeg(Latitude1, Longitude1, Latitude2,
Longitude2: double): double;
begin
Result := RadToDeg(GetBearingBetweenCoordsRad(Latitude1, Longitude1,
Latitude2, Longitude2));
end;
---------------------------------------
type Twinkel = record
grad : integer;
min : integer;
sek : integer;
end;
{==============================================================================}
{ Umrechnung von Geographischen Koordinaten in Gauss-Krueger-Koordinaten }
{ Formel: Grossmann,W., GeodAbbildungen, 1964, Seite 151 }
{ Parameter: geo.Breite (Grad.Min.Sek) in Altgrad : Twinkel }
{ geo.Laenge (Grad.Min.Sek) in Altgrad : Twinkel }
{ Zielsystemnummer (Meridiankennziffer) : longint }
{ Rechtswert (X) im Zielsystem : double }
{ Hochwert (Y) im Zielsystem : double }
{==============================================================================}
procedure GeoGk(br,la:Twinkel;sy:Longint;var x,y:double);
const
{26}
rho = 180 / pi;
var
brDezimal,laDezimal,rm,e2,c,bf,g,co,g2,g1,t,dl,fa,grad,min,sek :extended;
begin
{25}
e2 := 0.0067192188;
{27}
c := 6398786.849;
{in Dezimal}
{Breite}
brDezimal := br.grad + br.min / 60 + br.sek / 3600;
{Laenge}
laDezimal := la.grad + la.min / 60 + la.sek / 3600;
{64}
bf := brDezimal / rho;
{65}
g := 111120.61962 * brDezimal
-15988.63853 * sin(2*bf)
+16.72995 * sin(4*bf)
-0.02178 * sin(6*bf)
+0.00003 * sin(8*bf);
{70}
co := cos(bf);
{71}
g2 := e2 * (co * co);
{72}
g1 := c / sqrt(1+g2);
{73}
t := sin(bf) / cos(bf); {=tan(t)}
{74}
dl := laDezimal - sy * 3;
{77}
fa := co * dl / rho;
{78}
y := g
+ fa * fa * t * g1 / 2
+ fa * fa * fa * fa * t * g1 * (5 - t * t + 9 * g2) / 24;
{81}
rm := fa * g1
+ fa * fa * fa * g1 * (1 - t * t + g2) / 6
+ fa * fa * fa * fa * fa * g1 * (5 - 18 * t * t * t * t * t * t) / 120;
{84}
x := rm + sy * 1000000 + 500000;
end;
{==============================================================================}
{ Umrechnung von Geographischen Koordinaten in Gauss-Krueger-Koordinaten }
{ Formel: Grossmann,W., GeodAbbildungen, 1964, Seite 153 }
{ Parameter: Rechtswert : double }
{ Hochwert : double }
{ geo.Breite (Grad.Min.Sek) in Altgrad : Twinkel }
{ geo.Laenge (Grad.Min.Sek) in Altgrad : Twinkel }
{==============================================================================}
procedure GkGeo(rw,hw:double;var br,la:Twinkel);
const
{26}
rho= 180 / pi;
var
rm,e2,c,bI,bII,bf,co,g2,g1,t,fa,dl,gb,gl,grad,min,sek :extended;
mKen :integer;
begin
{25}
e2 := 0.0067192188;
{27}
c := 6398786.849;
{32}
mKen := trunc(rw / 1000000);
{33}
rm := rw - mKen * 1000000 - 500000;
{34}
bI := hw / 10000855.7646;
{35}
bII := bI * bI;
{36}
bf := 325632.08677 *bI *((((((0.00000562025
* bII - 0.00004363980)
* bII + 0.00022976983)
* bII - 0.00113566119)
* bII + 0.00424914906)
* bII - 0.00831729565)
* bII + 1);
{43}
bf := bf / 3600 / rho;
{44}
co := cos(bf);
{45}
g2 := e2 * (co * co);
{46}
g1 := c / sqrt(1 + g2);
{47}
t := sin(bf) / cos(bf); {=tan(bf)}
{51}
fa := rm / g1;
{52}
gb := bf
- fa * fa * t * (1 + g2) / 2
+ fa * fa * fa * fa * t * (5 + 3 * t * t + 6 * g2 - 6 * g2 * t * t) / 24;
{55}
gb := gb * rho;
{56}
dl := fa
- fa * fa * fa * (1 + 2 * t * t + g2) / 6
+ fa * fa * fa * fa * fa * (1 + 28 * t * t + 24 * t * t * t * t) / 120;
{59}
gl := dl *rho / co + mKen * 3;
{in grad.min.sek}
{Breite}
grad:=int(gb);
sek:=60*(gb-grad);
min:=int(sek);
sek:=60*(sek-min);
br.grad := trunc(grad);
br.min := trunc(min);
br.sek := trunc(sek);
{Laenge}
grad:=int(gl);
sek:=60*(gl-grad);
min:=int(sek);
sek:=60*(sek-min);
la.grad := trunc(grad);
la.min := trunc(min);
la.sek := trunc(sek);
end;