如何将经纬坐标转换为大地坐标 加急 ( 积分: 10 )

如何将经纬坐标转换为大地坐标 加急

如何将经纬坐标转换为大地坐标 加急

TGeoCoor = record
B, L, H: Double;
end;

TXYZCoor = record
X, Y, Z: Double;
end;


const
a84 = 6378137.0;
b84 = 6356752.31;
e1sqr = (a84 * a84 - b84 * b84) / (a84 * a84);
e2sqr = (a84 * a84 - b84 * b84) / (b84 * b84);

procedure WGS84TOXYZ(const GeoCoor: TGeoCoor
var XYZCoor: TXYZCoor);
var
N: Double;
begin
N := a84 / sqrt(1 - e1sqr * sqr(sin(GeoCoor.B)));

XYZCoor.X := (N + GeoCoor.H) * cos(GeoCoor.B) * cos(GeoCoor.L);
XYZCoor.Y := (N + GeoCoor.H) * cos(GeoCoor.B) * sin(GeoCoor.L);
XYZCoor.Z := (N + GeoCoor.H - N * e1sqr) * sin(GeoCoor.B);
end;

procedure XYZToWGS84(const XYZCoor: TXYZCoor
var GeoCoor: TGeoCoor);
var
p, T, sT, cT, N: Double;
begin
p := sqrt(XYZCoor.X * XYZCoor.X + XYZCoor.Y * XYZCoor.Y);
T := Get_atan(p * b84, XYZCoor.Z * a84);
sT := sin(T)
cT := cos(T);
GeoCoor.B := arctan((XYZCoor.Z + e2sqr * b84 * sT * sT * sT)
/ (p - e1sqr * a84 * cT * cT * cT));
GeoCoor.L := Get_atan(XYZCoor.X, XYZCoor.Y);
N := a84 / sqrt(1.0 - e1sqr * sin(GeoCoor.B) * sin(GeoCoor.B));
GeoCoor.H := p / cos(GeoCoor.B) - N;
end;