MATLAB实验大地坐标与空间直角坐标的换算程序设计(经典)

大地坐标与空间直角坐标的换算一、实验目的编写大地坐标与空间直角坐标相互转换的程序，并对格式化文件数据进行计算，验证程序。二、实验内容：1、大地坐标向空间直角坐标换算转换公式：BheNzLBhNyLBhNxsin])1([sincos)(coscos)(2（1）其中：L为经度，B为纬度，h为大地高，BeaN22sin1为卯酉圈曲率半径，abae22为第一偏心率，a为旋转椭球长半轴，b为短半轴。WGS84椭球参数：长半轴a=6378137扁率f=1/298.257223563根据上式创建以geo2xyz命名的函数，函数输入输出格式为[x,y,z]=geo2xyz(L,B,h)geo2xyz函数如下：function[x,y,z]=geo2xyz(L,B,h)a=6378137;%椭球长半轴f=1/298.257223563;%椭球扁率b=a*(1-f);%求椭球短半轴e=sqrt(a^2-b^2)/a;%椭球第一偏心率N=a./sqrt(1-(e^2)*(sin(B)).^2);%卯酉圈曲率半径%大地坐标换算为空间直角坐标x、y、zx=(N+h).*cos(B).*cos(L);y=(N+h).*cos(B).*sin(L);z=[N.*(1-e^2)+h].*sin(B);end度分秒转化为弧度函数如下：functionazimuth=dms2rad(dms)%度分秒转弧度函数dms_abs=abs(dms);dgree=fix(dms_abs);min_tem=(dms_abs-dgree)*100;min=fix(min_tem);second=(min_tem-min)*100;azimuth=(dgree+min/60+second/3600)*pi/180;ifdms0azimuth=-azimuth;elseazimuth=azimuth;endreturn2、空间直角坐标向大地坐标换算根据式（1）推导大地坐标向空间直角坐标转换公式：NByxhyxBNezBxyLcos)sinarctan()/arctan(22222注意计算纬度时需要用到迭代（条件：abs(B0-B)10.^-6），可用)arctan(22yxbazB作为初始值。创建以xyz2geo命名的函数，函数输入输出格式为[L,B,h]=xyz2geo(x,y,z)xyz2geo函数如下：function[L,B,h]=xyz2geo(x,y,z)a=6378137;%长半轴f=1/298.257223563;%第一扁率b=a*(1-f);%短半轴e=sqrt(a^2-b^2)/a;%第一偏心率%经度L计算L=atan(y./x).*180./pi;%求纬度，单位是度fori=1:4%for循环函数ifL(i)0L(i)=L(i)%满足条件，输出得到一个L矩阵elseL(i)=180+L(i)endend%最后输出得到循环后的L矩阵%纬度B计算B0=atan((a.*z)./(b.*sqrt(x.^2+y.^2)));%B初始值while1%迭代函数N=a./(sqrt(1-(e^2).*(sin(B0).^2)));%卯酉圈曲率半径NB=atan((z+N.*(e^2).*sin(B0))./(sqrt(x.^2+y.^2)));%求纬度B，单位为度ifabs(B0-B)10^-6break;endabs(B0-B)B0=B;endh=(sqrt(x.^2+y.^2))./cos(B)-N;%计算大地高hB=degree3dms(B.*180/pi);%度－度分秒(dd.mmss)end度转化为度分秒函数：functionB=degree3dms(du)%度－度分秒(dd.mmss)degree=fix(du);min=fix((du-degree)*60);second=(((du-degree)*60-min)*60);B=degree+min/100+second/10000;end3、实例计算验证首先将文件data1.txt中大地坐标转换为空间直角坐标，并将转换后的数据按照格式存贮在文件data2.txt中，data1.txt格式为：经度（dd.mmss）纬度（dd.mmss）大地高data2.txt格式为：xyztest程序如下：clear;clc;[filename,pathname]=uigetfile('*.*');%文件查找窗口file=fullfile(pathname,filename);%合并路径文件名A=importdata(file);dms=A.data(:,1:2);%提取A.data中所有行和第1第2列azimuth=dms2rad(dms);%度分秒转化成弧度[x,y,z]=geo2xyz(azimuth(:,1),azimuth(:,2),A.data(:,3));%输出得到xyz矩阵fid=fopen('data2.txt','wt');%写入文件路径fprintf(fid,'xyz\n');[m,n]=size(B);%输入矩阵m,n取行列号fori=1:mforj=1:nifj==nfprintf(fid,'%f\n',B(i,j));elsefprintf(fid,'%f\t',B(i,j));endendendfclose(fid);得到的data2和data1比较如下：然后将文件data2.txt中空间直角坐标转换为大地坐标，并将计算结果按照格式存贮在文件data3.txt中，data3.txt格式为：经度（dd.mmss）纬度（dd.mmss）大地高test程序如下：clear;clc;[filename,pathname]=uigetfile('*.*');%文件查找窗口file=fullfile(pathname,filename);%合并路径文件名A=importdata(file);[L,B,h]=xyz2geo(A.data(:,1),A.data(:,2),A.data(:,3));C=[degree3dms(L),B,h];%输出得到B矩阵fid=fopen('data3.txt','wt');%写入文件路径fprintf(fid,'经度（dd.mmss）纬度（dd.mmss）大地高\n');[m,n]=size(C);%输入矩阵m,n取行列号fori=1:mforj=1:nifj==nfprintf(fid,'%10.7f\n',C(i,j));elsefprintf(fid,'%10.7f\t',C(i,j));endendendfclose(fid);得到的data3和data1比较如下：