Cosine Formula
** cos A = cos B cos C + sin B sin C cos A
Sine Formula
Sin a Sin b Sin c
** ─── = ─── = ───
Sin A Sin B Sin C
球面間兩點的距離
** cos Z = sin ψa sinψb + cosψa cosψb cosλ
┌ ψa 、 ψb ── 兩點分別的緯度
└ λ ── 兩點間的經度差
Z
兩點間的距離S = ── × 地球實際周長
2 π
平面坐标系和地球经纬坐标系的变换
unit Coordinate;
{-------------------------------------------------------------------------------
经纬和平面坐标系的变换
{-------------------------------------------------------------------------------
{ Date:
{ Last update:
{ Coder:zjb
{ Platform:Win2K,Delphi 6
{------------------------------------------------------------------------------}
{
{
{------------------------------------------------------------------------------}
interface
uses
Math;
type
TGeographyConst = record
case Integer of
0: (dtShort, dtEasting, dtscale, dtLong, dtNorthing, dtCenter: Double);
1: (EarthParams: array[0..5] of Double);
end;
TDot = record
case Integer of
0: (Longitude, Latitude: Double);
1: (X, Y: Double);
end;
{ 经纬转换平面坐标系中的坐标 }
function JW2XY(Dot: TDot): TDot;
{ 平面坐标变换成经纬度 }
function XY2JW(Dot: TDot): TDot;
{ }
function LLToXY(Dot: TDot): TDot;
function XYToLL(Dot: TDot): TDot;
{ 初始化参数，在使用本单元前，必须进行初始化 }
function SetConst(Geogconst: TGeographyConst): Boolean;
implementation
{const
dsLongAxis = 6378245;
dsShortAxis = 6356863;
dsFalseEasting = 500000;
dsFalseNorthing = 0;
dsScale = 1;
dsCMRad= 114*Pi/180 ;
// dsCentralMeridian,
}
const
CDRadDeg: Double = 180 / Pi;
CDDegRad: Double = Pi / 180;
var
ArEarthConst: array[0..9] of Double;
dsLongAxis, dsShortAxis, dsFalseEasting, dsFalseNorthing,
dsScale, dsCMRad, dsDev: Double;
function ComputeR2R(t: integer): Double;
begin
if t <= 0 then Result := 1.0
else Result := (t shl 2 - 2) / t * ComputeR2R(t - 1);
end;
procedure ComputeEarthConst(LongRadius: double);
function CalcStepData(iStep: integer; rdBase: Double): Double;
var
i: integer;
rdSum, rdDelta, rdTmp: Double;
begin
rdSum := 0.0; i := iStep; rdDelta := 1;
while (rdDelta / (rdDelta + rdSum) > 1.0E-10) do
begin
rdTmp := Sqr(ComputeR2R(i));
rdDelta := (i shl 1 + 1) * rdTmp * Power(rdBase, i);
rdSum := rdSum + rdDelta; Inc(i);
end;
result := rdSum;
end;
var
dAe, dSum, dBase: double;
i: integer;
begin
dAe := LongRadius * (1.0 - dsDev);
dBase := dsDev / 16.0;
dSum := CalcStepData(0, dBase);
ArEarthConst[0] := dSum * dAe;
for i := 0 to 8 do
begin
dSum := CalcStepData(i + 1, dBase);
ArEarthConst[i + 1] := dSum * dAe * Power(4, i) / (i shl 1 + 1) / ComputeR2R(i);
end;
end;
function JW2XY(Dot: TDot): TDot;
begin
result := LLToXY(Dot);
end;
function LLToXY(Dot: TDot): TDot;
// 计算中心子午带的长度
function CalculateX0(rLat: Double): Double;
var
rdX0, rdSLat: Double;
i: integer;
begin
rdX0 := 0.0; rdSLat := Sin(rLat);
for i := 1 to 6 do
rdX0 := rdX0 - Power(rdSLat, i shl 1 - 1) * ArEarthConst[i];
result := rdX0 * Cos(rLat) + ArEarthConst[0] * rLat;
end;
var
rdX0, rdJD, rdWD, rdt, rdt2, rdm, rdN, rdDeta,
rdConst11, rdConst12, rdConst21, rdConst22: Double;
begin
if Dot.Latitude > Abs(89.999998) then
rdJD := dsCMRad * CDDegRad
else
rdJD := Dot.Longitude * CDDegRad;
rdWD := Dot.Latitude * CDDegRad;
rdX0 := CalculateX0(rdWD);
rdt := Tan(rdWD);
rdt2 := rdt * rdt;
rdm := (rdJD - dsCMRad) * Cos(rdWD);
rdN := dsLongAxis / Sqrt(1.0 - dsDev * Power(sin(rdWD), 2));
rdDeta := dsDev / (1 - dsDev) * Power(cos(rdWD), 2);
rdConst11 := (5.0 - rdt2 + (9.0 + 4.0 * rdDeta) * rdDeta) * Power(rdm, 4) / 24.0;
rdConst12 := (61.0 - (58.0 + rdt2) * rdt2) * Power(rdm, 6) / 720.0;
rdConst21 := (1 - rdt2 + rdDeta) * Power(rdm, 3) / 6.0;
rdConst22 := (5.0 - (18.0 + rdt2) * rdt2 + (14.0 - 58.0 * rdt2) * rdDeta) * Power(rdm, 5) / 120.0;
Result.Longitude := (rdN * (rdm + rdConst21 + rdConst22)) * dsScale + dsFalseEasting;
Result.Latitude := (rdX0 + rdN * rdt * (rdm * rdm * 0.5 + rdConst11 + rdConst12)) * dsScale + dsFalseNorthing;
end;
function SetConst(Geogconst: TGeographyConst): Boolean;
begin
try
with Geogconst do
begin
dsLongAxis := dtLong;
dsShortAxis := dtShort;
dsFalseEasting := dtEasting;
dsFalseNorthing := dtNorthing;
dsScale := dtscale;
dsCMRad := dtCenter * Pi / 180;
end;
dsDev := 1 - Sqr(1 - (dsLongAxis - dsShortAxis) / dsLongAxis);
ComputeEarthConst(dsLongAxis);
result := true;
except
result := false;
end;
end;
function XY2JW(Dot: TDot): TDot;
begin
result := XYToLL(Dot);
end;
function XYToLL(Dot: TDot): TDot;
// 计算低点的纬度
function CalculateLatitude0(rX: Double): Double;
var
rdValue, rdSValue, rdDelta, rdLat0: Double;
i: integer;
begin
rdValue := rX / ArEarthConst[0];
repeat
rdDelta := 0.0;
rdSValue := sin(rdValue);
for i := 1 to 5 do
rdDelta := rdDelta - Power(rdSValue, i shl 1 - 1) * ArEarthConst[i];
rdLat0 := (rX - rdDelta * Cos(rdValue)) / ArEarthConst[0];
rdDelta := Abs((rdLat0 - rdValue) / rdLat0);
rdValue := rdLat0;
until (rdDelta < 1.0E-10);
result := rdLat0;
end;
var
rdX, rdY: Double;
rdLat0, rdN, rdt, rdt2, rdV2, rdDeta, rdYN, rdYN2, rdDevL,
rdConst11, rdConst12, rdConst21, rdConst22: Extended;
begin
rdX := (Dot.Latitude - dsFalseNorthing) / dsScale;
rdY := (Dot.Longitude - dsFalseEasting) / dsScale;
rdLat0 := CalculateLatitude0(rdX);
rdt := tan(rdLat0);
rdt2 := rdt * rdt;
rdN := dsLongAxis / sqrt(1.0 - dsDev * Power(sin(rdLat0), 2));
rdDeta := dsDev / (1 - dsDev) * Power(cos(rdLat0), 2);
rdV2 := 1.0 + rdDeta;
rdYN := rdY / rdN;
rdYN2 := rdYN * rdYN;
rdConst11 := (5.0 + rdDeta + (3.0 - 9.0 * rddeta) * rdt2) / 12.0;
rdConst12 := (61.0 + (90.0 + 45.0 * rdt2) * rdt2) / 360.0;
rdConst21 := (1.0 + 2.0 * rdt2 + rdDeta) / 6.0;
rdConst22 := (5.0 + (28.0 + 24.0 * rdt2) * rdt2 + (6.0 + 8.0 * rdt2) * rdDeta) / 120.0;
rdDevL := rdYN * (1 - (rdConst21 + rdConst22 * rdYN2) * rdYN2) / cos(rdLat0);
Result.Longitude := CDRadDeg * (rdDevL + dsCMRad);
Result.Latitude := CDRadDeg * (rdLat0 - 0.5 * rdV2 * rdt * rdYN2 * (1 - (rdConst11 + rdConst12 * rdYN2) * rdYN2));
end;
end.