北京1954坐标转为经纬度坐标

种种原因最近（20200716）碰到了这个北京54坐标如何转换为经纬度坐标，网上找了不少内容，实在没有一个是比较贴心的，看到了使用Arcgis可以转换，但是测试之后还是果断放弃，真tm的费劲，果断放弃。

折腾两天之后发现一个非常sweet的软件，下面贴出来分享给各位，直接上菜吧：

Part--1：软件名字：CoordTrans.exe

此软件简单易用，看了基本差不多，如果想进行其他的坐标转换就在相应的模块中选择尝试，中央子午线的设置要参考自己坐标转换点的一个比较接近的经度，软件到此位置。

在使用的时候我碰到了一个问题就是，此软件在批量读取坐标点的时候遇到了问题，没有解决。

Part-2: 批量转换Matlab程序 

function BL_to_xy()
%北京54转经纬度坐标
format long
a=6378245;%长半轴
f = 0.00335232986925914;% 扁率
L0 = 81;%中央子午线经度
data = xlsread('修改为自己的文件位置');%读取需要转换的原始数据
[row,col] = size(data);
new_data = zeros(row,5);
new_data(:,3:5) = data(:,1:3); 
for i = 1:row
    Lat = data(i,2);
    Lon = data(i,1);
    PRO = GKZS(Lat,Lon,L0,a,f);
    new_data(i,1:2) = PRO;
end
xlswrite('修改为自己的文件位置',new_data);%转换后的结果保存
end

%% Sub function 此处的子程序参考网上的分享，没有标注出来源，
%各位看到此子程序的来源可以给提供一下，谢谢！！
function PRO = GKZS(Lat,Lon,L0,a,f)
% Lat: Latitude(rad) 
% Lon: longitude(rad)
% REF//程鹏飞,成英燕,文汉江,等.2000国家大地坐标系实用宝典[M].
%    //北京:测绘出版社,2008:144-148.
%L0是中央子午线，根据不同区块自己设定
Lat = Lat*pi/180;
Lon = Lon*pi/180;
MedLon = L0*pi/180; %中央子午线经度
Eth.R0 = a;%长半轴
Eth.f = f; % 扁率
Eth.e12 = 2*Eth.f - Eth.f*Eth.f; % 第一偏心率额e^2
Eth.e22 = Eth.e12/((1 - Eth.f)*(1 - Eth.f));% 第二偏心率额e'^2
%% 高斯投影正算公式
RN = Eth.R0/sqrt(1 - Eth.e12*sin(Lat)*sin(Lat));
Lon = Lon - MedLon;
Lon2 = Lon*Lon;
Lon4 = Lon2*Lon2;
tnLat = tan(Lat);
tn2Lat = tnLat*tnLat;
tn4Lat = tn2Lat*tn2Lat;
csLat = cos(Lat);
cs2Lat = csLat*csLat;
cs4Lat = cs2Lat*cs2Lat;
Eta2 = Eth.e22*cs2Lat;
NTBLP = RN*tnLat*cs2Lat*Lon2;
coe1 = (5 - tn2Lat + 9*Eta2 + 4*Eta2*Eta2)*cs2Lat*Lon2/24;
coe2 = (61 - 58*tn2Lat + tn4Lat)*cs4Lat*Lon4/720;
x = Merdian(Eth,Lat) + NTBLP*(0.5 + coe1 + coe2);
NBLP = RN*csLat*Lon;
coe3 = (1 - tn2Lat + Eta2)*cs2Lat*Lon2/6;
coe4 = (5 - 18*tn2Lat + tn4Lat + 14*Eta2 - 58*tn2Lat*Eta2)*cs4Lat*Lon4/120;
y = NBLP*(1 + coe3 + coe4) + 500000;%此处有修改，后续在笔记中加以体现
%
PRO=[x y];
end
%% sub function
function X0 = Merdian(Eth,Lat)
% REF//过家春.子午线弧长公式的简化及其泰勒级数解释[J].测绘学报,2014,43(2):125-130.
S0 = Eth.R0*(1 - Eth.e12);
e2 = Eth.e12;
e4 = e2*e2;
e6 = e4*e2;
e8 = e6*e2;
e10 = e8*e2;
e12 = e10*e2;
A1 = 1 + 3*e2/4 + 45*e4/64 + 175*e6/256 + 11025*e8/16384 + 43659*e10/65536 + 693693*e12/1048576;
B1 = 3*e2/8 + 15*e4/32 + 525*e6/1024 + 2205*e8/4096 + 72765*e10/131072 + 297297*e12/524288;
C1 = 15*e4/256 + 105*e6/1024 + 2205*e8/16384 + 10395*e10/65536 + 1486485*e12/8388608;
D1 = 35*e6/3072 + 105*e8/4096 + 10395*e10/262144 + 55055*e12/1048576;
E1 = 315*e8/131072 + 3465*e10/524288 + 99099*e12/8388608;
F1 = 693*e10/1310720 + 9009*e12/5242880;
G1 = 1001*e12/8388608;
X0 = S0*(A1*Lat - B1*sin(2*Lat) + C1*sin(4*Lat) - D1*sin(6*Lat) +...
    E1*sin(8*Lat) - F1*sin(10*Lat) + G1*sin(12*Lat));
end
%% Writed by themingyi 20200715

如有任何问题欢迎评论互动！