诚聘高新网站开发工程师山西百度推广开户
0. 码字不易，点赞加关注（公众号：WZZHHH，部分资料在公众号可以下载），使用请注明出处（根据我的研究方向，我会不断更新代码）。
1.计算PWV的方法一般采用有三种，
GNSS ERA5、ERA5 PWV和探空站PWV(RS PWV)进行相互比较。
2.GNSS PWV的获取
ZTD我用gamit计算得到，ZHD需要气压用Saastamoinen模型计算，该模型需要压力和站位信息。理想情况下，根据积分法，可以通过探空站的大气剖面得到Tm或者ERA5数据积分。因此，压力、温度和Tm是GNSS PWV检索的关键参数。
3.我的压力、温度和Tm获取方法
本文使用ERA5数据集通过水平和垂直方向的补偿得到GNSS站点垂直方向的ERA5数据集，得到GNSS站点处的温度、气压和Tm，并积分得到该位置的ERA5 PWV。
代码结果如图所示
4.ERA5的下载方法在前面的博文，自行学习下载
ERA5数据下载单层和气压层数据（温度、相对湿度、位势等）_gee怎么下载era5的nc数据-CSDN博客
5.下面是我的部分代码，具体获取请关注公众号：WZZHHH或者咸鱼：WZZHHH123
% nc提取的t、z、r、q等都是四维数据，时间是最后一维Edata(:,1) = T(:,:,:,p);Edata(:,3) = R(:,:,:,p);% 把位势z转化为海拔Edata(:,2) = Z(:,:,:,p)./9.80665;Edata(:,2) = 6371229.*Edata(:,2)./(6371229-Edata(:,2));% 得到站点的数据（两种情况：最底层高于GNSS站和最底层低于GNSS站）[g,Edatanew,GPS_T,GPS_pr,e] = ERAdata_T_P_PWV(Height,Edata,lat);% 由GNSS站点上方数据集（Edatanew）积分得到PWV和Tm[pwer,Tmm] = integral_pwv_Tm(Edatanew,g);% 得到站点的数据按照列是站点，行是时间进行排列ERAT(p,P) = GPS_T; ERAPR(p,P) = GPS_pr;pwera(p,P) = pwer.*1000; % 乘1000是mm Tm(p,P) = Tmm;es(p,P)=e;disp([List(I).name '的第' num2str(P) '列的第' num2str(p) '行'])endend% 将每个nc得到的数据纵向排列ERA_e = [ERA_e;es]; % 水汽压ERA_T = [ERA_T;ERAT]; % 温度ERA_P = [ERA_P;ERAPR]; % 气压ERA_pwv = [ERA_pwv;pwera]; % ERA5得到的pwvERA_Tm = [ERA_Tm;Tm]; % ERA5得到的Tm% 每完成一个，输出一次名字，便于找到aa = strcat(filen(1,end-8:end-3),' OK');disp(aa);