02.20 科研笔记：涡旋数量统计

原始数据：

FileZilla软件批量下载

ftp://ftp-access.aviso.altimetry.fr/duacs-experimental/dt-phy-grids/altimetry_arctic/version_02_00/multimission

G:\multimission

科研工具：

1. 涡旋的探测和追踪——Python

文件名及路径：G:\Aviso_multi_SLA\eddy_detect_recent.py； G:\Aviso_multi_SLA\eddy_track_recent.py；

G:\Aviso_multi_SLA\read_eddy_detect.py；G:\Aviso_multi_SLA\read_eddy_track.py

2.  涡旋统计——matlab

G:\Aviso_multi_SLA\number_statistics_2011_19.m

3. 可视化——GMT

在涡旋完成探测和追踪的基础上，我们将python中涡旋轨迹数据储存成npz，为了在matlab中读取数据，又将npz转为mat。

''' npz to mat '''
import scipy.io as io
io.savemat('eddy_track_2011——2019_Canada Basin.mat', mdict=np.load('eddy_track_2011——2019_Canada Basin.npz',allow_pickle=True))

将研究区域划分为的网格区域，对落在每个网格区域的涡旋数量、半径和振幅进行统计，最终得到涡旋的空间分布特征。具体操作函数如下：

function [xx,yy,M,N,Z_number,Z_R,Z_amp]=eddynumberstat(eddies,aa)
%% 读取涡旋半径、时间范围、经纬度信息，以及涡旋点的经纬度；半径路径和涡旋极性；
for i = 1:length(eddies)
    lon(i).lon=eddies{i}.lon;
    lat(i).lat=eddies{i}.lat;
    radius(i).radius=eddies{i}.scale;
    amp(i).amp=eddies{i}.amp;
    age(i).age=eddies{i}.age;
    type(i).type=eddies{i}.type;
end  

ind = zeros(length(eddies),1);
for j = 1:length(eddies)
    if(isempty(aa(j).rg))
        ind(j) = 0;
    else
        ind(j) = 1;
    end
end

%% find the begin point of the eddy 
ind1 = find(ind==1);
for kk = 1:length(ind1)
    ll_eddy(kk,1) = lon(ind1(kk)).lon(1);
    ll_eddy(kk,2) = lat(ind1(kk)).lat(1);
    r_eddy(kk).rg = radius(ind1(kk)).radius(1);
    amp_eddy(kk).rg = amp(ind1(kk)).amp(1);
%     mean_radius(kk,1) = mean(r_eddy(kk).rg);
end

%% Z region
xx = -180:1:-110;
yy = 65:1:83;
[xx,yy] = meshgrid(xx,yy);
[M,N] = size(xx);
Z_number = zeros(M,N);
Z_R = zeros(M,N);Z_amp = zeros(M,N);
for k = 1:length(ind1)
    lx = ll_eddy(k,1);
    ly = ll_eddy(k,2);
    for i = 1:19
        for j = 1:71
            if(abs(lx-xx(i,j))<0.5&abs(ly-yy(i,j))<0.5)
                Z_number(i,j) = Z_number(i,j)+1;
                Z_R(i,j) = Z_R(i,j) + r_eddy(k).rg; %% summation
                Z_amp(i,j)=Z_amp(i,j) + amp_eddy(k).rg;
            end
        end
    end
    disp(k)
end

end

主程序如下：

load('eddy_track_2011——2020_Canada Basin.mat');

%% 读取涡旋半径、时间范围、经纬度信息，以及涡旋点的经纬度；半径路径和涡旋极性；
for i = 1:length(eddies)
    lon(i).lon=eddies{i}.lon;
    lat(i).lat=eddies{i}.lat;
    radius(i).radius=eddies{i}.scale;
    amp(i).amp=eddies{i}.amp;
    age(i).age=eddies{i}.age;
    type(i).type=eddies{i}.type;
end  

%% 确定目标区域并检索
aim_eddy_area=[-180 -110;65 83];
for i = 1:length(eddies)
    aa(i).rg = find(lon(i).lon>aim_eddy_area(1,1)&lon(i).lon<aim_eddy_area(1,2)&lat(i).lat>aim_eddy_area(2,1)&lat(i).lat<aim_eddy_area(2,2)&age(i).age>5);    
    bb(i).rg=find(lon(i).lon>aim_eddy_area(1,1)&lon(i).lon<aim_eddy_area(1,2)&lat(i).lat>aim_eddy_area(2,1)&lat(i).lat<aim_eddy_area(2,2)&age(i).age>5&strcmp(type(i).type,'anticyclonic'));
    cc(i).rg=find(lon(i).lon>aim_eddy_area(1,1)&lon(i).lon<aim_eddy_area(1,2)&lat(i).lat>aim_eddy_area(2,1)&lat(i).lat<aim_eddy_area(2,2)&age(i).age>5&strcmp(type(i).type,'cyclonic'));
    disp(i)
end

[xx,yy,M,N,Z_number,Z_R,Z_amp]=eddynumberstat(eddies,aa);
[xx,yy,M,N,Z_number_a,Z_R_a,Z_amp_a]=eddynumberstat(eddies,bb);
[xx,yy,M,N,Z_number_c,Z_R_c,Z_amp_c]=eddynumberstat(eddies,cc);

代码解释：

函数中输入参数eddies为涡旋轨迹数据，aa为筛选条件

GMT可视化结果如下：