基于matlab的大地坐标与直角坐标间的转换_基于matlab的坐标转换

基于matlab的大地坐标与直角坐标间的转换由刀豆文库小编整理，希望给你工作、学习、生活带来方便，猜你可能喜欢“基于matlab的坐标转换”。

测量程序设计 实验报告

换算

实验名称：大地坐标与空间直角坐标的 实验四 大地坐标与空间直角坐标的换算

一、实验目的编写大地坐标与空间直角坐标相互转换的程序，并对格式化文件数据进行计算，验证程序。

二、实验内容：

1、大地坐标向空间直角坐标换算 转换公式：

x(Nh)cosBcosLy(Nh)cosBsinL（1）z[N(1e2)h]sinB其中：L为经度，B为纬度，h为大地高，Na1esinB22为卯酉圈曲率半径，ea2b2为第一偏心率，a为旋转椭球长半轴，b为短半轴。aWGS84椭球参数：长半轴 a=6378137 扁率 f = 1/298.257223563 根据上式创建以geo2xyz命名的函数，函数输入输出格式为 [x, y, z] = geo2xyz(L, B, h)

2、空间直角坐标向大地坐标换算

根据式（1）推导大地坐标向空间直角坐标转换公式：

Larctan(y/x)zNe2sinBBarctan()

22xyhx2y2NcosBaz注意计算纬度时需要用到迭代，可用Barctan()作为初始值。

22bxy创建以xyz2geo命名的函数，函数输入输出格式为 [L, B, h] = xyz2geo(x, y, z)

三、实验步骤

1、大地坐标向空间直角坐标换算 主程序：

%%大地坐标向空间直角坐标换算

%函数的输入输出格式为[x,y,z]=geo2xyz(L,B,h)[filename,pathname] = uigetfile('*.txt','请选择打开的数据文件');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');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(a^2-b^2)/a;N=a./(sqrt(1-e^2.*(sin(B)).^2));%N为卯酉圈半径率，e为第一偏心率 x=(N+h).*cos(B).*cos(L);y=(N+h).*cos(B).*sin(L);z=(N*(1-e^2)+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中的空间直角坐标系转换为大地坐标，并将计算结果按照格式存储在文件data3.txt中

%data3.txt格式为：经度（ddmm）纬度（ddmm）大地高 [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 = [pathname_out, filename_out];fid = fopen(fileout,'wt');fprintf(fid,' 经度(ddmm）纬度(ddmm）大地高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(a^2-b^2))/a;%第一偏心率 %% 经度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-(e^2).*(sin(B0).^2)));%卯酉圈曲率半径N B=atan2(z+N.*(e^2).*sin(B0),(sqrt(x.^2+y.^2)));if abs(B0-B)

H=(sqrt(x.^2+y.^2))./cos(B)-N;B=degree2dms(B.*180/pi);%纬度B End

function dms= degree2dms(jiaodu)%度－>度分秒（dd.mm）degree = fix(jiaodu);mimute = fix((jiaodu-degree)*60);second =((jiaodu-degree)*60-mimute)*60;dms = degree+mimute/100+second/10000;