大地坐标系转换工具 - 北京54西安80转WGS84四参数七参数计算器
大地坐标系转换工具
北京54/西安80 → WGS84/CGCS2000，支持四参数和七参数转换
使用说明
支持四参数和七参数进行坐标系转换，需要输入实测参数。
完整转换（投影坐标→WGS84）
计算原理
完整转换流程
从投影坐标（如北京54）转换为WGS84经纬度需要以下步骤：
- 高斯投影反算：投影坐标(X,Y) → 大地坐标(B,L)
- 大地转空间：大地坐标 → 空间直角坐标(X,Y,Z)
- 七参数转换：源椭球空间坐标 → 目标椭球空间坐标
- 空间转大地：空间直角坐标 → WGS84大地坐标
四参数（平面）
X₂ = ΔX + (1+K×10⁻⁶)(X₁cosθ - Y₁sinθ) Y₂ = ΔY + (1+K×10⁻⁶)(X₁sinθ + Y₁cosθ)
七参数（布尔莎模型）
X₂ = ΔX + (1+K×10⁻⁶)(X₁ - Rz×Y₁ + Ry×Z₁) Y₂ = ΔY + (1+K×10⁻⁶)(Rz×X₁ + Y₁ - Rx×Z₁) Z₂ = ΔZ + (1+K×10⁻⁶)(-Ry×X₁ + Rx×Y₁ + Z₁)
参数说明
- ΔX,ΔY,ΔZ：三个平移参数（米）
- Rx,Ry,Rz：三个旋转参数（秒）
- K：尺度因子（ppm）
代码实现
JavaScript
// 完整坐标转换：投影坐标 → WGS84经纬度
const ellipsoids = {
BJ54: { a: 6378245, f: 1 / 298.3 },
XA80: { a: 6378140, f: 1 / 298.257 },
WGS84: { a: 6378137, f: 1 / 298.257223563 }
};
// 高斯投影反算
function gaussInverse(x, y, L0, ellipsoid) {
const { a, f } = ellipsoid;
const e2 = 2 * f - f * f;
const Bf = x / (a * (1 - e2 / 4));
const Wf = Math.sqrt(1 - e2 * Math.sin(Bf) ** 2);
const Nf = a / Wf;
const Mf = a * (1 - e2) / Wf ** 3;
const tf = Math.tan(Bf);
const B = Bf - (Nf * tf / Mf) * (y ** 2 / (2 * Nf ** 2));
const L = L0 + y / (Nf * Math.cos(Bf));
return { B: B * 180 / Math.PI, L: L * 180 / Math.PI };
}
// 大地坐标→空间直角坐标
function blhToXyz(B, L, H, ellipsoid) {
const { a, f } = ellipsoid;
const e2 = 2 * f - f * f;
const N = a / Math.sqrt(1 - e2 * Math.sin(B * Math.PI / 180) ** 2);
const x = (N + H) * Math.cos(B * Math.PI / 180) * Math.cos(L * Math.PI / 180);
const y = (N + H) * Math.cos(B * Math.PI / 180) * Math.sin(L * Math.PI / 180);
const z = (N * (1 - e2) + H) * Math.sin(B * Math.PI / 180);
return { x, y, z };
}
// 七参数转换
function sevenParam(x, y, z, dx, dy, dz, rx, ry, rz, k) {
const rxr = rx * Math.PI / 648000;
const ryr = ry * Math.PI / 648000;
const rzr = rz * Math.PI / 648000;
const m = 1 + k * 1e-6;
return {
x: dx + m * (x - rzr * y + ryr * z),
y: dy + m * (rzr * x + y - rxr * z),
z: dz + m * (-ryr * x + rxr * y + z)
};
}
// 空间直角坐标→大地坐标
function xyzToBlh(x, y, z, ellipsoid) {
const { a, f } = ellipsoid;
const e2 = 2 * f - f * f;
const L = Math.atan2(y, x);
let B = Math.atan2(z, Math.sqrt(x ** 2 + y ** 2));
for (let i = 0; i < 5; i++) {
const N = a / Math.sqrt(1 - e2 * Math.sin(B) ** 2);
const H = Math.sqrt(x ** 2 + y ** 2) / Math.cos(B) - N;
B = Math.atan2(z, Math.sqrt(x ** 2 + y ** 2) * (1 - e2 * N / (N + H)));
}
return { B: B * 180 / Math.PI, L: L * 180 / Math.PI };
}
// 使用示例：北京54投影坐标转WGS84
const projX = 4426641.2, projY = -2148744.5, L0 = 117 * Math.PI / 180;
const bl = gaussInverse(projX, projY, L0, ellipsoids.BJ54);
const xyz1 = blhToXyz(bl.B, bl.L, 0, ellipsoids.BJ54);
const xyz2 = sevenParam(xyz1.x, xyz1.y, xyz1.z, -29.7, -130.6, -95.8, 0.0025, 0.0033, 0, 2.4);
const wgs84 = xyzToBlh(xyz2.x, xyz2.y, xyz2.z, ellipsoids.WGS84);
console.log('WGS84:', wgs84);Java
public class DatumTransform {
static class Ellipsoid {
double a, f;
Ellipsoid(double a, double f) { this.a = a; this.f = f; }
}
static Ellipsoid BJ54 = new Ellipsoid(6378245, 1.0/298.3);
static Ellipsoid WGS84 = new Ellipsoid(6378137, 1.0/298.257223563);
// 高斯投影反算
static double[] gaussInverse(double x, double y, double L0, Ellipsoid e) {
double e2 = 2 * e.f - e.f * e.f;
double Bf = x / (e.a * (1 - e2 / 4));
double Wf = Math.sqrt(1 - e2 * Math.sin(Bf) * Math.sin(Bf));
double Nf = e.a / Wf;
double B = Bf - (Nf * Math.tan(Bf) / (e.a * (1 - e2) / Math.pow(Wf, 3))) *
(y * y / (2 * Nf * Nf));
double L = L0 + y / (Nf * Math.cos(Bf));
return new double[] {Math.toDegrees(B), Math.toDegrees(L)};
}
// 大地坐标→空间直角坐标
static double[] blhToXyz(double B, double L, double H, Ellipsoid e) {
double e2 = 2 * e.f - e.f * e.f;
double Br = Math.toRadians(B), Lr = Math.toRadians(L);
double N = e.a / Math.sqrt(1 - e2 * Math.sin(Br) * Math.sin(Br));
return new double[] {
(N + H) * Math.cos(Br) * Math.cos(Lr),
(N + H) * Math.cos(Br) * Math.sin(Lr),
(N * (1 - e2) + H) * Math.sin(Br)
};
}
// 七参数转换
static double[] sevenParam(double x, double y, double z,
double dx, double dy, double dz, double rx, double ry, double rz, double k) {
double rxr = rx * Math.PI / 648000;
double ryr = ry * Math.PI / 648000;
double rzr = rz * Math.PI / 648000;
double m = 1 + k * 1e-6;
return new double[] {
dx + m * (x - rzr * y + ryr * z),
dy + m * (rzr * x + y - rxr * z),
dz + m * (-ryr * x + rxr * y + z)
};
}
}Python
import math
class Ellipsoid:
def __init__(self, a, f):
self.a = a
self.f = f
BJ54 = Ellipsoid(6378245, 1/298.3)
WGS84 = Ellipsoid(6378137, 1/298.257223563)
def gauss_inverse(x, y, L0, ellipsoid):
"""高斯投影反算"""
e2 = 2 * ellipsoid.f - ellipsoid.f ** 2
Bf = x / (ellipsoid.a * (1 - e2 / 4))
Wf = math.sqrt(1 - e2 * math.sin(Bf) ** 2)
Nf = ellipsoid.a / Wf
Mf = ellipsoid.a * (1 - e2) / Wf ** 3
tf = math.tan(Bf)
B = Bf - (Nf * tf / Mf) * (y ** 2 / (2 * Nf ** 2))
L = L0 + y / (Nf * math.cos(Bf))
return math.degrees(B), math.degrees(L)
def blh_to_xyz(B, L, H, ellipsoid):
"""大地坐标→空间直角坐标"""
e2 = 2 * ellipsoid.f - ellipsoid.f ** 2
N = ellipsoid.a / math.sqrt(1 - e2 * math.sin(math.radians(B)) ** 2)
x = (N + H) * math.cos(math.radians(B)) * math.cos(math.radians(L))
y = (N + H) * math.cos(math.radians(B)) * math.sin(math.radians(L))
z = (N * (1 - e2) + H) * math.sin(math.radians(B))
return x, y, z
def seven_param(x, y, z, dx, dy, dz, rx, ry, rz, k):
"""七参数转换"""
rxr = rx * math.pi / 648000
ryr = ry * math.pi / 648000
rzr = rz * math.pi / 648000
m = 1 + k * 1e-6
return (
dx + m * (x - rzr * y + ryr * z),
dy + m * (rzr * x + y - rxr * z),
dz + m * (-ryr * x + rxr * y + z)
)
def xyz_to_blh(x, y, z, ellipsoid):
"""空间直角坐标→大地坐标"""
e2 = 2 * ellipsoid.f - ellipsoid.f ** 2
L = math.atan2(y, x)
B = math.atan2(z, math.sqrt(x ** 2 + y ** 2))
for _ in range(5):
N = ellipsoid.a / math.sqrt(1 - e2 * math.sin(B) ** 2)
H = math.sqrt(x ** 2 + y ** 2) / math.cos(B) - N
B = math.atan2(z, math.sqrt(x ** 2 + y ** 2) * (1 - e2 * N / (N + H)))
return math.degrees(B), math.degrees(L)
# 使用示例
B, L = gauss_inverse(4426641.2, -2148744.5, math.radians(117), BJ54)
xyz1 = blh_to_xyz(B, L, 0, BJ54)
xyz2 = seven_param(*xyz1, -29.7, -130.6, -95.8, 0.0025, 0.0033, 0, 2.4)
wgs84 = xyz_to_blh(*xyz2, WGS84)
print(f"WGS84: B={wgs84[0]:.8f}°, L={wgs84[1]:.8f}°")