MATLAB程序实现经纬度转换成平面坐标

 近期搜了下经纬度坐标转换成直角坐标的程序和原理啥的，哎，真难啊，其实道理都懂吧，下面直接贴代码，

主要是加了个for循环而言，优化输出效果。

%%%MATLAB程序实现经纬度转换成平面尔坐标：
M_PI=3.14159265358979323846;
L = 6381372 * M_PI * 2; %地球周长  
W = L; % 平面展开后，x轴等于周长  
H = L / 2; % y轴约等于周长一半  
mill = 2.3; % 米勒投影中的一个常数，范围大约在正负2.3之间  
position=[120.7015202,36.37423; 
120.6987175,36.37457569;
120.6997954,36.37591239;
120.70691,36.37579616;
120.7056165,36.37248342;
120.7031731,36.37753964;
120.6928965,36.37800457;
120.6943337,36.37521499;
120.6973521,36.37876006;
120.6962022,36.37643544;
120.7011609,36.37905063;
120.6939026,36.38021291;
120.6983582,36.38056159;
120.7025263,36.38120084;
120.6914592,36.38201441;
120.6960585,36.38247931;
120.7005141,36.38276987;
120.6998673,36.37079794;
120.6928965,36.37079794;
120.6964897,36.36824059;
120.6969209,36.37143727;
120.7052571,36.36899618;
120.7088504,36.37021674;
120.7087066,36.36731063;
120.7130185,36.36829872;
120.6896626,36.36661314;
120.6937588,36.36242812;
120.6993643,36.38741865;
120.7129466,36.37201847;
120.7002266,36.36428816;]%30个经纬度坐标点
n=size(position,1);

%%lon=120.7015202;%经度
%%lat=36.37423;%纬度
new_position=[];
for i =1:n
    lon=position(i,1);
    lat=position(i,2);
    x = lon * M_PI / 180; % 将经度从度数转换为弧度  
    y = lat * M_PI / 180; %将纬度从度数转换为弧度  
    y1 = 1.25 * log(tan(0.25 * M_PI + 0.4 * y)); % 米勒投影的转换  
    % 弧度转为实际距离  
    dikaerX = (W / 2) + (W / (2 * M_PI)) * x ; %笛卡尔坐标x
    dikaerY = (H / 2) - (H / (2 * mill)) * y1 ;%笛卡尔坐标y
    new_position(i,1)=dikaerX;
    new_position(i,2)=dikaerY;
    fprintf('第%d个点的',i)
    fprintf('坐标是=(%f %f);',new_position(i,1),new_position(i,2))
    fprintf('\n')
end

参考文章链接：<https://blog.csdn.net/qq_32693445/article/details/79597213?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-3.nonecase&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-3.nonecase>