1、测量程序设计 实验报告 实验名称：大地坐标与空间直角坐标的 换算 实验四 大地坐标与空间直角坐标的换算 一、实验目的 编写大地坐标与空间直角坐标相互转换的程序，并对格式化文件数据 进行计算，验证程序。 二、实验内容： 1、大地坐标向空间直角坐标换算 转换公式： （1）BheNzLyxsin)1(co2 其中：L 为经度，B 为纬度，h 为大地高， 为卯酉圈曲率半径，ea2si1 为第一偏心率，a 为旋转椭球长半轴，b 为短半轴。be2 WGS84 椭球参数：长半轴 a = 6378137 扁率 f = 1/298.257223563 根据上式创建以 geo2xyz 命名的函数，函数输入输出格式
2、为 x, y, z = geo2xyz (L, B, h) 2、空间直角坐标向大地坐标换算 根据式（1）推导大地坐标向空间直角坐标转换公式： NByxhez xyLcos)sinart()/c22 注意计算纬度时需要用到迭代，可用 作为初始值。)arctn(2yxbz 创建以 xyz2geo 命名的函数，函数输入输出格式为 L, B, h = xyz2geo (x, y, z) 三、实验步骤 1、 大地坐标向空间直角坐标换算 主程序： %大地坐标向空间直角坐标换算 %函数的输入输出格式为x,y,z=geo2xyz(L,B,h) filename,pathname = uigetfile(*.t
3、xt,请选择打开的数据文件); file = pathname, 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);
4、 fprintf(fid,%15.7f %15.7f %15.7fn,A); close(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);
5、 z=(N*(1-e2)+h).*sin(B); end function rad=dms2rad(jiaodu) %度分秒弧度（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; end 2、 空间直角坐标向大地坐标换算 主程序： % 将文件 data.2.txt 中的空间直角坐标系转换为大地坐标，并将计算结果按 照格式存储
6、在文件 data3.txt 中 %data3.txt 格式为：经度（ddmmss） 纬度（ddmmss） 大地高 filename,pathname=uigetfile(*.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 = pathn
7、ame_out, filename_out; fid = fopen(fileout,wt); fprintf(fid, 经度(ddmmss） 纬度(ddmmss） 大地高n); 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)/a a=6378137; %长半轴 b=a*(1-f); %短半轴 e=(sqrt(a2-b2)/a; %第一偏心率 %
8、经度 L 的计算 L=atan2(y,x); L=rad2dms(L); % 纬度 B 的计算 B0=atan2(a*z),(b.*sqrt(x.2+y.2); % B 初始值 while 1 N=a./(sqrt(1-(e2).*(sin(B0).2); %卯酉圈曲率半径 N B=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;
