1.文中适用精度对比GNSS PWV,采用ERA5相对湿度、气压、温度以积分方法得到PWV。
2.也可用于ERA5插值温度、气压的精度对比。
2.由于GNSS与ERA5并不并址,需要进行水平和垂直方向的补偿。
3.水平方向补偿:在最接近GPS站经纬度坐标的4个网格点,利用反距离加权水平内插到每个水平压力级的GPS站位置 ,得到GPS测站点垂直剖面数据。
4.接下来,在垂直方向上将气象数据内插(外推)到每个GPS天线位置上。
5.下面是部分代码:
% 利用ERA5气压层数据积分计算PWV
%ERA5数据使用位势、温度、相对湿度
clc;clear;
% 加载站点的经纬度
[station_x,~] = xlsread('F:\predict_pm\xls\air_surface.xlsx');
List = dir('F:\predict_pm\ERA5\bj\allyear\leves\');%气压层数据文件夹
filename = 'F:\predict_pm\result\ERA5_data\21年数据.mat';%生成文件存放位置
% 经纬度边界(左上角点的位置)ERA5边界,ERA5分辨率0.25.
lon_t = 114.5; lat_t = 41.25;
resolution = 0.25;
ERA_PR=[];ERA_T=[];pwvera5=[];esat=[];
for I = 3:size(List,1)
% z:位势、t:温度、r:相对湿度
filen = fullfile(List(I).folder,List(I).name);
year = str2num(List(I).name(1:4));
a = ncinfo(filen);
%读取时间
Time_num = ncread(filen,'time');
times = length(Time_num);
%读取经纬度
longitude = ncread(filen,'longitude');
latitude = ncread(filen,'latitude');
%判断边界是否正确
if longitude(1,1)~=lon_t||latitude(1,1)~=lat_t
warning('边界出错!!!');
end
% 读取位势、温度、相对湿度、气压层
z = ncread(filen,'z');
t = ncread(filen,'t');
r = ncread(filen,'r');
levels = ncread(filen,'level');
% 提前分配空间
pwera = zeros(times,size(station_x,1));
ERAPR = zeros(times,size(station_x,1));
ERAT = zeros(times,size(station_x,1));
rh = zeros(times,size(station_x,1));
es = zeros(times,size(station_x,1));
% 根据站点数量分别计算
for P = 1:size(station_x,1)
% 读取站点经纬度、海拔
lon = station_x(P,1);
lat = station_x(P,2);
Height = station_x(P,3);