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

大地坐标与空间直角坐标的换算

一、实验目的

编写大地坐标与空间直角坐标相互转换的程序，并对格式化文件数据进

行计算，验证程序。

二、实验内容：

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

转换公式：

x(Nh)cosBcosL

y(Nh)cosBsinL（1）

2

z[N(1e)h]sinB

其中：L为经度，B为纬度，h为大地高，Na为卯酉圈曲率半径，

1e2sin2B

a2b2

e为第一偏心率，a为旋转椭球长半轴，b为短半轴。

a

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、z

x=(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;

ifdms0

azimuth=-azimuth;

else

azimuth=azimuth;

end

return

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

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

Larctan(y/x)

zNe2sinB

Barctan()

22

xy

22

xy

hN

cosB

注意计算纬度时需要用到迭代（条件：abs(B0-B)10.^-6），可用

az

Barctan()作为初始值。

bx2y2

创建以xyz2geo命名的函数，函数输入输出格式为

[L,B,h]=xyz2geo(x,y,z)

xyz2geo函数如下：