海洋热浪(Marine Heat Waves，MHW)是指在一定海域内发生的海表温度异常偏高的现象。本文针对2023年的中国近海海洋热浪绘制海洋热浪的发生天数，最大强度和平均强度的直观图。
本文选取的中国近海区域为10°～40°N，105°～130°E，主要包括渤海、黄海、东海和南海以及台湾东侧的太平洋的部分海区。
本文选取的数据为美国国家海洋和大气管理局的最佳插值海平面温度（OISST）日分辨率海温数据集（V2）。该数据集使用了改进的高分辨率雷达（AVHRR）红外卫星海温资料以及浮标和船舶等原位观测资料。日平均海温数据的空间分辨率为0.25°×0.25°，时间跨度为1981年6月至2024年5月。由于0.25的分辨率对于我的课题太精细了所以我的数据是已经经过线性插值处理过的，分辨率为1°×1°。
本文采用Hobdy等关于海洋热浪的定义，即海洋热浪是指在一定海域内发生的日海表温度至少连续5 天超过当地季节阈值（即气候基准期内同期日海表温度的第90百分位）的事件，其持续时间可达数月，空间范围可延伸至数千千米。
本文海洋热浪所采用的气候基准期为1993-2022年，计算当天气候态值及阈值。
发生总天数和频次
start_date=1993;
date=2023;
current_date=start_date;
for i=1:30
fn1=["E:\sst\sst.day.mean."+num2str(current_date)+".nc"];
sst1=ncread(fn1,'sst');
sst1=sst1(420:4:520,400:4:520,:);
sst(i,:,:,:)=sst1;
current_date=start_date+1;
end
sst(sst<=0)=NaN;
path0="C:\sst.day.mean.2023.nc";
sst0=ncread(path0,'sst');
sst1=sst0(420:4:520,400:4:520,:);
sst1(sst1<=0)=NaN;
day2=zeros(26,31,365);
for j=1:26
for k=1:31
for i=1:365
daya=sst(:,j,k,i);
day1=sort(daya);
if sst1(j,k,i)>=day1(27)
day2(j,k,i)=1;
else
day2(j,k,i)=0;
end
end
end
end
A1=day2;
for j=1:26
for k=1:31
count=0;
for i = 1:365
if A1(j,k,i)>0 % 如果元素大于0
count=count+1; % 增加计数器
else % 如果元素不大于0
if count<5 % 如果之前有连续5个以上大于0的元素
A1(j,k,i-count+1:i)=0;
end
if count>=5
A1(j,k,i-count+2:i)=0;
end
count = 0; % 重置计数器和起始索引
end
end
end
end
count = sum(A1(:) == 1);
%8176
A0=day2;
for j=1:26
for k=1:31
for i = 1:365
if A0(j,k,i)>0 % 如果元素大于0
count=0;
while i<=365 && A0(j,k,i)>0
count=count+1;
i=i+1;
end
if count<5
A0(j,k,(i-count):(i-1))=0;
else
A0(j,k,(i-count+1):(i-1))=0;
end
end
end
end
end
A=sum(A1,3);
B=sum(A0,3);
B=flipud(B);
B=rot90(rot90(rot90(B)));
lon=104.875:1:129.875;
lat=9.875:1:39.875;
figure
m_proj('mercator','lon',[105 130],'lat',[10 40]);
[X1,Y1]=meshgrid(105:1:130,10:1:40);
m_contourf(X1,Y1,B)
m_coast('patch',[.6 .6 .6],'edgecolor','k');
m_grid('linestyle','none','xtick',105:2.5:130,'ytick',10:5:40,'fontsize',8,'tickdir','out');
colormap(m_colmap('jet'));
c = colorbar;
c.Title.String='发生频次/次数';
set(c,'Ticks',0:1:15);
set(gcf,'color','w','position',[100 100 1000 500]);
title(['发生频次'],'fontsize',12,'color','b');
A=flipud(A);
A=rot90(rot90(rot90(A)));
lon=104.875:1:129.875;
lat=9.875:1:39.875;
figure
m_proj('mercator','lon',[105 130],'lat',[10 40]);
[X1,Y1]=meshgrid(105:1:130,10:1:40);
m_contourf(X1,Y1,A,10)
m_coast('patch',[.6 .6 .6],'edgecolor','k');
m_grid('linestyle','none','xtick',105:2.5:130,'ytick',10:5:40,'fontsize',8,'tickdir','out');
colormap(m_colmap('jet'));
c = colorbar;
c.Title.String='发生天数/天';
set(c,'Ticks',0:10:250)
set(gcf,'color','w','position',[100 100 1000 500]);
title(['发生天数'],'fontsize',12,'color','b');
得到的图像为
四季的发生天数(季节划分是3-5月为春季，6-8月为夏季，9-11月为秋季，12月至第二年2月为冬季)
start_date=1993;
date=2023;
current_date=start_date;
%10-40N,105-130E
for i=1:30
fn1=["E:\大三下\卫星海洋学\课设\数据\sst\sst.day.mean."+num2str(current_date)+".nc"];
sst1=ncread(fn1,'sst');
sst1=sst1(420:4:520,400:4:520,:);
sst(i,:,:,:)=sst1;
current_date=start_date+1;
end
sst(sst<=0)=NaN;
path0="E:\sst\sst.day.mean.2023.nc";
sst0=ncread(path0,'sst');
sst1=sst0(420:4:520,400:4:520,:);
sst1(sst1<=0)=NaN;
path2="E:\sst\sst.day.mean.2024.nc";
sst2=ncread(path2,'sst');
sst2=sst2(420:4:520,400:4:520,:);
sst2(sst2<=0)=NaN;
day2=zeros(26,31,92);
for j=1:26
for k=1:31
for i=60:151
daya=sst(:,j,k,i);
day1=sort(daya);
if sst1(j,k,i)>=day1(27)
day2(j,k,i)=1;
else
day2(j,k,i)=0;
end
end
end
end
A1=day2;
for j=1:26
for k=1:31
count=0;
for i = 1:92
if A1(j,k,i)>0 % 如果元素大于0
count=count+1; % 增加计数器
else % 如果元素不大于0
if count<5 % 如果之前有连续5个以上大于0的元素
A1(j,k,i-count+1:i)=0;
end
if count>=5
A1(j,k,i-count+2:i)=0;
end
count = 0; % 重置计数器和起始索引
end
end
end
end
count = sum(A1(:) == 1);
A1=sum(A1,3);
A1=flipud(A1);
A1=rot90(rot90(rot90(A1)));
lon=104.875:1:129.875;
lat=9.875:1:39.875;
figure
subplot(2,2,1)
m_proj('mercator','lon',[105 130],'lat',[10 40]);
[X1,Y1]=meshgrid(105:1:130,10:1:40);
m_contourf(X1,Y1,A1,10)
m_coast('patch',[.6 .6 .6],'edgecolor','k');
m_grid('linestyle','none','xtick',105:2.5:130,'ytick',10:5:40,'fontsize',8,'tickdir','out');
colormap(m_colmap('jet'));
c = colorbar;
c.Title.String='发生天数/天';
set(c,'Ticks',0:10:250)
set(gcf,'color','w','position',[100 100 1000 500]);
title(['春季发生天数'],'fontsize',12,'color','b');
day2=zeros(26,31,92);
for j=1:26
for k=1:31
for i=152:243
daya=sst(:,j,k,i);
day1=sort(daya);
if sst1(j,k,i)>=day1(27)
day2(j,k,i)=1;
else
day2(j,k,i)=0;
end
end
end
end
A2=day2;
for j=1:26
for k=1:31
count=0;
for i = 1:92
if A2(j,k,i)>0 % 如果元素大于0
count=count+1; % 增加计数器
else % 如果元素不大于0
if count<5 % 如果之前有连续5个以上大于0的元素
A2(j,k,i-count+1:i)=0;
end
if count>=5
A2(j,k,i-count+2:i)=0;
end
count = 0; % 重置计数器和起始索引
end
end
end
end
count = sum(A2(:) == 1);
A2=sum(A2,3);
A2=flipud(A2);
A2=rot90(rot90(rot90(A2)));
lon=104.875:1:129.875;
lat=9.875:1:39.875;
subplot(2,2,2)
m_proj('mercator','lon',[105 130],'lat',[10 40]);
[X1,Y1]=meshgrid(105:1:130,10:1:40);
m_contourf(X1,Y1,A2,10)
m_coast('patch',[.6 .6 .6],'edgecolor','k');
m_grid('linestyle','none','xtick',105:2.5:130,'ytick',1:5:40,'fontsize',8,'tickdir','out');
colormap(m_colmap('jet'));
c = colorbar;
c.Title.String='发生天数/天';
set(c,'Ticks',0:10:250)
set(gcf,'color','w','position',[100 100 1000 500]);
title(['夏季发生天数'],'fontsize',12,'color','b');
day2=zeros(26,31,91);
for j=1:26
for k=1:31
for i=244:344
daya=sst(:,j,k,i);
day1=sort(daya);
if sst1(j,k,i)>=day1(27)
day2(j,k,i)=1;
else
day2(j,k,i)=0;
end
end
end
end
A3=day2;
for j=1:26
for k=1:31
count=0;
for i = 1:91
if A3(j,k,i)>0 % 如果元素大于0
count=count+1; % 增加计数器
else % 如果元素不大于0
if count<5 % 如果之前有连续5个以上大于0的元素
A3(j,k,i-count+1:i)=0;
end
if count>=5
A3(j,k,i-count+2:i)=0;
end
count = 0; % 重置计数器和起始索引
end
end
end
end
count = sum(A3(:) == 1);
A3=sum(A3,3);
A3=flipud(A3);
A3=rot90(rot90(rot90(A3)));
lon=104.875:1:129.875;
lat=9.875:1:39.875;
subplot(2,2,3)
m_proj('mercator','lon',[105 130],'lat',[10 40]);
[X1,Y1]=meshgrid(105:1:130,10:1:40);
m_contourf(X1,Y1,A3,10)
m_coast('patch',[.6 .6 .6],'edgecolor','k');
m_grid('linestyle','none','xtick',105:2.5:130,'ytick',10:5:40,'fontsize',8,'tickdir','out');
colormap(m_colmap('jet'));
c = colorbar;
c.Title.String='发生天数/天';
set(c,'Ticks',0:10:250)
set(gcf,'color','w','position',[100 100 1000 500]);
title(['秋季发生天数'],'fontsize',12,'color','b');
day2=zeros(26,31,32);
day3=zeros(26,31,60);
for j=1:26
for k=1:31
for i=334:365
daya=sst(:,j,k,i);
day1=sort(daya);
if sst1(j,k,i)>=day1(27)
day2(j,k,i)=1;
else
day2(j,k,i)=0;
end
end
end
end
for j=1:26
for k=1:31
for i=1:60
daya=sst(:,j,k,i);
day1=sort(daya);
if sst2(j,k,i)>=day1(27)
day3(j,k,i)=1;
else
day3(j,k,i)=0;
end
end
end
end
A4=cat(3,day2,day3);
for j=1:26
for k=1:31
count=0;
for i = 1:92
if A4(j,k,i)>0 % 如果元素大于0
count=count+1; % 增加计数器
else % 如果元素不大于0
if count<5 % 如果之前有连续5个以上大于0的元素
A4(j,k,i-count+1:i)=0;
end
if count>=5
A4(j,k,i-count+2:i)=0;
end
count = 0; % 重置计数器和起始索引
end
end
end
end
count = sum(A4(:) == 1);
A4=sum(A4,3);
A4=flipud(A4);
A4=rot90(rot90(rot90(A4)));
lon=104.875:1:129.875;
lat=5.875:1:39.875;
subplot(2,2,4)
m_proj('mercator','lon',[105 130],'lat',[10 40]);
[X1,Y1]=meshgrid(105:1:130,10:1:40);
m_contourf(X1,Y1,A4,10)
m_coast('patch',[.6 .6 .6],'edgecolor','k');
m_grid('linestyle','none','xtick',105:2.5:130,'ytick',10:5:40,'fontsize',8,'tickdir','out');
colormap(m_colmap('jet'));
c = colorbar;
c.Title.String='发生天数/天';
set(c,'Ticks',0:10:250)
set(gcf,'color','w','position',[100 100 1000 500]);
title(['冬季发生天数'],'fontsize',12,'color','b');
得到的图像为 
平均强度和最大强度
clear;clc;
start_date=1993;
date=2023;
current_date=start_date;
for i=1:30
fn1=["E:\sst\sst.day.mean."+num2str(current_date)+".nc"];
sst1=ncread(fn1,'sst');
sst1=sst1(420:4:520,400:4:520,:);
sst(i,:,:,:)=sst1;
current_date=start_date+1;
end
sst=mean(sst,1);
sst=reshape(sst,26,31,365);
sst(sst<=0)=NaN;
path0="C:\sst\sst.day.mean.2023.nc";
sst0=ncread(path0,'sst');
sst1=sst0(420:4:520,420:4:480,:);
sst1(sst1<=0)=NaN;
C=sst1-sst;
idx=find(A1==1); % 找到值为 1 的元素的索引
[row,col,page]=ind2sub(size(A1),idx); % 将索引转换为行列坐标
a=zeros(26,31,365);
for i=1:length(row)
a(row(i),col(i),page(i))=C(row(i),col(i),page(i));
end
a0=sum(a,3)/365;
a0=flipud(a0);
a0=rot90(rot90(rot90(a0)));
lon=104.875:1:129.875;
lat=9.875:1:39.875;
figure
m_proj('mercator','lon',[105 130],'lat',[10 30]);
[X1,Y1]=meshgrid(105:1:130,10:1:40);
m_contourf(X1,Y1,a0)
m_coast('patch',[.6 .6 .6],'edgecolor','k');
m_grid('linestyle','none','xtick',105:2.5:130,'ytick',10:5:40,'fontsize',8,'tickdir','out');
colormap(m_colmap('jet'));
c = colorbar;
c.Title.String='平均强度/℃';
set(gcf,'color','w','position',[100 100 1000 500]);
title(['平均强度'],'fontsize',12,'color','b');
a1=zeros(26,16);
for i=1:26
for j=1:16
a1(i,j)=max(a(i,j,:));
end
end
a1=flipud(a1);
a1=rot90(rot90(rot90(a1)));
lon=104.875:1:129.875;
lat=14.875:1:29.875;
figure
m_proj('mercator','lon',[105 130],'lat',[15 30]);
[X1,Y1]=meshgrid(105:1:130,15:1:30);
m_contourf(X1,Y1,a1)
m_coast('patch',[.6 .6 .6],'edgecolor','k');
m_grid('linestyle','none','xtick',105:2.5:130,'ytick',15:2.5:30,'fontsize',8,'tickdir','out');
colormap(m_colmap('jet'));
c = colorbar;
c.Title.String='最大强度/℃';
set(gcf,'color','w','position',[100 100 1000 500]);
title(['最大强度'],'fontsize',12,'color','b');
得到的图像为 
四季的强度为
clear
start_date=1993;
date=2023;
current_date=start_date;
for i=1:30
fn1=["E:\sst\sst.day.mean."+num2str(current_date)+".nc"];
sst1=ncread(fn1,'sst');
sst1=sst1(420:4:520,400:4:520,:);
sst(i,:,:,:)=sst1;
current_date=start_date+1;
end
sst=mean(sst,1);
sst=reshape(sst,26,31,365);
sst(sst<=0)=NaN;
path0="C:\sst\sst.day.mean.2023.nc";
sst0=ncread(path0,'sst');
sst1=sst0(420:4:520,400:4:520,:);
sst1(sst1<=0)=NaN;
path2="E:\sst\sst.day.mean.2024.nc";
sst2=ncread(path2,'sst');
sst2=sst2(420:4:520,400:4:520,:);
sst2(sst2<=0)=NaN;
day2=zeros(26,31,365);
for j=1:26
for k=1:31
for i=1:30
daya=sst(:,j,k,i);
day1=sort(daya);
if sst1(j,k,i)>=day1(27)
day2(j,k,i)=1;
else
day2(j,k,i)=0;
end
end
end
end
A1=day2;
for j=1:26
for k=1:31
count=0;
for i = 1:92
if A1(j,k,i)>0 % 如果元素大于0
count=count+1; % 增加计数器
else % 如果元素不大于0
if count<5 % 如果之前有连续5个以上大于0的元素
A1(j,k,i-count+1:i)=0;
end
if count>=5
A1(j,k,i-count+2:i)=0;
end
count = 0; % 重置计数器和起始索引
end
end
end
end
count = sum(A1(:) == 1);
A1=sum(A1,3);
A1=flipud(A1);
A1=rot90(rot90(rot90(A1)));
%3-5月为春季，6-8月为夏季，9-11月为秋季，12月至第二年2月为冬季
sst1_1=sst1(:,:,60:151);
sst1_2=sst1(:,:,152:243);
sst1_3=sst1(:,:,244:334);
sst1_41=sst1(:,:,334:365);
sst1_42=sst2(:,:,1:60);
sst1_4=cat(3,sst1_41,sst1_42);
sst_1=sst(:,:,1:92);
sst_2=sst(:,:,93:184);
sst_3=sst(:,:,184:274);
sst_4=sst(:,:,274:365);
c1=sst1_1-sst_1;
c2=sst1_2-sst_2;
c3=sst1_3-sst_3;
c4=sst1_4-sst_4;
c1(c1<=0)=0;
c2(c2<=0)=0;
c3(c3<=0)=0;
c4(c4<=0)=0;
A1_1=A1(:,:,1:92);
A1_2=A1(:,:,93:184);
A1_3=A1(:,:,184:274);
A1_4=A1(:,:,274:365);
idx1=find(A1_1==1); % 找到值为 1 的元素的索引
[row1,col1,page1]=ind2sub(size(A1_1),idx1); % 将索引转换为行列坐标
a0=zeros(26,16,91);
for i=1:length(row1)
a0(row1(i),col1(i),page1(i))=c1(row1(i),col1(i),page1(i));
end
a1=sum(a0,3)/91;
a1=flipud(a1);
a1=rot90(rot90(rot90(a1)));
idx2=find(A1_2==1); % 找到值为 1 的元素的索引
[row2,col2,page2]=ind2sub(size(A1_2),idx2); % 将索引转换为行列坐标
a0=zeros(26,31,91);
for i=1:length(row2)
a0(row2(i),col2(i),page2(i))=c2(row2(i),col2(i),page2(i));
end
a2=sum(a0,3)/91;
a2=flipud(a2);
a2=rot90(rot90(rot90(a2)));
idx3=find(A1_3==1); % 找到值为 1 的元素的索引
[row3,col3,page3]=ind2sub(size(A1_3),idx3); % 将索引转换为行列坐标
a0=zeros(26,16,91);
for i=1:length(row3)
a0(row3(i),col3(i),page3(i))=c3(row3(i),col3(i),page3(i));
end
a3=sum(a0,3)/91;
a3=flipud(a3);
a3=rot90(rot90(rot90(a3)));
idx4=find(A1_4==1); % 找到值为 1 的元素的索引
[row4,col4,page4]=ind2sub(size(A1_4),idx4); % 将索引转换为行列坐标
a0=zeros(26,31,92);
for i=1:length(row4)
a0(row4(i),col4(i),page4(i))=c4(row4(i),col4(i),page4(i));
end
a4=sum(a0,3)/92;
a4=flipud(a4);
a4=rot90(rot90(rot90(a4)));
lon=104.875:1:129.875;
lat=14.875:1:29.875;
figure
subplot(2,2,1)
m_proj('mercator','lon',[105 130],'lat',[15 30]);
[X1,Y1]=meshgrid(105:1:130,15:1:30);
m_contourf(X1,Y1,a1)
m_coast('patch',[.6 .6 .6],'edgecolor','k');
m_grid('linestyle','none','xtick',105:2.5:130,'ytick',15:2.5:30,'fontsize',8,'tickdir','out');
colormap(m_colmap('jet'));
c = colorbar;
c.Title.String='春季平均强度/℃';
set(gcf,'color','w','position',[100 100 1000 500]);
title(['春季平均强度'],'fontsize',12,'color','b');
subplot(2,2,2)
m_proj('mercator','lon',[105 130],'lat',[15 30]);
[X1,Y1]=meshgrid(105:1:130,15:1:30);
m_contourf(X1,Y1,a2)
m_coast('patch',[.6 .6 .6],'edgecolor','k');
m_grid('linestyle','none','xtick',105:2.5:130,'ytick',15:2.5:30,'fontsize',8,'tickdir','out');
colormap(m_colmap('jet'));
c = colorbar;
c.Title.String='夏季平均强度/℃';
set(gcf,'color','w','position',[100 100 1000 500]);
title(['夏季平均强度'],'fontsize',12,'color','b');
subplot(2,2,3)
m_proj('mercator','lon',[105 130],'lat',[10 40]);
[X1,Y1]=meshgrid(105:1:130,15:1:30);
m_contourf(X1,Y1,a3)
m_coast('patch',[.6 .6 .6],'edgecolor','k');
m_grid('linestyle','none','xtick',105:2.5:130,'ytick',10:5:40,'fontsize',8,'tickdir','out');
colormap(m_colmap('jet'));
c = colorbar;
c.Title.String='秋季平均强度/℃';
set(gcf,'color','w','position',[100 100 1000 500]);
title(['秋季平均强度'],'fontsize',12,'color','b');
subplot(2,2,4)
m_proj('mercator','lon',[105 130],'lat',[15 30]);
[X1,Y1]=meshgrid(105:1:130,15:1:30);
m_contourf(X1,Y1,a4)
m_coast('patch',[.6 .6 .6],'edgecolor','k');
m_grid('linestyle','none','xtick',105:2.5:130,'ytick',15:2.5:30,'fontsize',8,'tickdir','out');
colormap(m_colmap('jet'));
c = colorbar;
c.Title.String='冬季平均强度/℃';
set(gcf,'color','w','position',[100 100 1000 500]);
title(['冬季平均强度'],'fontsize',12,'color','b');
得到的图像为
本文读取nc文件的步长为4是因为在插值的时候把值定义为0，0，0，a，a为我插值的值。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
