Hi,I tray RUN this code but i get tis error ( below the code)
5 views (last 30 days)
Show older comments
% 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';
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);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in an_example2 (line 53)
event_line(sst_full,MHW,mclim,m90,[1 2],1982,[2015 9 1],[2016 5 1]);
1 Comment
Walter Roberson
on 6 Oct 2024
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}' ');'])
I'm out. There are very very seldom good reasons to use eval()
Answers (3)
Torsten
on 6 Oct 2024
Moved: Image Analyst
on 6 Oct 2024
Variables and/or functions are missing to run your code.
According to the error message, m90_plot is an array of size 12784, but the number x1-datenum(data_start,1,1)+1 is bigger than 12784. Thus you get an access error.
0 Comments
Image Analyst
on 6 Oct 2024
Edited: Image Analyst
on 6 Oct 2024
@B Your index is bigger than 12784, which is the max index m90_plot can take. Put these lines in and tell us what you see in the command window.
size(m90_plot)
fprintf('data_start = %g\n', data_start);
fprintf('datenum(data_start,1,1) = %g\n', datenum(data_start,1,1));
index = x1-datenum(data_start,1,1)+1
fprintf('x1-datenum(data_start,1,1)+1 = %g\n', index);
m90_plot is an array so the index must be an integer. Is index an integer or a floating point number.
datenum most likely returns a huge number. Look
datenum(2024, 10, 6)
Not sure what you're attempting to do there but it doesn't seem right unless x1 is also a datenum value.
Also, don't use eval. Learn how to use fprintf
Instead of
eval(['mean_here=mean_' metric_used{i} ';']);
eval(['t_here=trend_' metric_used{i} ';']);
do
fprintf('mean_here = mean_%s\n', metric_used{i});
fprintf('t_here = trend_%s\n', metric_used{i});
0 Comments
B
on 8 Oct 2024
1 Comment
Image Analyst
on 8 Oct 2024
Thank who? Who is your "Answer" supposed to be a comment to? I think you did not even see my Answer otherwise you would have answered when I said:
Put these lines in and tell us what you see in the command window.
size(m90_plot)
fprintf('data_start = %g\n', data_start);
fprintf('datenum(data_start,1,1) = %g\n', datenum(data_start,1,1));
index = x1-datenum(data_start,1,1)+1
fprintf('x1-datenum(data_start,1,1)+1 = %g\n', index);
So go back up and read my answer and put your reply in the comment below it.
See Also
Categories
Find more on Matrix Indexing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!