原始数据:
FileZilla软件批量下载
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为筛选条件