matlab读取gamit解算的met文件，获取ZTD、ZWD、PWV、ZHD等数据（附matlab代码）

0. 码字不易，点赞加关注（公众号：WZZHHH，后期代码将在公众号可以下载），使用请注明出处（根据我的研究方向，我会不断更新代码）。

00.提供gamit安装服务，可咨询咸鱼：WZZHHH23（咸鱼平台更放心）

1. 读取gamit解算的met文件，根据前缀met_读取相应文件下met文件，met文件的列分别为1为year、2为doy（年积日）、3为小时、4为分钟、5为秒、6为ZTD、7为ZWD、8不知道、9为PWV、11为气压、12为温度、13为ZHD。

2. 下图为具体的met文件，对照上面列，找到自己需要的。

3. 我按照年积日存放的met文件，存放方式不同的，记得修改代码。我的存放方式如下图。

4. 下面是我的代码，可能没大佬写得好，凑合看吧，能用。

我用的ZTD数据是第六列，其他数据可在第84行更改列即可。

function load_ZTD_file
% 读取gamit解算的met文件
% 我按照年积日存放的met文件
% 存放方式不同的，记得修改
% 根据前缀met—读取相应文件下ztd文件清单，再根据需要更改列
% 我用的ZTD数据是第六列。其他数据可在第84行更改
% met文件的列分别为1为year、2为doy（年积日）、3为小时、4为分钟、5为秒、
% 6为ZTD、7为ZWD、8不知道、9为PWV、11为气压、12为温度、13为ZHD。

clc,clear;
% met文件路径，2021年数据为例
ypath = 'D:\paper_write\ZTD\2021\';  
% 提取数据后，存放路径
save_path= 'D:\paper_write\paper_code\2\mat\ztd\2021ztddata.mat';
%% ---------------------找到存在met的站点名字---------------------------

% 遍历文件夹，看总共有几个站点
ydata = dir(ypath);

station = {};  % station 站点
for i=3:size(ydata,1)

    % 遍历年积日文件夹
    dpath = [ypath ydata(i).name '\'];
    ddata = dir(dpath);
    for ii=3:size(ddata,1)

        % i=3是第一天，第一天初步建立站点名称
        if i==3
            station = [station;{ddata(ii).name(5:8)}];
        else

            % 后面增加的站点加入到station中
            if ~contains(ddata(ii).name(5:8),station)
                station = [station;{ddata(ii).name(5:8)}];
            end
        end
    end

end

%% -----------------------------提取gamit解算的met文件------------------

% 这个部分本可以和上面合在一起，为了便于理解，，所以分开写

ztd_all = {};
parfor i = 3:size(ydata,1)
    dpath = [ypath ydata(i).name '\'];

    % 避免文件夹有其他文件，只读取met_开头
    List = dir(fullfile(dpath,'/','met_*'));

    for I = 1:size(List,1)

        % 站点名称
        stn = List(I).name(5:8);

        % 读取met文件
        filename = [dpath,List(I).name];
        fid = fopen(filename,'r');

        % 提取met文件的ztd
        ztd = [];

        % 逐行读取
        while ~feof(fid)
            line = fgetl(fid);

            % 依次读取每行，met文件带*的删除，剩下就是想要的数据
            if ~contains(line,'*')
                ztd = [ztd;str2num(line)]; %#ok<ST2NM>
            end
        end

        % 我使用的是1小时时间分辨率
        % 避免个别时间可能会缺失，缺少数据为nan
        ztd = hours(ztd);
        fclose(fid);

        % cellstr：char转cell、num2cell：double转cell
        stn = repmat(cellstr(stn),24,1);

        % 第六列是ztd，我只需要这个
        ztd_all = [ztd_all;[stn num2cell([ztd(:,1:3) ztd(:,6)])]];
    end
end
% ztd_2：把ztd相同测站放在一起，相同测站按照日期排列更好计算
% ztd_2就是排列好的ztd_all

% 按照station名称排列站点数据
ztd = [];
for I =1:size(station)
    flag = contains(ztd_all(:,1),(station(I,:)));
    ztd_1 = ztd_all(flag,:);
    ztd = [ztd;ztd_1];
end

save(save_path,'ztd');
end

function [ztd] = hours(ztd)
% 25行证明数据齐全，没有的为nan
if size(ztd,1)==25
    % 第25行为下一天的0时，对其删除
    ztd = ztd(1:24,:);
else

    % 年月日时，后面为0的数组
    ztdz = nan(24,14);
    yearz = repmat(ztd(1,1:2),24,1);
    hour = (0:1:23).';
    ztdz = [yearz hour ztdz];

    % 将存在的数据进行填充
    % ztd第3列为小时0-23，把存在数据的小时填充到ztdz中
    for i=1:size(ztd,1)
        ztdz(ztd(i,3)+1,:) = ztd(i,:);
    end
    ztd = ztdz;
end

end