高斯投影转换工具 - 大地2000(CGCS2000)高斯-克吕格投影坐标与WGS84经纬度互转
高斯投影转换工具
大地2000(CGCS2000)高斯-克吕格投影坐标与WGS84经纬度互转，支持三度带和六度带，支持批量转换
使用说明
本工具支持大地2000(CGCS2000)高斯-克吕格投影坐标与WGS84经纬度的相互转换。CGCS2000与WGS84坐标系差异在厘米级，本工具视为等效。投影坐标X为北方向（纵坐标），Y为东方向（横坐标，含500km假东偏移）。
度带设置
反算：高斯投影坐标 → WGS84经纬度
批量转换
每行一个坐标，格式：X,Y（X为纵坐标，Y为横坐标，含带号）
投影说明
高斯-克吕格投影（Gauss-Krüger Projection）是一种横轴等角切椭圆柱投影，由中国测绘工作广泛使用，是 CGCS2000（大地2000）坐标系的标准投影方式。
投影特点
- 中央经线投影后为直线，且长度不变
- 离中央经线越远，变形越大
- 等角投影，角度关系保持不变
- 采用分带投影，限制变形范围
分带方式
- 三度带：从东经1.5°起，每隔3°为一带，中央经线为 3°、6°、9°...135°。带号 n = L₀ / 3
- 六度带：从东经0°起，每隔6°为一带，中央经线为 3°、9°、15°...135°。带号 n = (L₀ + 3) / 6
坐标约定
- X（纵坐标）：北方向，赤道起算的子午线弧长
- Y（横坐标）：东方向，中央经线起算，加500km假东偏移避免负值
- 带号：Y值前可加带号，如 Y = 39432827 表示第39带，实际东偏 432827 - 500000 = -67173米
CGCS2000 椭球参数
- 长半轴 a = 6378137.0 米
- 扁率倒数 1/f = 298.257222101
- 与 WGS84 椭球参数差异极小（厘米级），实际应用中视为等效
代码实现
JavaScript
// CGCS2000 高斯-克吕格投影正反算
const CGCS2000 = { a: 6378137.0, f: 1 / 298.257222101 };
const FALSE_EASTING = 500000;
// 正算：经纬度 → 高斯投影坐标
function gaussForward(B, L, L0, ellipsoid = CGCS2000) {
const { a, f } = ellipsoid;
const e2 = 2 * f - f * f;
const e12 = e2 / (1 - e2);
const Br = B * Math.PI / 180;
const l = (L - L0) * Math.PI / 180;
const sinB = Math.sin(Br), cosB = Math.cos(Br), tanB = Math.tan(Br);
const N = a / Math.sqrt(1 - e2 * sinB * sinB);
const t = tanB, eta2 = e12 * cosB * cosB;
// 子午线弧长
const e4 = e2 * e2, e6 = e4 * e2;
const m0 = 1 - e2/4 - 3*e4/64 - 5*e6/256;
const m2 = 3/8 * (e2 + e4/4 + 15*e6/128);
const m4 = 15/256 * (e4 + 3*e6/4);
const m6 = 35/3072 * e6;
const Xm = a * (m0*Br - m2*Math.sin(2*Br) + m4*Math.sin(4*Br) - m6*Math.sin(6*Br));
const l2=l*l, l3=l2*l, l4=l2*l2, l5=l4*l, l6=l4*l2;
const t2=t*t, t4=t2*t2;
const x = Xm + N*t/2*l2 + N*t/24*(5-t2+9*eta2+4*eta2*eta2)*l4
+ N*t/720*(61-58*t2+t4+270*eta2-330*t2*eta2)*l6;
const y = N*l + N/6*(1-t2+eta2)*l3
+ N/120*(5-18*t2+t4+14*eta2-58*t2*eta2)*l5;
return { x, y: y + FALSE_EASTING };
}
// 反算：高斯投影坐标 → 经纬度
function gaussInverse(x, y, L0, ellipsoid = CGCS2000) {
const { a, f } = ellipsoid;
const e2 = 2 * f - f * f;
const e12 = e2 / (1 - e2);
const e4 = e2 * e2, e6 = e4 * e2;
const m0 = 1 - e2/4 - 3*e4/64 - 5*e6/256;
const mu = x / (a * m0);
const e1 = (1 - Math.sqrt(1 - e2)) / (1 + Math.sqrt(1 - e2));
// 底点纬度 Bf
const Bf = mu
+ (3*e1/2 - 27*e1**3/32) * Math.sin(2*mu)
+ (21*e1**2/16 - 55*e1**4/32) * Math.sin(4*mu)
+ (151*e1**3/96) * Math.sin(6*mu)
+ (1097*e1**4/512) * Math.sin(8*mu);
const sinBf = Math.sin(Bf), cosBf = Math.cos(Bf), tanBf = Math.tan(Bf);
const Wf = Math.sqrt(1 - e2 * sinBf * sinBf);
const Nf = a / Wf;
const Mf = a * (1 - e2) / Wf**3;
const tf = tanBf, etaf2 = e12 * cosBf * cosBf;
const y2=y*y, y3=y2*y, y4=y2*y2, y5=y4*y, y6=y4*y2;
const tf2=tf*tf, tf4=tf2*tf2;
const B = Bf - (Nf*tf/Mf) * (
y2/(2*Nf**2)
- (5+3*tf2+etaf2-9*etaf2*tf2)*y4/(24*Nf**4)
+ (61+90*tf2+45*tf4-256*etaf2-36*etaf2*tf2)*y6/(720*Nf**6)
);
const L = L0*Math.PI/180 + y/(Nf*cosBf)
- (1+2*tf2+etaf2)*y3/(6*Nf**3*cosBf)
+ (5+28*tf2+24*tf4+6*etaf2+8*etaf2*tf2)*y5/(120*Nf**5*cosBf);
return { B: B*180/Math.PI, L: L*180/Math.PI };
}
// 使用示例
// 正算：WGS84经纬度 → CGCS2000高斯投影坐标（3度带，中央经线117°）
const fwd = gaussForward(39.90923, 116.397428, 117);
console.log(`X=${fwd.x.toFixed(4)}, Y=${fwd.y.toFixed(4)}`);
// 反算：CGCS2000高斯投影坐标 → WGS84经纬度
// Y需先去除带号和假东偏移：39432827 → 432827 - 500000 = -67173
const inv = gaussInverse(4418501.3, -67173, 117);
console.log(`经度=${inv.L.toFixed(8)}°, 纬度=${inv.B.toFixed(8)}°`);
Java
public class GaussKrugerProjection {
// CGCS2000 椭球参数
static final double A = 6378137.0;
static final double F = 1.0 / 298.257222101;
static final double FALSE_EASTING = 500000;
// 正算：经纬度 → 高斯投影坐标
public static double[] gaussForward(double B, double L, double L0) {
double e2 = 2 * F - F * F;
double e12 = e2 / (1 - e2);
double Br = Math.toRadians(B);
double l = Math.toRadians(L - L0);
double sinB = Math.sin(Br), cosB = Math.cos(Br), tanB = Math.tan(Br);
double N = A / Math.sqrt(1 - e2 * sinB * sinB);
double t = tanB, eta2 = e12 * cosB * cosB;
// 子午线弧长
double e4 = e2 * e2, e6 = e4 * e2;
double m0 = 1 - e2/4 - 3*e4/64 - 5*e6/256;
double m2 = 3.0/8 * (e2 + e4/4 + 15*e6/128);
double m4 = 15.0/256 * (e4 + 3*e6/4);
double m6 = 35.0/3072 * e6;
double Xm = A * (m0*Br - m2*Math.sin(2*Br) + m4*Math.sin(4*Br) - m6*Math.sin(6*Br));
double l2=l*l, l3=l2*l, l4=l2*l2, l5=l4*l, l6=l4*l2;
double t2=t*t, t4=t2*t2;
double x = Xm + N*t/2*l2 + N*t/24*(5-t2+9*eta2+4*eta2*eta2)*l4
+ N*t/720*(61-58*t2+t4+270*eta2-330*t2*eta2)*l6;
double y = N*l + N/6*(1-t2+eta2)*l3
+ N/120*(5-18*t2+t4+14*eta2-58*t2*eta2)*l5;
return new double[]{x, y + FALSE_EASTING};
}
// 反算：高斯投影坐标 → 经纬度
public static double[] gaussInverse(double x, double y, double L0) {
double e2 = 2 * F - F * F;
double e12 = e2 / (1 - e2);
double e4 = e2 * e2, e6 = e4 * e2;
double m0 = 1 - e2/4 - 3*e4/64 - 5*e6/256;
double mu = x / (A * m0);
double e1 = (1 - Math.sqrt(1 - e2)) / (1 + Math.sqrt(1 - e2));
// 底点纬度 Bf
double Bf = mu
+ (3*e1/2 - 27*Math.pow(e1,3)/32) * Math.sin(2*mu)
+ (21*e1*e1/16 - 55*Math.pow(e1,4)/32) * Math.sin(4*mu)
+ (151*Math.pow(e1,3)/96) * Math.sin(6*mu)
+ (1097*Math.pow(e1,4)/512) * Math.sin(8*mu);
double sinBf = Math.sin(Bf), cosBf = Math.cos(Bf), tanBf = Math.tan(Bf);
double Wf = Math.sqrt(1 - e2 * sinBf * sinBf);
double Nf = A / Wf;
double Mf = A * (1 - e2) / Math.pow(Wf, 3);
double tf = tanBf, etaf2 = e12 * cosBf * cosBf;
double y2=y*y, y3=y2*y, y4=y2*y2, y5=y4*y, y6=y4*y2;
double tf2=tf*tf, tf4=tf2*tf2;
double B = Bf - (Nf*tf/Mf) * (
y2/(2*Nf*Nf)
- (5+3*tf2+etaf2-9*etaf2*tf2)*y4/(24*Math.pow(Nf,4))
+ (61+90*tf2+45*tf4-256*etaf2-36*etaf2*tf2)*y6/(720*Math.pow(Nf,6))
);
double L = Math.toRadians(L0) + y/(Nf*cosBf)
- (1+2*tf2+etaf2)*y3/(6*Math.pow(Nf,3)*cosBf)
+ (5+28*tf2+24*tf4+6*etaf2+8*etaf2*tf2)*y5/(120*Math.pow(Nf,5)*cosBf);
return new double[]{Math.toDegrees(B), Math.toDegrees(L)};
}
public static void main(String[] args) {
// 正算
double[] fwd = gaussForward(39.90923, 116.397428, 117);
System.out.printf("X=%.4f, Y=%.4f%n", fwd[0], fwd[1]);
// 反算（Y已去除假东偏移）
double[] inv = gaussInverse(4418501.3, -67173, 117);
System.out.printf("经度=%.8f°, 纬度=%.8f°%n", inv[1], inv[0]);
}
}Python
import math
# CGCS2000 椭球参数
A = 6378137.0
F = 1 / 298.257222101
FALSE_EASTING = 500000
def gauss_forward(B: float, L: float, L0: float) -> tuple[float, float]:
"""正算：经纬度 → 高斯投影坐标"""
e2 = 2 * F - F ** 2
e12 = e2 / (1 - e2)
Br = math.radians(B)
l = math.radians(L - L0)
sinB, cosB, tanB = math.sin(Br), math.cos(Br), math.tan(Br)
N = A / math.sqrt(1 - e2 * sinB ** 2)
t, eta2 = tanB, e12 * cosB ** 2
# 子午线弧长
e4, e6 = e2 ** 2, e2 ** 3
m0 = 1 - e2/4 - 3*e4/64 - 5*e6/256
m2 = 3/8 * (e2 + e4/4 + 15*e6/128)
m4 = 15/256 * (e4 + 3*e6/4)
m6 = 35/3072 * e6
Xm = A * (m0*Br - m2*math.sin(2*Br) + m4*math.sin(4*Br) - m6*math.sin(6*Br))
l2, l3, l4, l5, l6 = l**2, l**3, l**4, l**5, l**6
t2, t4 = t**2, t**4
x = (Xm + N*t/2*l2 + N*t/24*(5-t2+9*eta2+4*eta2**2)*l4
+ N*t/720*(61-58*t2+t4+270*eta2-330*t2*eta2)*l6)
y = (N*l + N/6*(1-t2+eta2)*l3
+ N/120*(5-18*t2+t4+14*eta2-58*t2*eta2)*l5)
return x, y + FALSE_EASTING
def gauss_inverse(x: float, y: float, L0: float) -> tuple[float, float]:
"""反算：高斯投影坐标 → 经纬度"""
e2 = 2 * F - F ** 2
e12 = e2 / (1 - e2)
e4, e6 = e2 ** 2, e2 ** 3
m0 = 1 - e2/4 - 3*e4/64 - 5*e6/256
mu = x / (A * m0)
e1 = (1 - math.sqrt(1 - e2)) / (1 + math.sqrt(1 - e2))
# 底点纬度 Bf
Bf = (mu
+ (3*e1/2 - 27*e1**3/32) * math.sin(2*mu)
+ (21*e1**2/16 - 55*e1**4/32) * math.sin(4*mu)
+ (151*e1**3/96) * math.sin(6*mu)
+ (1097*e1**4/512) * math.sin(8*mu))
sinBf, cosBf, tanBf = math.sin(Bf), math.cos(Bf), math.tan(Bf)
Wf = math.sqrt(1 - e2 * sinBf ** 2)
Nf = A / Wf
Mf = A * (1 - e2) / Wf ** 3
tf, etaf2 = tanBf, e12 * cosBf ** 2
y2, y3, y4, y5, y6 = y**2, y**3, y**4, y**5, y**6
tf2, tf4 = tf**2, tf**4
B = Bf - (Nf*tf/Mf) * (
y2/(2*Nf**2)
- (5+3*tf2+etaf2-9*etaf2*tf2)*y4/(24*Nf**4)
+ (61+90*tf2+45*tf4-256*etaf2-36*etaf2*tf2)*y6/(720*Nf**6)
)
L = (math.radians(L0) + y/(Nf*cosBf)
- (1+2*tf2+etaf2)*y3/(6*Nf**3*cosBf)
+ (5+28*tf2+24*tf4+6*etaf2+8*etaf2*tf2)*y5/(120*Nf**5*cosBf))
return math.degrees(B), math.degrees(L)
if __name__ == "__main__":
# 正算
x, y = gauss_forward(39.90923, 116.397428, 117)
print(f"X={x:.4f}, Y={y:.4f}")
# 反算（Y已去除假东偏移）
B, L = gauss_inverse(4418501.3, -67173, 117)
print(f"经度={L:.8f}°, 纬度={B:.8f}°")