基于matlab的大地坐标与直角坐标间的转换精编版

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

二、实验内容：
1、大地坐标向空间直角坐标换算
转换公式：
B
h e N z L B h N y L
B h N x sin ])1([sin cos )(cos cos )(2+-=+=+= （1）
其中：L 为经度，B 为纬度，h 为大地高，B e a
N 22sin 1-=为卯酉圈曲率半径，
a
b a e 2
2-=为第一偏心率，a 为旋转椭球长半轴，b 为短半轴。

WGS84椭球参数：长半轴 a = 6378137
扁率 f = 1/298.257223563
根据上式创建以geo 2xyz 命名的函数，函数输入输出格式为
[x, y, z] = geo 2xyz (L, B, h)
2、空间直角坐标向大地坐标换算
根据式（1）推导大地坐标向空间直角坐标转换公式： N B
y x h y
x B Ne z B x y L -+=++==cos )sin arctan()
/arctan(2
2222 注意计算纬度时需要用到迭代，可用)arctan(22y
x b az B +=作为初始值。

创建以xyz2geo 命名的函数，函数输入输出格式为
[L, B, h] = xyz 2geo (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 z\n'); fprintf(fid,'%15.7f %15.7f %15.7f\n',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格式为：经度（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 = [pathname_out, filename_out];
fid = fopen(fileout,'wt');
fprintf(fid,' 经度(ddmmss）纬度(ddmmss）大地高\n'); fprintf(fid,'%10.7f %10.7f %7.3f\n',[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)<10^-6
break;
end
abs(B0-B)
B0=B;
end
%% 大地高H的计算
H=(sqrt(x.^2+y.^2))./cos(B)-N;
B=degree2dms(B.*180/pi);%纬度B
End
function dms= degree2dms(jiaodu)
%度－>度分秒（dd.mmss）
degree = fix(jiaodu);
mimute = fix((jiaodu-degree)*60); second = ((jiaodu-degree)*60-mimute)*60; dms = degree+mimute/100+second/10000;。