Index exceeds the number of array elements. Index must not exceed 12784.

[Smabmn]

When I run the example this error appears

Index exceeds the number of array elements. Index must not exceed 12784.

Error in event_line (line 84)
y1=m90_plot(x1-datenum(data_start,1,1)+1);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

[ZijieZhaoMMHW]

@Smabmn
Dear Smabmn,

May I ask how you run the code? I run the example on my local machine and there is no error reported.

Cheers,
Zijie

[Smabmn]

i used 2024b
I followed this steps
1 add to my path
open an_example 2.mat
Run

[ZijieZhaoMMHW]

@Smabmn
Are you using my data or your own data? did you modify the code in an_example 2.mat?

[Smabmn]

your data.

%% An example about how to apply m_mhw to real-world data
% Here we provide an example about how to apply m_mhw to real-world data.
addpath(genpath('C:\Users\Marine Physics\MATLAB Drive\Heat_Wave'));
%% 1. Loading data

% Load NOAA OI SST V2 data
sst_full=NaN(32,32,datenum(2016,12,31)-datenum(1982,1,1)+1);
for i=1982:2016;
file_here=['sst_' num2str(i)];
load(file_here);
eval(['data_here=sst_' num2str(i) ';'])
sst_full(:,:,(datenum(i,1,1):datenum(i,12,31))-datenum(1982,1,1)+1)=data_here;
end

% This data includes SST in [147-155E, 45-37S] in resolution of 0.25 from
% 1982 to 2016.

load('lon_and_lat');
size(sst_full); %size of data
datenum(2016,12,31)-datenum(1982,1,1)+1 % The temporal length is 35 years.
%% 2. Detecting MHWs and MCSs

% Here we detect marine heatwaves off eastern Tasmania based on the
% traditional definition of MHWs (Hobday et al. 2016). We detected MHWs
% during 1993 to 2016 for climatologies and thresholds in 1982 to 2005.
tic
[MHW,mclim,m90,mhw_ts]=detect(sst_full,datenum(1982,1,1):datenum(2016,12,31),datenum(1982,1,1),datenum(2005,12,31),datenum(1993,1,1),datenum(2016,12,31));; %take about 30 seconds.
toc

% Additionally, we also detect MCSs during 1982 to 2005 based on the same
% climatologies.

[MCS,~,m10,mcs_ts]=detect(sst_full,datenum(1982,1,1):datenum(2016,12,31),datenum(1982,1,1),datenum(2005,12,31),datenum(1993,1,1),datenum(2016,12,31),'Event','MCS','Threshold',0.1);

% Have a look of these two data.
MHW(1:5,:);
MCS(1:5,:);

% You could see that the properties mhw_onset and mhw_end are in a
% strange format. This is due to the fact that they are originally
% constructed in numeric format. We could change it to date format by
% following steps.

datevec(num2str(MHW{1:5,:}),'yyyymmdd') % vector
datestr(datevec(num2str(MHW{1:5,:}),'yyyymmdd')) % string
datenum(num2str(MHW{1:5,:}),'yyyymmdd') % number

%% 3. Visualizing MHW/MCS time series

% Have a look of MHW events in grid (1,2) from Sep 2015 to Apr 2016
figure('pos',[10 10 1000 1000]);
event_line(sst_full,MHW,mclim,m90,[1 2],1982,[2015 9 1],[2016 5 1]);

% Have a look of MCS events in grid (1,2) during 1994
figure('pos',[10 10 1000 1000]);
event_line(sst_full,MCS,mclim,m10,[1 2],1982,[1994 1 1],[1994 12 31],'Event','MCS','Color',[0.5 0.5 1]);

% %% 4. Mean states and trends

% Please note that this section requires the toolbox m_map

% Now we would like to know the mean states and annual trends of MHW
% frequency, i.e. how many MHW events would be detected per year and how it
% changes with time.

[mean_freq,annual_freq,trend_freq,p_freq]=mean_and_trend_new(MHW,mhw_ts,1982,'Metric','Frequency');

% These four outputs separately represent the total mean, annual mean,
% annual trend and associated p value of frequency.

% This function could detect mean states and trends for six different
% variables (Frequency, mean intensity, max intensity, duration and total
% MHW/MCs days).

metric_used={'Frequency','MeanInt','MaxInt','CumInt','Duration','Days'};

