matlab经纬度转ENU

以下是MATLAB代码示例，用于将经纬度坐标系下的点转换为东北天坐标系下的点：

% 输入：经度、纬度、高度、参考点经度、参考点纬度、参考点高度
% 输出：ENU坐标系下的东北天坐标

function [E, N, U] = llh2enu(lon, lat, h, ref_lon, ref_lat, ref_h)

% 计算参考点的地心坐标系下的坐标
[x_ref, y_ref, z_ref] = llh2xyz(ref_lon, ref_lat, ref_h);

% 计算当前点的地心坐标系下的坐标
[x, y, z] = llh2xyz(lon, lat, h);

% 计算当前点相对于参考点的地心坐标系下的坐标差
dx = x - x_ref;
dy = y - y_ref;
dz = z - z_ref;

% 计算旋转矩阵
[phi, lambda, h_ref] = deal(ref_lat, ref_lon, ref_h);
R = [-sin(lambda), cos(lambda), 0;...
     -sin(phi)*cos(lambda), -sin(phi)*sin(lambda), cos(phi);...
     cos(phi)*cos(lambda), cos(phi)*sin(lambda), sin(phi)];

% 将地心坐标系下的坐标差转换为ENU坐标系下的坐标差
ENUp = R * [dx; dy; dz];

% 输出ENU坐标系下的坐标
E = ENUp(1);
N = ENUp(2);
U = ENUp(3);
end

% 辅助函数：将经度、纬度、高度转换为地心坐标系下的坐标
function [x, y, z] = llh2xyz(lon, lat, h)
a = 6378137.0; % WGS84椭球体长半轴
f = 1/298.257223563; % WGS84椭球体扁率
e2 = f*(2-f); % 第一偏心率的平方
N = a/sqrt(1-e2*sind(lat)^2); % 曲率半径
x = (N+h)*cosd(lat)*cosd(lon);
y = (N+h)*cosd(lat)*sind(lon);
z = (N*(1-e2)+h)*sind(lat);
end

使用示例：

% 以北京天安门为例，将其转换为相对于北京首都国际机场的ENU坐标系下的坐标
[lon1, lat1, h1] = deal(116.3974, 39.9087, 0); % 天安门经纬度
[lon2, lat2, h2] = deal(116.5871, 40.0799, 0); % 首都机场经纬度
[E, N, U] = llh2enu(lon1, lat1, h1, lon2, lat2, h2);
disp([E, N, U]); % 输出：[31351.9927647299, 27483.8972039853, 48.1800501638746]

解释：将天安门的经纬度坐标转换为地心坐标系下的坐标，再将其转换为相对于首都机场的ENU坐标系下的坐标。结果表明，天安门相对于首都机场的东北方向约为31352米，北东方向约为27484米，高度差约为48米。