# coding=utf-8
import math
"""2000国家大地坐标系(cgcs2000)椭球常数"""
a = 6378137.0 # 椭球长半轴
b = 6356752.31414036 # 椭球短半轴
f = 1.0 / 298.257222101 # 椭球扁率
e1 = 0.0066943800229 # 椭球第一偏心率 e^2
ee = 0.00673949677548 # 椭球第二偏心率 e`^2
c = 6399593.62586 # 极点子午圈曲率半径
def xytoBl():
m0 = a * (1 - e1)
m2 = 1.5 * e1 * m0
m4 = 5.0 / 4.0 * e1 * m2
m6 = 7.0 / 6.0 * e1 * m4
m8 = 9.0 / 8.0 * e1 * m6
a0 = m0 + 0.5 * m2 + 3.0 / 8.0 * m4 + 5.0 / 16.0 * m6 + 35.0 / 128.0 * m8
a2 = 0.5 * m2 + 0.5 * m4 + 15.0 / 32.0 * m6 + 7.0 / 16.0 * m8
a4 = 1.0 / 8.0 * m4 + 3.0 / 16.0 * m6 + 7.0 / 32.0 * m8
a6 = 1.0 / 32.0 * m6 + 1.0 / 16.0 * m8
a8 = 1.0 / 128.0 * m8
X = float(input("请输入x坐标值:"))
Y = float(input("请输入y坐标值:")) - 36000000 - 500000
l0 = float(input("请输入中央子午线:")) / 180.0 * math.pi
# 迭代求取低点纬度
B = [0] * 20
f = [0] * 20
B[0] = float(X) / a0
f[0] = -(a2 / 2.0 * math.sin(2 * B[0]) + a4 / 4.0 * math.sin(4.0 * B[0]) - a6 / 6.0 * math.sin(
6.0 * B[0]) + a8 / 8.0 * math.sin(8.0 * B[0]))
deltad = B[1] - B[0]
for i in range(10):
f[i] = -(a2 / 2.0 * math.sin(2 * B[i]) + a4 / 4.0 * math.sin(4.0 * B[i]) - a6 / 6.0 * math.sin(
6.0 * B[i]) + a8 / 8.0 * math.sin(8.0 * B[i]))
B[i + 1] = (X - f[i]) / a0
deltad = B[i + 1] - B[i]
if deltad < pow(10, -10):
i += 1
break
"""参数"""
tf = math.tan(B[i])
mf = c / pow(math.sqrt(1 + ee * pow(math.cos(B[i]), 2)), 3)
nf = c / math.sqrt(1 + ee * pow(math.cos(B[i]), 2))
n = math.sqrt(ee) * math.cos(B[i])
# longitude纬度计算公式
latitude = B[i] - (tf / (2.0 * mf) * Y * (Y / nf)) * (1.0 - (1.0 / 12.0) * (
5.0 + 3.0 * pow(tf, 2) + pow(n, 2) - 9.0 * pow(n, 2) * pow(tf, 2)) *
pow(Y / nf, 2) + 1.0 / 360.0 * (
61.0 + 90.0 * pow(tf, 2) + 45.0 * pow(tf, 4)) *
pow(Y / nf, 4))
# longitude 纬度计算公式
longitude = 1.0 / math.cos(B[i]) * (Y / nf) * (
1.0 - 1.0 / 6.0 * (1.0 + 2.0 * pow(tf, 2) + pow(n, 2)) * pow(Y / nf, 2) + 1.0 / 120.0 * (
5.0 + 28.0 * pow(tf, 2) + 24.0 * pow(tf, 4) + 6.0 * pow(n, 2) + 8.0 * pow(n, 2) * pow(tf, 2)) * pow(Y / nf,
4))
# 将经纬度的弧度转换为度
lat = latitude * 180 / 3.14159265358979
long = (longitude + l0) * 180 / 3.14159265358979
return lat, long
2000坐标系高斯反算;
计算过程中的角度均为弧度;
最后将弧度转换为角度, 角度 = 弧度 * 180 / pi