gamit生成的met文件的处理方式：2.找到未处理的站点

1.前文讲解了分类met文件的gamit生成的met文件的处理方式：1.移动到年积日文件夹下-CSDN博客）

2.这里将分类的met文件进行识别，找到未处理的文件，并将其复制到相应文件夹下，为后期gamit处理缺失数据做准备。

clc,clear;

%%  ---------------------改此处即可--------------------------------
year = '2021'; % 数据年份
ztd_path = ['D:\paper_write\ZTD\',year,'\'];     % met文件地址
igs_path = ['F:\A-GNSS数据\B-IGS\',year,'\'];    % 精密星历地址
brdc_path= ['F:\A-GNSS数据\B-brdc\',year,'\'];   % 广播星历地址
d_path   = ['F:\A-GNSS数据\B-陆态网\',year,'\']; % d文件地址（o文件也可）
err_path = 'D:\paper_write\data\';               % 缺失met数据存放地址

%%  --------------------------------------------------------------


ztd = dir(ztd_path);
for i = 3:size(ztd,1)

    % met和d文件的年积日文件夹
    ztdday_path = [ztd_path,ztd(i).name,'\'];
    dday_path = [d_path,ztd(i).name,'\'];

    % 有d文件，没对应的met文件，则这个站点就是没处理的站点
    % 逐个d文件进行判断
    dfile = dir(dday_path);
    zfile = dir(ztdday_path);

    % 将met文件的站点排列
    z_name={};
    for ii = 3:size(zfile)
        z_name = [z_name;{zfile(ii).name}];
    end

    for ii = 3:size(dfile)

        dfile_name = dfile(ii).name(1:4); % 站点名称
        dfile_day = dfile(ii).name(5:7);  % d文件的年积日

        % 根据d文件名称组合成met文件的名称，判断是否有一样的文件存在
        dname = ['met_',dfile_name,'.',year(3:4),dfile_day];

        
        
        % 如果不存在met文件，将移动对应的d文件，精密星历和广播星历
        % 为再次解算做准备
        if ~contains(dname,z_name)
            % 复制brdc文件
            copyfile([brdc_path,'brdc',dfile_day,'0.',year(3:4),'n'],err_path);

            % 复制igs文件
            gps_week = ymd2gpsweek(year,dfile_day); % 将年积日转化为gps周
            copyfile([igs_path,'igs',gps_week,'.sp3'],err_path);

            % 复制d文件（o文件更好，也不用转换了）
            copyfile([dday_path,dfile(ii).name],err_path);
        end
    end


end
%% -----------------------------函数------------------------------------

% 年月日转换为GPS周
function gps_week = ymd2gpsweek(year,doy)

year = str2double(year);
doy = str2double(doy);
[year,mon,day]=doy2ymd(year, doy);

mjd = ymd2mjd(year, mon, day);
week = floor((mjd - 44244.0) / 7.0);
day = floor(mjd - 44244.0 - week * 7.0);
gps_week = [num2str(week),num2str(day)];

end

% 年积日转换为年月日
function [year,mon,day]=doy2ymd(year, doy)
day = doy;
mon = 0;
monthdays = [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31];
for i=1:13
    monthday = monthdays(i);
    if i==3 && leapyear(year) == 1
        monthday = monthday + 1;
    end
    if day > monthday
        day = day - monthday;
    else
        mon = i - 1;
        break
    end
end
end

% 年月日转化为儒略日
function mjd=ymd2mjd(year, mon, day)
if mon <= 2
    mon = mon + 12;
    year = year - 1;
end
mjd = 365.25 * year - rem(365.25 * year , 1.0) - 679006.0;
mjd = mjd + floor(30.6001 * (mon + 1)) + 2.0 - ...
    floor(year / 100.0) + floor(year / 400) + day;
end