GNSS与探空站不并址的解决（文末代码和压缩包分享，自行下载）

Part1   这个是免费的，点个赞来点动力，不然没动力啊，铁铁。

里面有文章出处，代码在最后，还有适合小白的压缩包分享，（纯手敲，有个别错字请担待，语文从小不好）。

Part2   有关Matlab 获取代码关注公众号WZZHHH，或者咸鱼关注：WZZHHH123

对IGRAv2进行质量控制得到PWV和Tm的matlab代码；

使用Matlab计算IGRAv2探空站的Tm和PWV；

提取探空站IGRAv2全部数据；

ERA5 解算合集（温度、气压、PWV、水汽压和 Tm）代码获取；

远程安装gamit等服务，博主乐于助人，心地善良，但是不可能免费滴。

Part3   探空站和GNSS站点不并址

根据论文（DOI: 10.1175/JCLI-D-16-0591.1），解决GNSS和探空站不并址如何对比经度的问题。

对于这种比较：

（1）要求水平间距小于50 km

（2）高程差小于200 m。

Part4   如何使用Matlab实现

下面是我随便选探空站。

这是选的GNSS站

根据经纬度计算每个GNSS站与探空站的距离，再计算其高度差，匹配合适的站点即可。运行代码可以获得如下图所示的匹配站点。

按照【GNSS站ID、探空站ID、距离】存为excel表格。

Part5   Matlab代码展示

如下图是压缩包的内容，方便小白学习（虽然你是小白，但是请尊重我，别什么都问我【譬如下载数据这些百度都能解决的】，信不信我顺着网线刀你）。

压缩包地址再公众号的文章（WZZHHH）获取，免费的哦。公众号回复：GNSS与探空站不并址，即可得到

下面是代码展示：

% 基于matlab2022写的
% 微信公众号和CSDN：WZZHHH
% 根据经纬度，计算GNSS和探空站之间的距离
% 并根据论文（DOI: 10.1175/JCLI-D-16-0591.1）
% 解决GNSS和探空站不并址的对比问题
% 1.水平间距小于50 km
% 2.高程差小于200 m

clc,clear;
% 根据经纬度，计算两个点之间的距离
% GNSS站点的经纬度高程

[station_gnss,name_gnss] = xlsread('GNSS站点的ID经纬度高程.xlsx');
[station_sound,name_sound] = xlsread('探空站的经纬度高程');

suitable_station={};

for ii = 1:size(station_sound,1)

    % 探空站ID、经纬度、高度
    sound_id = name_sound{ii,1};
    lat1 = station_sound(ii,1);
    lon1 = station_sound(ii,2);
    h_sound = station_sound(ii,3);

    for i = 1:size(station_gnss,1)
        % 经纬度
        lat = station_gnss(i,1);
        lon = station_gnss(i,2);
        
        % 计算两个站之间的距离
        [c1km,c2km]=distances([lat lon], [lat1 lon1]);
        dist = (c1km+c2km)/2;

        % i=1时候，最近的为第一个站点
        % n为距离最近的GNSS站在所有GNSS站的所在行
        if i==1
            mindise=dist;
            gnss_id = name_gnss{i,1};   % GNSS的站点名
            h_gnss = station_gnss(i,3); % GNSS高度
            n=i;
        else
            % i~=1时，dist小于上一轮最小值时，最近值为dist
            if dist<mindise
                mindise=dist;
                gnss_id = name_gnss{i,1};   % GNSS的站点名
                h_gnss = station_gnss(i,3); % GNSS高度
                n=i;
            end
        end
    end

    if mindise<50
        if h_sound-h_gnss<=200&&h_sound-h_gnss>=-200
            % 【'GNSS站点名','探空站','两者距离'】
            suitable_station=[suitable_station;{gnss_id,sound_id,num2str(mindise)}];
        end
    end
clear dist;
end
header = {'GNSS','探空站','两者距离'};
xlswrite('GNSS和探空站距离和高度满足的点.xlsx',[header;suitable_station])

function [d1km,d2km]=distances(latlon1,latlon2)
radius=6371;
lat1=latlon1(1)*pi/180;
lat2=latlon2(1)*pi/180;
lon1=latlon1(2)*pi/180;
lon2=latlon2(2)*pi/180;
deltaLat=lat2-lat1;
deltaLon=lon2-lon1;
a=sin((deltaLat)/2)^2 + cos(lat1)*cos(lat2) * sin(deltaLon/2)^2;
c=2*atan2(sqrt(a),sqrt(1-a));
d1km=radius*c;    %Haversine distance

x=deltaLon*cos((lat1+lat2)/2);
y=deltaLat;
d2km=radius*sqrt(x*x + y*y); %Pythagoran distance
end