read_o.m

可以用法来读取全球定位系统GPS N文件和O文件的程序

function [OHeader,ObsData]=Read_O(q)
%  global DateT 
%OHeader.
     %interval     历元间隔
     %TimeOBeg       历元开始时间
     %TimeOEnd       历元结束时间
     %ApproXYZ.X  ApproXYZ.Y  ApproXYZ.Z  测站的近似坐标
     %Ant_H  Ant_E  Ant_N   天线相位中心的偏移量
     %NumOSignal     接收到的卫星信号个数
     %TypeOSignal       接收到的卫星信号的类型
%ObsData(i).        i表示历元
     %TimeObs      观测历元的时间
     %Int_timeObs  观测历元的时间（取整）
     %SatNum       接收到的卫星数
     %SatPRN(k)    卫星随机码    k表示其维数，其值等SatNum
     %DataL1(k)       L1数据
     %DataL2(k)       L2数据
     %DataC1(k)       C1数据
     %DataP1(k)       P1数据
     %DataP2(k)       P2数据
 
     
%打开O文件
[filename,pathname]=uigetfile('*.??O','打开O文件');
fid=fopen(strcat(pathname,filename),'rt');
if fid==-1
    msgbox('The filename or pathname is not correct','warning','warn');
    return;
end

%读取O文件的头文件
OHeader.interval=0;
OHeader.TimeOBeg=0;
OHeader.TimeOEnd=0;

while (1)
    line=fgets(fid);
    if (line==-1)
        break;
    end
    if (findstr(line,'APPROX POSITION XYZ')~=0)
        Temp=sscanf(line,'%f%f%f');
        OHeader.ApproXYZ.X=Temp(1);
        OHeader.ApproXYZ.Y=Temp(2);
        OHeader.ApproXYZ.Z=Temp(3);
    end
    if (findstr(line,'DELTA H/E/N')~=0)
        Temp=sscanf(line,'%f%f%f');
        OHeader.Ant_H=Temp(1);
        OHeader.Ant_E=Temp(2);
        OHeader.Ant_N=Temp(3);
    end
    if (findstr(line,'TYPES OF OBSERV')~=0)
        OHeader.NumOSignal=str2num(line(5:6));
        for i=1:OHeader.NumOSignal
            ii=11+(i-1)*6;
             if (line(ii:ii+1)=='L1')   %0-L1;1-L2;2-C1;3-P1;4-P2;5-D1;6-D2;
                OHeader.TypeOSignal(i)=0;
            elseif ((line(ii:ii+1)=='L2'))
                OHeader.TypeOSignal(i)=1;
            elseif ((line(ii:ii+1)=='C1'))
                OHeader.TypeOSignal(i)=2;
            elseif ((line(ii:ii+1)=='P1'))
                OHeader.TypeOSignal(i)=3;
            elseif ((line(ii:ii+1)=='P2'))
                OHeader.TypeOSignal(i)=4;
            elseif ((line(ii:ii+1)=='D1'))
                OHeader.TypeOSignal(i)=5;
            elseif ((line(ii:ii+1)=='D2'))
                OHeader.TypeOSignal(i)=6;
            end
        end
    end
    if (findstr(line,'INTERVAL')~=0)
        OHeader.interval=str2num(line(5:11));
    end
    if (findstr(line,'TIME OF FIRST OBS')~=0)
        Temp=str2num(line(3:42));
        OHeader.TimeOBeg=Time_MJDtime(Temp(1),Temp(2),Temp(3),Temp(4),Temp(5),round(Temp(6)));
    end
    if (findstr(line,'TIME OF LAST OBS')~=0)
        Temp=str2num(line(3:42));
        OHeader.TimeOEnd=Time_MJDtime(Temp(1),Temp(2),Temp(3),Temp(4),Temp(5),round(Temp(6)));
    end
    if (findstr(line,'END OF HEADER')~=0)
        break;
    end
end

%读取观测数据
num_Epoch=1;
MJD_1s=(Time_MJDtime(2007,12,01,00,00,01)-Time_MJDtime(2007,12,01,00,00,00));
while (1)
    
    %读取第一行
    line=fgets(fid);
    if (line==-1)
        break;
    end
    Temp=str2num(line(2:32));
     
    if(length(Temp)~=8)
        break;
    end
    if(Temp(1)>80)
        Temp(1)=Temp(1)+1900;
    else
        Temp(1)=Temp(1)+2000;
    end
    Temp_dclock=str2num(line(69:end));
    if length(Temp_dclock)==0
        ObsData(num_Epoch).TimeObs=Time_MJDtime(Temp(1),Temp(2),Temp(3),Temp(4),Temp(5),Temp(6));
    else
        ObsData(num_Epoch).TimeObs=Time_MJDtime(Temp(1),Temp(2),Temp(3),Temp(4),Temp(5),Temp(6))-Temp_dclock*MJD_1s;
    end
    ObsData(num_Epoch).Int_timeObs=Time_MJDtime(Temp(1),Temp(2),Temp(3),Temp(4),Temp(5),round(Temp(6)));
    ObsData(num_Epoch).SatNum=Temp(8);
    for i=1:ObsData(num_Epoch).SatNum
        ii=31+i*3;
        ObsData(num_Epoch).SatPRN(i)=str2num(line(ii:ii+1));
    end
    if(length(line)<70)
        DateT=0.0;
    else
        DateT=str2num(line(69:end)); 
    end
    
    %读取下面几行
     if (OHeader.NumOSignal>5)  
        for i=1:ObsData(num_Epoch).SatNum
            line=fgets(fid);
            line1=fgets(fid);
            line=[line line1];
            Temp=Read_Epoch_Obs(OHeader.NumOSignal,OHeader.TypeOSignal,line,DateT); 
            ObsData(num_Epoch).DataL1(i)=Temp(1);
            ObsData(num_Epoch).DataL2(i)=Temp(2);
            ObsData(num_Epoch).DataC1(i)=Temp(3);
            ObsData(num_Epoch).DataP1(i)=Temp(4);
            ObsData(num_Epoch).DataP2(i)=Temp(5);
        end
    else
        for i=1:ObsData(num_Epoch).SatNum  
            line=fgets(fid);
            Temp=Read_Epoch_Obs(OHeader.NumOSignal,OHeader.TypeOSignal,line,DateT);  
            ObsData(num_Epoch).DataL1(i)=Temp(1);
            ObsData(num_Epoch).DataL2(i)=Temp(2);
            ObsData(num_Epoch).DataC1(i)=Temp(3);
            ObsData(num_Epoch).DataP1(i)=Temp(4);
            ObsData(num_Epoch).DataP2(i)=Temp(5);
        end
    end
    num_Epoch=num_Epoch+1;
end

%SiteData.ObsData=ObsData;
fclose(fid);