
clear all;
close all;
clc;

%%

folder_path = '/Volumes/Qiao/cyclone_track_era5/track_density/mat_files/';
file_list = dir(fullfile(folder_path, '*syno.mat'));
file_list0 = dir(fullfile(folder_path, '*meso.mat'));

%------------------读SLP----------------
datapath = ['/Volumes/Qiao/cyclone_track_era5/era5_interp2erai/input_era5/pmsl_2007.era5.nc'];  % 构建文件名，假设文件名格式为 data_YYYY.nc
lon0=ncread(datapath,'longitude');
lat0=ncread(datapath,'latitude');
lon=repmat(lon0,[1 721]);
lat=repmat(lat0',[1440 1]);
lon(lon < 0) = lon(lon < 0) + 360;
time=ncread(datapath,'time');
msl=ncread(datapath,'msl');%Pa 'air_pressure_at_mean_sea_level'
refTime = datetime(1900, 1, 1, 0, 0, 0);
realTime = refTime + hours(time);
%MSL=msl(:,:,904);%2007-8-14--18:00
MSL=msl(:,:,905);%2007-8-15--00:00
clear msl time refTime lon0 lat0

%----------------读winds-----------------
file=['/Volumes/Qiao/ERA5/win_winds/era5_uwind_win_',num2str(2007),'.nc'];
filev=['/Volumes/Qiao/ERA5/win_winds/era5_vwind_win_',num2str(2007),'.nc'];
u0=ncread(file,'u10');
v0=ncread(filev,'v10');
Lon0=ncread(file,'longitude');
Lat0=ncread(file,'latitude');
Lon=double(repmat(Lon0,[1 121]));
Lon(Lon < 0) = Lon(Lon < 0) + 360;
Lat=double(repmat(Lat0',[601 1]));
Time=ncread(file,'time');
% 转换ERA5时间为MATLAB的datetime格式   输出转换后的真实日期和时间
refTime = datetime(1900, 1, 1, 0, 0, 0);
RealTime = refTime + hours(Time);
%--------选%2007-8-14--18:00
% u=u0(:,:,1795);
% v=v0(:,:,1795);
%2007-8-15--00:00
u=u0(:,:,1801);
v=v0(:,:,1801);
clear u0 v0 Lon0 Lat0 Time

%% -------------------------画图------------------------------------
%cdata=rgb2mat('/Users/wangxiaoqiao/Downloads/software/MATLAB/rgb_files/cmp_flux.rgb',50);
cdata=rgb2mat('/Users/wangxiaoqiao/Downloads/software/MATLAB/rgb_files/GMT_no_green.rgb',10);

str0=["2003","2004","2005","2006","2007","2008","2009","2010","2011","2012","2013","2014","2015","2016","2017","2018","2019"];
img=gcf;%设置画布figure的大小，使其在屏幕上的显示位置是以（0.2，0.2）为原点，长0.64，宽0.32。同gca一样，仍然是左边界，下边界为0，上边界，右边界为1。
set(gcf,'unit','normalized','position',[0.01,0.01,0.8,0.98]);
%[ha, pos] = tight_subplot(Nh, Nw, gap, marg_h, marg_w);设置子图的间距
%（行数，列数，[子图间的高度间隔 子图间的宽度间隔],[整个图形距底边的高度间距 距顶边高度间距],[整个图形距左边的间距 距右边的间距])
ha = tight_subplot(5,4,[0.036 -0.26],[.10 .04],[.0001 .001]);

%创建一个a-i的数组，便于命名
str=["a","b","c","d","e","f","g","h","i","j","k","l","m","n","o","p","q","r"];

for i = 1:17
    subplot(ha(i));
    file_name = file_list(i).name;
    file_name0 = file_list0(i).name;
    
    % 执行load函数
    data = load(fullfile(folder_path, file_name));
    data0 = load(fullfile(folder_path, file_name0));
    
    % 可以在此处对数据进行处理
    xc=data.xc;yc=data.yc;density=data.density;
    xc0=data0.xc;yc0=data0.yc;density0=data0.density;
    
    Density=density+density0;
    Density(Density==0) = nan;
    
    %计算空间累积的数值
    x1=find(xc>=150 & xc<=220 & yc>=-90 & yc<=-60);
    lon1=xc(x1);
    lat1=yc(x1);
    
    lonx1=[170 189 189 170 170];
    laty1=[-84 -84 -78 -78 -84];
%     lonx1=[170 205 205 170 170];
%     laty1=[-84 -84 -78 -78 -84];
    in1=inpolygon(lon1,lat1,lonx1,laty1);%in=1表示该点在多边形区域内部
    in1=double(in1);
    %对应区域做空间累积
    tmp1=Density(x1);
    tmp2=in1.*tmp1;
    y1=find(tmp2==0);
    tmp2(y1)=[];
    density_sum=nansum(tmp2);
    %
    
    m_proj('Miller', 'long',[159 210], 'lat',[-84 -69]);hold on;
    m_pcolor(xc,yc,Density);hold on
    m_usercoast('/Users/wangxiaoqiao/Downloads/software/MATLAB/topodata/whole_Antarctic/topodata.mat','linewi',2,'color',[.4 .4 .4]);hold on
    
    
    if i==1|i==9|i==13
        m_grid('linewi',1,'tickdir','in','tickstyle','dd','xtick',160:15:210,'xticklables',[],...
            'ytick',-84:5:-69,'fontsize',14, 'ticklen',0.007,'FontName','sans serif');hold on;
    elseif i==5
        m_grid('linewi',1,'tickdir','in','tickstyle','dd','xtick',160:15:210,'xticklables',[],...
            'ytick',-84:5:-69,'fontsize',14, 'ticklen',0.007,'FontName','sans serif');hold on;
        m_line(lonx1,laty1,'Color', 'r', 'LineWidth', 2); hold on;
        
    elseif i==2|i==3|i==4|i==6|i==7|i==8|i==10|i==11|i==12
        m_grid('linewi',1,'tickdir','in','tickstyle','dd','xtick',160:15:210,'xticklables',[],...
            'ytick',-84:5:-69,'yticklables',[],'fontsize',14, 'ticklen',0.007,'FontName','sans serif');hold on;
    elseif i==14|i==15|i==16
        m_grid('linewi',1,'tickdir','in','tickstyle','dd','xtick',160:15:210,...
            'ytick',-84:5:-69,'yticklables',[],'fontsize',14, 'ticklen',0.007,'FontName','sans serif');hold on;
    elseif i==17
        m_grid('linewi',1,'tickdir','in','tickstyle','dd','xtick',160:15:210,...
            'ytick',-84:5:-69,'fontsize',14, 'ticklen',0.007,'FontName','sans serif');hold on;
        
    end
    
    m_text(160,-67.5,str(i),'Fontsize',16,'Color','k','FontWeight','bold',...
        'HorizontalAlignment','left','FontName','sans serif');hold on
    box on
    
    colormap(cdata);
    caxis([0 12]);
    m_text(178,-67.5,str0(i),'Fontsize',15,'Color','k',...
        'HorizontalAlignment','left','FontName','sans serif');hold on
    m_text(202,-67.5,['(' num2str(density_sum) ')'],'Fontsize',15,'Color','r',...
        'HorizontalAlignment','left','FontName','sans serif','FontWeight','bold');hold on
end
    
colormap(cdata);
cb=colorbar('h');
caxis([0 10])
set(cb,'YTick',0:2:10,'position',[0.54,0.16,0.32,0.016],'tickdir','in','ticklen',0.007,'Fontsize',15,'FontName','sans serif');
set(get(cb,'title'),'string','Cyclone track density (number of cyclones)','position',[190 -58],'Fontsize',16,'FontName','sans serif');

ax1=subplot(ha(18));
cdata0=rgb2mat('/Users/wangxiaoqiao/Downloads/software/MATLAB/rgb_files/amwg_blueyellowred.rgb',40);
    m_proj('Miller', 'long',[159 210], 'lat',[-84 -69]);hold on;
[C,h]=m_contourf(lon,lat,MSL./100,300);hold on
set(h,'LineColor','none');
    m_usercoast('/Users/wangxiaoqiao/Downloads/software/MATLAB/topodata/whole_Antarctic/topodata.mat','linewi',2,'color',[.4 .4 .4]);hold on
        m_grid('linewi',1,'tickdir','in','tickstyle','dd','xtick',160:15:210,...
            'ytick',-84:5:-69,'yticklables',[],'fontsize',14, 'ticklen',0.007,'FontName','sans serif');hold on;
m_text(160,-67.5,str(18),'Fontsize',15,'Color','k','FontWeight','bold',...
    'HorizontalAlignment','left','FontName','sans serif');
m_vec(31,Lon(1:18:end,1:6:end),Lat(1:18:end,1:6:end),u(1:18:end,1:6:end),v(1:18:end,1:6:end),'k','Shaftwidth',0.8,'headlength',4.5,'centered','yes','edgeclip','on');hold on;
m_vec(31,[164],[-82],[5],[0],'r','Shaftwidth',1.5,'headlength',4.5,'centered','yes','edgeclip','on');hold on;
m_text(162,-83,'5 m s^{-1}','Fontsize',14,...
    'FontWeight','bold','Color','r',...
    'FontName','sans serif');

% m_text(155,-57,"Sea level pressure",'Fontsize',16,'Color','k','FontWeight','bold',...
%     'HorizontalAlignment','left','FontName','sans serif');hold on

box on

colormap(ax1,cdata0);
caxis(ax1,[950 980]);
cb0=colorbar;hold on
set(cb0,'YTick',950:10:980,'tickdir','in','ticklen',0.007,'Fontsize',13,'FontName','sans serif','position',[0.89,0.12,0.01,0.18]);hold on
set(get(cb0,'title'),'string','[hPa]','Fontsize',14,'FontName','sans serif');

%%
print(img,'-djpeg','-r300','/Volumes/Expansion/case_2007/0.figure_modified/cyclone_density/trackdenisty_win_2025_0305.jpg');
close all
