|
这篇文章主要实现XYZ转BLH及XYZ转NEU两个功能，是为了简化一下有需要的同志的工作，原理方面不过多赘述，直接将代码附上
function [B,L,H] = xyz2blh(X,Y,Z)
%function:
% 将坐标由空间直角坐标系转换为地心坐标系
% input:
% X、Y、Z-->空间直角坐标系坐标
% output:
% B、L、H-->地心坐标系坐标
%所用常数(CGCS2000坐标系)
a=6378137;
% e2=0.00669438002290; %第一偏心率的平方
e2=0.0066943799013; %第一偏心率的平方(WGS84)
%1.计算大地经度L
if (X==0) && (Y>0)
L=90;
elseif X==0 && Y1e-10
t_B0=t_B1;
t_B1=(a*e2*t_B0/sqrt(1+t_B0^2-e2*t_B0^2)+Z)/sqrt(X^2+Y^2);
end
B=atan2(t_B1,1);
%3.计算大地高
N=a/sqrt(1-e2*(sin(B))^2);
H=sqrt(X^2+Y^2)/cos(B)-N;
考虑到由空间直角坐标系与北天东坐标系间的转换也比较常用，特补充记录一下
function p1_NEU = xyz2neu(p0_XYZ,p1_XYZ)
%function:
% 将p1点坐标由空间直角坐标系转换为站心坐标系(站心坐标系原点为p)
% input:
% p0_XYZ-->测站点坐标
% p1_XYZ-->空间点坐标
% output:
% N、E、U-->地心坐标系坐标
% 1.构造旋转矩阵R
x0=p0_XYZ(1); y0=p0_XYZ(2); z0=p0_XYZ(3);
x1=p1_XYZ(1); y1=p1_XYZ(2); z1=p1_XYZ(3);
[B,L,~]=xyz2blh(x0,y0,z0); %调用函数得到纬度B、经度L
R=[-sin(B)*cos(L) -sin(B)*sin(L) cos(B)
-sin(L) cos(L) 0
cos(B)*cos(L) cos(B)*sin(L) sin(B)];
% 2.实现转换
p1_NEU=R*[x1-x0;y1-y0;z1-z0];
end
|