for i=1:6;
eval(['[mean_' metric_used{i} ',annual_' metric_used{i} ',trend_' metric_used{i} ',p_' metric_used{i} ']=mean_and_trend_new(MHW,mhw_ts,1982,' '''' 'Metric' '''' ',' 'metric_used{i}' ');'])
end

% plot mean and trend

% It could be detected that, as a global hotspot, the oceanic region off
% eastern Tasmania exhibits significant positive trends of MHW metrics.

figure('pos',[10 10 1500 1500]);
m_proj('mercator','lat',[-45 -37],'lon',[147 155]);

for i=1:6;
subplot(2,6,i);
eval(['mean_here=mean_' metric_used{i} ';']);
eval(['t_here=trend_' metric_used{i} ';']);
m_pcolor(lon_used,lat_used,mean_here');
shading interp
m_coast('patch',[0.7 0.7 0.7]);
m_grid;
colormap(jet);
s=colorbar('location','southoutside');
title(metric_used{i},'fontname','consolas','fontsize',12);

subplot(2,6,i+6);
eval(['mean_here=mean_' metric_used{i} ';']);
eval(['t_here=trend_' metric_used{i} ';']);
m_pcolor(lon_used,lat_used,t_here');
shading interp
m_coast('patch',[0.7 0.7 0.7]);
m_grid;
colormap(jet);
s=colorbar('location','southoutside');
title(['Trend-' metric_used{i}],'fontname','consolas','fontsize',12);

end

%% 5. Applying cluster algoirthm to MHW - A kmeans example.

% We get so many MHWs now.... Could we distinguish them into different
% gropus based on their metrics?

% Change it to matrix;
MHW_m=MHW{:,:};

% Extract mean, max, cumulative intensity and duration.
MHW_m=MHW_m(:,[3 4 5 7]);

[data_for_k,mu,sd]=zscore(MHW_m);

% Determine suitable groups of kmeans cluster.
index_full=[];
cor_full=[];
for i=2:20;
k=kmeans(data_for_k,i,'Distance','cityblock','maxiter',200);
k_full=[];
for j=1:i;
k_full=[k_full;nanmean(data_for_k(k==j,:))];
end

k_cor=k_full(k,:);
k_cor=k_cor(:);

[c,p]=corr([data_for_k(:) k_cor]);

index_full=[index_full;2];
cor_full=[cor_full;c(1,2)];

end

figure('pos',[10 10 1500 1500]);
subplot(1,2,1);
plot(2:20,cor_full,'linewidth',2);
hold on
plot(9*ones(1000,1),linspace(0.6,1,1000),'r--');
xlabel('Number of Groups','fontsize',16,'fontweight','bold');
ylabel('Correlation','fontsize',16,'fontweight','bold');
title('Correlation','fontsize',16);
set(gca,'xtick',[5 9 10 15 20],'fontsize',16);

subplot(1,2,2);
plot(3:20,diff(cor_full),'linewidth',2);
hold on
plot(9*ones(1000,1),linspace(-0.02,0.14,1000),'r--');
xlabel('Number of Groups','fontsize',16,'fontweight','bold');
ylabel('First difference of Correlation','fontsize',16,'fontweight','bold');
title('First Difference of Correlation','fontsize',16);
set(gca,'fontsize',16);

% Use 9 groups.

k=kmeans(data_for_k,9,'Distance','cityblock','maxiter',200);

k_9=[];
prop_9=[];
for i=1:9;
data_here=data_for_k(k==i,:);
data_here=nanmean(data_here);
data_here=data_here.*sd+mu;
k_9=[k_9;data_here];
prop_9=[prop_9;nansum(k==i)./size(data_for_k,1)];
end

loc_x=[1.5 1.5 1.5 2.5 2.5 2.5 3.5 3.5 3.5];
loc_y=[1.5 2.5 3.5 1.5 2.5 3.5 1.5 2.5 3.5];
text_used={['1: ' num2str(round(prop_9(1)*100)) '%'],['2: ' num2str(round(prop_9(2)*100)) '%'],['3: ' num2str(round(prop_9(3)*100)) '%'],...
['4: ' num2str(round(prop_9(4)*100)) '%'],['5: ' num2str(round(prop_9(5)*100)) '%'],['6: ' num2str(round(prop_9(6)*100)) '%'],...
['7: ' num2str(round(prop_9(7)*100)) '%'],['8: ' num2str(round(prop_9(8)*100)) '%'],['9: ' num2str(round(prop_9(9)*100)) '%']};

figure('pos',[10 10 1500 1500]);

h=subplot(2,2,1);
data_here=k_9(:,1);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
pcolor(1:4,1:4,data_here);
set(h,'ydir','reverse');
axis off
colormap(jet);
text(loc_x,loc_y,text_used,'fontsize',16,'horiz','center','fontweight','bold');
colorbar;
title('Durations','fontsize',16,'fontweight','bold');

h=subplot(2,2,2);
data_here=k_9(:,2);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
pcolor(1:4,1:4,data_here);
axis off
set(h,'ydir','reverse');
colormap(jet);
text(loc_x,loc_y,text_used,'fontsize',16,'horiz','center','fontweight','bold');
colorbar
title('MaxInt','fontsize',16,'fontweight','bold');

h=subplot(2,2,3);
data_here=k_9(:,3);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
pcolor(1:4,1:4,data_here);
axis off
set(h,'ydir','reverse');
colormap(jet);
text(loc_x,loc_y,text_used,'fontsize',16,'horiz','center','fontweight','bold');
colorbar
title('MeanInt','fontsize',16,'fontweight','bold');

h=subplot(2,2,4);
data_here=k_9(:,4);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
[x,y]=meshgrid(1:4,1:4);
pcolor(1:4,1:4,data_here);
set(h,'ydir','reverse');
axis off
colormap(jet);
text(loc_x,loc_y,text_used,'fontsize',16,'horiz','center','fontweight','bold');
colorbar
title('CumInt','fontsize',16,'fontweight','bold');

% Their associated SSTA patterns

% Calculate SSTA

time_used=datevec(datenum(1982,1,1):datenum(2016,12,31));
m_d_unique=unique(time_used(:,2:3),'rows');

ssta_full=NaN(size(sst_full));

for i=1:size(m_d_unique);
date_here=m_d_unique(i,:);
index_here=find(time_used(:,2)==date_here(1) & time_used(:,3)==date_here(2));
ssta_full(:,:,index_here)=sst_full(:,:,index_here)-nanmean(sst_full(:,:,index_here),3);
end

sst_1993_2016=ssta_full(:,:,(datenum(1993,1,1):datenum(2016,12,31))-datenum(1982,1,1)+1);
time_used=MHW{:,:};
time_used=time_used(:,1:2);
start_full=datenum(num2str(time_used(:,1)),'yyyymmdd')-datenum(1993,1,1)+1;
end_full=datenum(num2str(time_used(:,2)),'yyyymmdd')-datenum(1993,1,1)+1;

for i=1:9;
start_here=start_full(k==i);
end_here=end_full(k==i);

index_here=[];

for j=1:length(start_here);
    period_here=start_here(j):end_here(j);
    index_here=[index_here;period_here(:)];
end

eval(['sst_' num2str(i) '=nanmean(sst_1993_2016(:,:,index_here),3);'])

end

color_used=hot;
color_used=color_used(end:-1:1,:);
figure('pos',[10 10 1500 1500]);

plot_index=[1 4 7 2 5 8 3 6 9];

for i=1:9;
eval(['plot_here=sst_' num2str(i) ';']);
subplot(3,3,plot_index(i));
eval(['data_here=sst_' num2str(i) ';'])

m_contourf(lon_used,lat_used,data_here',-3:0.1:3,'linestyle','none');


if i~=3;
m_grid('xtick',[],'ytick',[]);
else;
    m_grid('linestyle','none');
end
m_gshhs_h('patch',[0 0 0]);
colormap(color_used);
caxis([0 2]);

title(['Group (' num2str(i) '):' num2str(round(prop_9(i)*100)) '%'],'fontsize',16);

end

hp4=get(subplot(3,3,9),'Position');
s=colorbar('Position', [hp4(1)+hp4(3)+0.025 hp4(2) 0.025 0.85],'fontsize',14);
s.Label.String='^{o}C';

[ZijieZhaoMMHW]

@Smabmn
May I ask what the differences between an_example2.m and my original an_example.m?

Cheers,
Zijie

[Smabmn]

No differences,
but in an_example2.m , { m_proj('mercator','lat',[-45 -37],'lon',[147 155]);} not writing in original an_example.m
and ,i haven't used my data for the code yet.

[ZijieZhaoMMHW]

@Smabmn
Thanks. But the output from your run (shown in the attached image) looks very weird as it doesn't corresponding to any date num during 1982-2016. Must be something unexpected happening. I will try to run the code on the latest MATLAB version to test if it's a version issue but I truly doubt it as MATLAB hardly change its underlying code logic through version updating. I suggest you run the code again following my example step by step to see how it works
Cheers,
Zijie