【MATLAB】根据SP3文件绘制星下点轨迹

【MATLAB】根据SP3文件绘制星下点轨迹

XYZ与BLH转换的代码参考：ground track.py
@Pan Li. Email:lipan.whu@gmail.com
@Jiahuan Hu. Email:hhu@whu.edu.cn

本方法通过读取精密星历中的卫星坐标进行直接的星下点轨迹绘制，并实现了对GNSS系统的选择，示例所给的SP3文件便给出了GPS、GLONASS、Galileo、BDS和QZSS，可进行系统的自由选择，代码如下所示：

%%%%根据SP3文件绘制星下点轨迹
h = geoshow('landareas.shp', 'FaceColor', 'c');
grid on
hold on
xlabel('Longitude');
ylabel('Latitude');
set(gca,'Ylim',[-90,90],'ytick',[-90:30:90]);
set(gca,'yticklabel',{'90°S','60°','30°','0°','30°','60°','90°N'});
set(gca,'Xlim',[-180,180],'xtick',[-180:30:180]);
set(gca,'xticklabel',{'180°W','150°','120°','90°','60°','30°','0°','30°','60°','90°','120°','150°','180°E'});
set(gca,'Box','on');
set(gca,'FontSize',10,'Fontname', 'Times New Roman','Fontweight', 'bold');
set(gca,'GridAlpha',1,'GridLineStyle','--');
title('BDS Satllite track','FontSize',14,'Fontweight', 'bold');

%sp3文件读取
file='igs20730.sp3';
fid = fopen(file,'rt');
while 1			   %头文件读取
   line = fgetl(fid);
   answer = findstr(line(2),'##');
   if  ~isempty(answer)
       line = fgetl(fid);
       nsat=str2num(line(4:6));
       sys = line(10:60);
       for i=1:5
            line = fgetl(fid);
            sys = [sys,line(10:60)];
       end
   end;
   for i=1:5
       ss='GRECJ';
       sns(i)=length(findstr(sys,ss(i)));
   end
   answer = findstr(line(1),'/*');
   if  ~isempty(answer), break; end;        %到头文件尾，则跳出循环
end
line = fgetl(fid);
line = fgetl(fid);
line = fgetl(fid);
for i=1:24*4
    line = fgetl(fid);
    for j=1:nsat
        line = fgetl(fid);
        X(i,j)=str2num(line(6:18));
        Y(i,j)=str2num(line(20:32));
        Z(i,j)=str2num(line(34:46));
        [B0,L0]=XYZtoBLH(X(i,j),Y(i,j),Z(i,j));
        B(i,j)=B0;
        L(i,j)=L0;
    end
end

%各系统对应卫星列数
G=1:sns(1);
R=sns(1)+1:sns(1)+sns(2);
E=sns(1)+sns(2)+1:sns(1)+sns(2)+sns(3);
C=sns(1)+sns(2)+sns(3)+1:sns(1)+sns(2)+sns(3)+sns(4);
J=sns(1)+sns(2)+sns(3)+sns(4)+1:sns(1)+sns(2)+sns(3)+sns(4)+sns(5);

% for j=[G,R,E,C]
for j=[C]  %%%%在此选择GNSS系统（GRECJ）%%%%
    for i=1:95
        if L(i,j)-L(i+1,j)<6&& L(i,j)-L(i+1,j)>-6
           plot(L(i:i+1,j)*57.3,B(i:i+1,j)*57.3,'MarkerSize',8,'linewidth',2,'color',[18 53 85]/255);
           hold on
       end
    end
end

%WGS84坐标转换到大地经纬度
function [b,l] = XYZtoBLH(x,y,z)
v=6378137.0;
e2=(1.0/298.257223563)*(2-(1.0/298.257223563));
r2=x*x+y*y;
b=atan(z/(sqrt(r2)));
l=atan2(y,x);
h=sqrt(r2+z*z)-v;
end

若选择的系统为北斗，则最终绘制的结果如下：

若采用m_map提供的相关底图，可得到类似如下结果：