type
// 笛卡尔坐标系
tagCRDCARTESIAN = record
X, Y, z: double;
end;
CRDCARTESIAN = tagCRDCARTESIAN;
// 大地坐标系
tagCRDGEODETIC = record
longitude: double; // 经度
latitude: double; // 纬度
height: double; // 大地高,可设为0
end;
CRDGEODETIC = tagCRDGEODETIC;
function GeodeticToCartesian(pcg: CRDGEODETIC): CRDCARTESIAN;
var
b: double; // 纬度度数
L: double; // 经度度数
L0: double; // 中央经线度数
L1: double; // L-L0
t: double; // tanB
m: double; // ltanB
N: double; // 卯酉圈曲率半径
q2: double;
X: double; // 高斯平面纵坐标
Y: double; // 高斯平面横坐标
s: double; // 赤道至纬度B的经线弧长
f: double; // 参考椭球体扁率
e2: double; // 椭球第一偏心率
a: double; // 参考椭球体长半轴
// b:double; //参考椭球体短半轴
a1: double;
a2: double;
a3: double;
a4: double;
b1: double;
b2: double;
b3: double;
b4: double;
c0: double;
c1: double;
c2: double;
c3: double;
Datum, prjno, zonewide: Integer;
IPI: double;
begin
Datum := 84; // 投影基准面类型：北京54基准面为54，西安80基准面为80，WGS84基准面为84
prjno := 0; // 投影带号
zonewide := 3;
IPI := 0.0174532925199433333333; // 3.1415926535898/180.0
b := pcg.latitude; // 纬度
L := pcg.longitude; // 经度
if (zonewide = 6) then
begin
prjno := trunc(L / zonewide) + 1;
L0 := prjno * zonewide - 3;
end
else
begin
prjno := trunc((L - 1.5) / 3) + 1;
L0 := prjno * 3;
end;
if (Datum = 54) then
begin
a := 6378245;
f := 1 / 298.3;
end
else if (Datum = 84) then
begin
a := 6378137;
f := 1 / 298.257223563;
end;
L0 := L0 * IPI;
L := L * IPI;
b := b * IPI;
e2 := 2 * f - f * f; // (a*a-b*b)/(a*a);
L1 := L - L0;
t := Tan(b);
m := L1 * cos(b);
N := a / sqrt(1 - e2 * sin(b) * sin(b));
q2 := e2 / (1 - e2) * cos(b) * cos(b);
a1 := 1 + 3 / 4 * e2 + 45 / 64 * e2 * e2 + 175 / 256 * e2 * e2 * e2 + 11025 /
16384 * e2 * e2 * e2 * e2 + 43659 / 65536 * e2 * e2 * e2 * e2 * e2;
a2 := 3 / 4 * e2 + 15 / 16 * e2 * e2 + 525 / 512 * e2 * e2 * e2 + 2205 /
2048 * e2 * e2 * e2 * e2 + 72765 / 65536 * e2 * e2 * e2 * e2 * e2;
a3 := 15 / 64 * e2 * e2 + 105 / 256 * e2 * e2 * e2 + 2205 / 4096 * e2 * e2 *
e2 * e2 + 10359 / 16384 * e2 * e2 * e2 * e2 * e2;
a4 := 35 / 512 * e2 * e2 * e2 + 315 / 2048 * e2 * e2 * e2 * e2 + 31185 /
13072 * e2 * e2 * e2 * e2 * e2;
b1 := a1 * a * (1 - e2);
b2 := -1 / 2 * a2 * a * (1 - e2);
b3 := 1 / 4 * a3 * a * (1 - e2);
b4 := -1 / 6 * a4 * a * (1 - e2);
c0 := b1;
c1 := 2 * b2 + 4 * b3 + 6 * b4;
c2 := -(8 * b3 + 32 * b4);
c3 := 32 * b4;
s := c0 * b + cos(b) * (c1 * sin(b) + c2 * sin(b) * sin(b) * sin(b) + c3 * sin
(b) * sin(b) * sin(b) * sin(b) * sin(b));
X := s + 1 / 2 * N * t * m * m + 1 / 24 * (5 - t * t + 9 * q2 + 4 * q2 * q2)
* N * t * m * m * m * m + 1 / 720 * (61 - 58 * t * t + t * t * t * t)
* N * t * m * m * m * m * m * m;
Y := N * m + 1 / 6 * (1 - t * t + q2) * N * m * m * m + 1 / 120 *
(5 - 18 * t * t + t * t * t * t - 14 * q2 - 58 * q2 * t * t)
* N * m * m * m * m * m;
Y := Y + 1000000 * prjno + 500000;
Result.X := X;
Result.Y := Y - 38000000;
Result.z := 0;
end;
python笛卡尔坐标系_将GPS的WGS84坐标（大地坐标系）转换为平面坐标（笛卡尔坐标系）的方法...
