大地坐标与空间直角坐标互相转换

引言

坐标转换代码，由大地坐标（经、纬度、高）与空间直角坐标（XYZ）互相转换。MATLAB代码，下面代码默认转换的是WGS84坐标系下的空间直角和大地坐标，可以根据需要选择椭球参数。

Matlab代码

clc; clear; close all

% X Y Z B L H
% -2036837.469 5210777.883 3052481.957 28.77899495 111.349998 56.69428242
% 表2，测试数据出自文献：韩晓爽, 潘超, 刘宇哲,等. 星载激光测高仪脚点坐标解算方法及机载验证[J]. 遥测遥控, 2020, 41(3):10.

pos = [-2036837.469, 5210777.883, 3052481.957];
[lat,lon,h]=xyz2ell(pos);
LAT = rad2deg(lat);
LON = rad2deg(lon);
[LAT LON h]

LAT = 28.77899495;
LON = 111.349998;
h = 56.69428242;
[X,Y,Z] = BLh2XYZ(LAT,LON,h);
[X,Y,Z]

BLH2XYZ

function [X,Y,Z] = BLh2XYZ(LAT,LON,h) 

% Example:
%(BLh)WGS84-(XYZ)WGS84
% LAT = 40.9987167395335;
% LON = 39.7652393428761;
% h = 51.403;

refell = 1;

switch refell
    case 1
        % IERS 2003 numerical standards
        % ellipsoid parameters for xyz2ellip.m
        a_tidefree = 6378136.6; %m      Equatorial radius of the Earth
        f_tidefree = 1/298.25642;     % Flattening factor of the Earth
        a = a_tidefree;  %m      Equatorial radius of the Earth
        f = f_tidefree;       % Flattening factor of the Earth
    case 2
        % GRS 80 (http://www.bkg.bund.de/nn_164850/geodIS/EVRS/EN/References/...
        %  Definitions/Def__GRS80-pdf,templateId=raw,property=publication...
        %  File.pdf/Def_GRS80-pdf.pdf)
        a_grs80    = 6378137;
        f_grs80    = 0.00335281068118;
        a = a_grs80;   %m      Equatorial radius of the Earth
        f = f_grs80;        % Flattening factor of the Earth
    case 3
        % WGS84
        a=6378137;
        f=1/298.25722356;
    case 4
        % Hayford
        a=6378388;
        f=1/297;
end

b = a-f*a;

lat=LAT*pi/180;
lon=LON*pi/180;
e2=(a^2-b^2)/a^2;
N=a/sqrt(1-e2*(sin(lat)^2));
X=(N+h)*cos(lat)*cos(lon);
Y=(N+h)*cos(lat)*sin(lon);
Z=(N*(1-e2)+h)*sin(lat);

XYZ2BLH

% ************************************************************************
%   Description:
%   Transformation from Cartesian coordinates X,Y,Z to ellipsoidal 
%   coordinates lam,phi,elh. based on Occam subroutine transf.
%
%   Input:										
%      pos = [x,y,z]                 [m,m,m]
%               can be a matrix: number of rows = number of stations
%
%   Output:
%      coor_ell = [lat,lon,h]      [rad,rad,m]
% 
%   External calls: 	
%      global   a_...              Equatorial radius of the Earth [m]     
%               f_...              Flattening factor of the Earth
%
%   Coded for VieVS: 
%   17 Dec 2008 by Lucia Plank
%
%   Revision:
%   
% *************************************************************************
function [lat,lon,h]=xyz2ell(pos)

% choose reference ellipsoid
%       1 ...... tide free
%       2 ...... GRS80
%       3 ...... WGS84
refell =3;

switch refell
    case 1
        % IERS 2003 numerical standards (tide free crust)
        % ellipsoid parameters for xyz2ellip.m
        a = 6378136.6; %m      Equatorial radius of the Earth
        f = 1/298.25642;     % Flattening factor of the Earth
    case 2
        % GRS 80
        % (http://www.bkg.bund.de/nn_164850/geodIS/EVRS/EN/References/...
        %  Definitions/Def__GRS80-pdf,templateId=raw,property=publication...
        %  File.pdf/Def_GRS80-pdf.pdf)
        a = 6378137;  %m      Equatorial radius of the Earth
        f = 0.00335281068118;       % Flattening factor of the Earth
    case 3
        % WGS84
        a=6378137;
        f=1/298.25722356;
    case 4
        % Hayford
        a=6378388;
        f=1/297;
end

e2=2*f-f^2;

lon=angle(pos(:,1)+1i*pos(:,2));

lat=angle(sqrt(pos(:,1).^2+pos(:,2).^2)+1i*pos(:,3));
for j=1:10
  N=a./sqrt(1-e2*sin(lat).*sin(lat));
  h=sqrt(pos(:,1).^2+pos(:,2).^2)./cos(lat)-N;
  lat=angle(sqrt(pos(:,1).^2+pos(:,2).^2).*((1-e2)*N+h)+1i*pos(:,3).*(N+h));
end

%lat=cart2phigd(pos);