1、测量程序设计实验报告实验名称：大地坐标与空间直角坐标的换算实验四 大地坐标与空间直角坐标的换算一、实验目的编写大地坐标与空间直角坐标相互转换的程序，并对格式化文件数据进行计算，验证程序。二、实验内容：1、大地坐标向空间直角坐标换算转换公式：（1）BheNzLyxsin)1(co2其中：L 为经度，B 为纬度，h 为大地高， 为卯酉圈曲率半径，ea2si1为第一偏心率，a 为旋转椭球长半轴，b 为短半轴。be2WGS84 椭球参数：长半轴 a = 6378137扁率 f = 1/298.257223563根据上式创建以 geo2xyz 命名的函数，函数输入输出格式为x, y, z = geo2x
2、yz (L, B, h)2、空间直角坐标向大地坐标换算根据式（1）推导大地坐标向空间直角坐标转换公式：NByxhezxyLcos)sinart()/c22注意计算纬度时需要用到迭代，可用 作为初始值。)arctn(2yxbz创建以 xyz2geo 命名的函数，函数输入输出格式为L, B, h = xyz2geo (x, y, z)三、实验步骤1、 大地坐标向空间直角坐标换算主程序：%大地坐标向空间直角坐标换算%函数的输入输出格式为x,y,z=geo2xyz(L,B,h)filename,pathname = uigetfile(*.txt,请选择打开的数据文件);file = pathname
3、, filename;data = importdata(file);L=data.data(:,1);B=data.data(:,2);h=data.data(:,3);x,y,z=geo2xyz(L,B,h);A=x,y,z;A=A;filename_out,pathname_out = uiputfile(*.txt,请选择要输出数据文件);fileout = pathname_out, filename_out;fid = fopen(fileout,wt);fprintf(fid, x y zn);fprintf(fid,%15.7f %15.7f %15.7fn,A);close(
4、all);函数：function x,y,z=geo2xyz(L,B,h)%大地坐标经纬度转换成空间直角坐标B=dms2rad(B);L=dms2rad(L);a=6378137;%a 是长半轴f=1/298.257223563;%f 是扁率b=a-a*f;e=sqrt(a2-b2)/a;N=a./(sqrt(1-e2.*(sin(B).2);%N 为卯酉圈半径率，e 为第一偏心率x=(N+h).*cos(B).*cos(L);y=(N+h).*cos(B).*sin(L);z=(N*(1-e2)+h).*sin(B);endfunction rad=dms2rad(jiaodu)%度分秒弧度
5、（rad）degree = fix(jiaodu);mimute = fix(jiaodu-degree)*100);second = (jiaodu-degree-mimute/100)*10000;degree = degree+mimute/60+second/3600;rad=degree/180*pi;end2、 空间直角坐标向大地坐标换算主程序：% 将文件 data.2.txt 中的空间直角坐标系转换为大地坐标，并将计算结果按照格式存储在文件 data3.txt 中%data3.txt 格式为：经度（ddmmss） 纬度（ddmmss） 大地高filename,pathname=u
6、igetfile(*.txt,请选择打开的数据文件);file=pathname,filename;data=importdata(file);x=data.data(:,1);y=data.data(:,2);z=data.data(:,3);L,B,H=xyz2geo(x,y,z);filename_out,pathname_out = uiputfile(*.txt,请选择要输出数据文件);fileout = pathname_out, filename_out;fid = fopen(fileout,wt);fprintf(fid, 经度(ddmmss） 纬度(ddmmss） 大地高n
7、);fprintf(fid,%10.7f %10.7f %7.3fn,L,B,H);fclose(all);函数：function L,B,H=xyz2geo(x,y,z)%将直角坐标转换为大地坐标% 已知:WGS-84 椭球参数f=1/298.257223563;%扁率 f=(a-b)/aa=6378137; %长半轴b=a*(1-f); %短半轴e=(sqrt(a2-b2)/a; %第一偏心率% 经度 L 的计算L=atan2(y,x);L=rad2dms(L);% 纬度 B 的计算B0=atan2(a*z),(b.*sqrt(x.2+y.2); % B 初始值while 1N=a./(sqrt(1-(e2).*(sin(B0).2); %卯酉圈曲率半径 NB=atan2(z+N.*(e2).*sin(B0),(sqrt(x.2+y.2);if abs(B0-B)度分秒（dd.mmss）degree = fix(jiaodu);mimute = fix(jiaodu-degree)*60);second = (jiaodu-degree)*60-mimute)*60;dms = degree+mimute/100+second/10000;
