Matlab中空间直角坐标系中三维速度Vxyz转Ven平面速度

这篇博客介绍了如何使用Matlab进行空间直角坐标XYZ到大地坐标BLH的转换,以及三维速度Vxyz转换为Ven平面速度的过程。通过读取数据文件,应用XYZ2BLH和VXYZ2NEU两个辅助函数,最终将结果写入新的文件。整个转换过程中涉及到了坐标转换和矢量运算的知识。
摘要由CSDN通过智能技术生成

关于空间直角坐标系中三维速度Vxyz转Ven平面速度,需要完成连个步骤的工作:①空间直角坐标XYZ转为大地坐标BLH;②
三维速度Vxyz转Ven平面速度;
一、空间直角坐标XYZ转为大地坐标BLH
在这里插入图片描述
二、三维速度Vxyz转Ven平面速度
在这里插入图片描述
三、正式Matlab代码
(1)工程文件下设置的文件目录,data下存放输入文件,result下存生成结果文件:
在这里插入图片描述
(2)data文件夹下存放的输入文件:
在这里插入图片描述
(3)输入文件的格式,第一列为站点名,第二列至第四列为对应XYZ坐标,第五列至第六列为对应Vx,Vy,Vz速度值:
在这里插入图片描述
(4)编写Matlab代码:
1、main_Vxyz2Ven.m

clc,clear

%///读入文件数据
file_address = ['data/','raw_end.txt'];                                    %加载文件路径
%[data_name,X,Y,Z,Vx,Vy,Vz]=textread(file_address,'%s%f%f%f%f%f%f');       %第二种打开文件方式

fid = fopen(file_address, 'r');                                            %以只读的方式打开文件
if (fid  == -1 )
    disp('Error opening file');
else
    c = textscan(fid, '%s %f %f %f %f %f %f');                             % 读取文件
    fclose(fid);
    data_name=c{1};      %第一列为站点
    X=c{2};              %第二列为坐标X
    Y=c{3};              %第三列为坐标Y
    Z=c{4};              %第二列为坐标Z
    Vx=c{5};             %第三列为速度Vx
    Vy=c{6};             %第二列为速度Vy
    Vz=c{7};             %第三列为速度Vz
    
    
    N1=0;
    N1=length(X);
    %///调用函数将原XYZ坐标转为BL大地坐标
    [B,L] = XYZtoBLH(X,Y,Z);
    
    
    %///XYZ方向速度转为NE平面速度速度
    [Venu] = VXYZ2NEU(Vx,Vy,Vz,B,L,N1);
    
    
    %///转为度
    B_END=B*180/pi;
    L_END=L*180/pi;
    
    
    %///结果写入文件
    cd Result_data
    fid2=fopen('endresult.txt','w');
    for i=1:N1
        data_NAME = data_name(i);
        data_NAME = cell2mat(data_NAME);
        fprintf(fid2,' %5s',data_NAME(1:4));
        fprintf(fid2,' %8.3f',L_END(i));             %经度
        fprintf(fid2,' %8.3f',B_END(i));             %纬度
        fprintf(fid2,' %8.6f',Venu(i,2));            %Vn
        fprintf(fid2,' %8.6f',Venu(i,1));            %Ve
        if i ~= N1
            fprintf(fid2,'\n');
        end
    end
    fclose(fid2);
    cd ..
    
    disp('*******转换程序运行结束*******');
end

2、XYZ2BLH.m

function [B,L] = XYZ2BLH(X,Y,Z)
%///WGS84坐标转换到大地经纬度
v=6378137;
e2=(1.0 / 298.257223563) * (2 - (1.0 / 298.257223563));

L = atan2(Y, X);
B0 = atan(Z ./sqrt(X .* X + Y .* Y));
while 1
    N = v ./ sqrt(1 - e2 .* sin(B0) .* sin(B0));
    H = Z ./ sin(B0) - N .* (1 - e2);
    B = atan(Z .* (N+H) ./ (sqrt(X .* X + Y .* Y) .* (N .* (1-e2)+H)));
    if abs(B - B0) < 1e-6
        break;
    end
	B0 = B;
	N = v ./ sqrt(1 - e2 .* sin(B) .* sin(B));
end

end

3、VXYZ2VNEU.m

function [Venu] = VXYZ2NEU(Vx,Vy,Vz,B_end,L_end,N1)
%B_end纬度,L_end经度
for num=1:N1
    B = B_end(num);
    L = L_end(num);
    V = [Vx(num);Vy(num);Vz(num)];
    T = [-sin(L) cos(L) 0;-sin(B)*cos(L) -sin(B)*sin(L) cos(B);cos(B)*cos(L) cos(B)*sin(L) sin(B)];
    Venu(num,:)=T*V;
end
end


评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值