直角坐标与大地坐标转换计算matlab小程序

实验目的

1. 掌握直角坐标系及大地坐标系坐标系统的概念和基本知识。
2. 了解直角坐标系及大地坐标系下坐标的转换方法。

实验内容及步骤

1. 运行主程序以获取实时导航数据。

2. 运行本实验程序，在“GPS 时刻”框中记录所有实时 GPS 时间

   

3. 选择一个 GPS 时刻，查看该时刻接收机位置在 ECEF 和 WGS84 坐标系下对应的坐标，并按照附表一对数据进行记录。

实验报告

1. 按附表格式整理实验数据。
2. 编写一个小程序，输入直角坐标系坐标，自动计算大地坐标系坐标；输入大地坐标系坐标，自动计算直角坐标系坐标。

附表一

直角坐标系坐标	大地坐标系坐标
X: -2191120.12	纬度: 22.24900022
Y: 5171136.65	经度: 114.02696478
Z: 3000863.42	高度: 43.08

直角坐标系坐标	大地坐标系坐标
X: -2193324.43	纬度: 26.24900189
Y: 5172237.06	经度: 113502700356
Z: 3000224.74	高度: 44.53

直角坐标系坐标	大地坐标系坐标
X: -2192223.08	纬度: 24.24899735
Y: 5175532.98	经度: 118.02700715
Z: 3002261.86	高度: 40.39

实现三组直角坐标与大地坐标转换计算的流程图

实现过程描述

1. 开始：程序开始执行。
2. 输入直角坐标 (x, y, z)：提示用户输入直角坐标，例如 [x, y, z] 格式。
3. 计算大地坐标系 (经度, 纬度, 高度)：将输入的直角坐标转换为大地坐标系的经度、纬度和高度。
4. 显示大地坐标系：将计算得到的大地坐标系的经度、纬度和高度显示给用户。
5. 输入大地坐标 (纬度, 经度, 高度)：提示用户输入大地坐标，例如 [纬度, 经度, 高度] 格式。
6. 计算直角坐标 (x, y, z)：将输入的大地坐标转换为直角坐标系的 x、y 和 z 坐标。
7. 显示直角坐标：将计算得到的直角坐标系的 x、y 和 z 坐标显示给用户。
8. 重复输入输出：判断是否已经完成三次输入输出，若未完成，则返回到输入直角坐标步骤；若已完成，则结束程序。
9. 结束：程序执行结束。

实现三组直角坐标与大地坐标转换计算的程序运行结果

代码

function coordinate_conversion()
 
    % 循环三次，每次交替输入直角坐标和大地坐标
    for i = 1:3
        % 输入直角坐标，计算大地坐标
        disp(['请输入第 ', num2str(i), ' 组直角坐标系坐标 (x, y, z):']);
        input_str = input('格式应为[x, y, z]: ', 's');
        coords = extractNumbers(input_str);
        
        if numel(coords) < 3
            disp('输入坐标数目不足，需要三个数字(x, y, z).');
            continue;
        end
        
        x = coords(1);
        y = coords(2);
        z = coords(3);
        
        % 计算大地坐标系坐标 (纬度, 经度, 高度)
        [latitude, longitude, altitude] = cart2geo(x, y, z);
        
        % 显示计算得到的大地坐标系坐标
        disp(['第 ', num2str(i), ' 组直角坐标系坐标转换为大地坐标系坐标 (纬度, 经度, 高度):']);
        disp(['纬度: ', num2str(latitude, '%.6f')]);
        disp(['经度: ', num2str(longitude, '%.6f')]);
        disp(['高度: ', num2str(altitude, '%.6f')]);
        
        % 输入大地坐标，计算直角坐标
        disp(['请输入第 ', num2str(i), ' 组大地坐标系坐标 (纬度, 经度, 高度):']);
        geo_str = input('格式应为[纬度, 经度, 高度]: ', 's');
        geo_coords = extractNumbers(geo_str);
        
        if numel(geo_coords) < 3
            disp('输入坐标数目不足，需要三个数字(纬度, 经度, 高度).');
            continue;
        end
        
        latitude = geo_coords(1);
        longitude = geo_coords(2);
        altitude = geo_coords(3);


        
        % 计算直角坐标系坐标 (x, y, z)
        [x_calc, y_calc, z_calc] = geo2cart(longitude, latitude, altitude);
        
        % 显示计算得到的直角坐标系坐标
        disp(['第 ', num2str(i), ' 组大地坐标系坐标转换为直角坐标系坐标 (x, y, z):']);
        disp(['x: ', num2str(x_calc, '%.2f')]);
        disp(['y: ', num2str(y_calc, '%.2f')]);
        disp(['z: ', num2str(z_calc, '%.2f')]);
    end
 
end
 
function nums = extractNumbers(str)
    % 使用正则表达式提取字符串中的数字
    expr = '[-+]?\d*\.?\d+';
    nums = str2double(regexp(str, expr, 'match'));
end
 
function [latitude, longitude, altitude] = cart2geo(x, y, z)
    % 实现直角坐标系到大地坐标系的转换
    a = 6378137.0;  % WGS84椭球长半轴 (单位: 米)
    f = 1 / 298.257223563;  % WGS84椭球扁率
    e2 = f * (2 - f);  % 第一偏心率的平方
    
    % 计算经度 (longitude) 和纬度 (latitude)
    longitude = atan2(y, x) * 180 / pi;
    p = sqrt(x^2 + y^2);
    latitude = atan2(z, p * (1 - e2)) * 180 / pi;
    
    % 计算大地高 (altitude)
    N = a / sqrt(1 - e2 * sind(latitude)^2);
    altitude = p / cosd(latitude) - N;
end
 
function [x, y, z] = geo2cart(longitude, latitude, altitude)
    % 实现大地坐标系到直角坐标系的转换
    a = 6378137.0;  % WGS84椭球长半轴 (单位: 米)
    f = 1 / 298.257223563;  % WGS84椭球扁率
    e2 = f * (2 - f);  % 第一偏心率的平方
    
    % 将经度和纬度转换为弧度
    longitude_rad = longitude * pi / 180;
    latitude_rad = latitude * pi / 180;
    
    % 计算子午线曲率半径
    N = a / sqrt(1 - e2 * sind(latitude)^2);
    
    % 计算直角坐标系坐标 (x, y, z)
    x = (N + altitude) * cosd(latitude) * cosd(longitude);
    y = (N + altitude) * cosd(latitude) * sind(longitude);
    z = (N * (1 - e2) + altitude) * sind(latitude);
end

---