diff --git a/detect.m b/detect.m index ee20dc2..f915e9b 100644 --- a/detect.m +++ b/detect.m @@ -50,6 +50,9 @@ % mhw_end - A numeric value in format of datennum(yyyy,mm,dd) indicating the end year for the period across % which MHW/MCS events are detected. % +% mhw_peak - A numeric value in format of datennum(yyyy,mm,dd) indicating the date of maximum intensity for the period across +% which MHW/MCS events are detected. +% % 'Event' - Default is 'MHW'. % - 'MHW' - detecting MHW events. % - 'MCS' - detecting MCS events. @@ -82,7 +85,8 @@ % corresponds to a particular event and each column corresponds to a % metric. Specified metrics are: % - mhw_onset - onset date of each event. -% - mhw_end - end date of each event. +% - mhw_end - end date of each event. +% - mhw_peak - peak date of each event. % - mhw_dur - duration of each event. % - int_max - maximum intensity of each event. % - int_mean - mean intensity of each event. @@ -99,38 +103,47 @@ % (m-by-n-by-(datenum(MHW_end,1,1)-datenum(MHW_start)+1)) containing % spatial intensity of MHW/MCS in each day. +tb = 0;% temporarily for linux as toolbox not available + +if license('test','Statistics_Toolbox')==0 || tb == 0 + vEvent = 'MHW'; + vThreshold = 0.9; + vWindowHalfWidth = 5; + vsmoothPercentileWidth = 31; + vminDuration = 5; + vmaxGap = 2; + ClimTemp = temp; + ClimTime = time; + vpercentile = 'python';%Because matlab version of percentile needs toolbox + + paramNames = {'vEvent','vThreshold','vWindowHalfWidth','vsmoothPercentileWidth','vminDuration',... + 'vmaxGap','ClimTemp','ClimTime','vpercentile'}; + defaults = {'MHW',0.9,5,31,5,2,temp,time,'python'}; + + varargin = reshape(varargin,2,length(varargin)/2); + + for i = 1:length(defaults) + if any(ismember(varargin(1,:),paramNames{i})) + feval(@()assignin('caller',paramNames{i},varargin{2,ismember(varargin(1,:),paramNames{i})})) + end + end + +else + paramNames = {'vEvent','vThreshold','vWindowHalfWidth','vsmoothPercentileWidth','vminDuration',... + 'vmaxGap','vClimTemp','vClimTime','vpercentile'}; + defaults = {'MHW',0.9,5,31,5,2,temp,time,'matlab'}; + -% vEvent = 'MHW'; -% vThreshold = 0.9; -% vWindowHalfWidth = 5; -% vsmoothPercentileWidth = 31; -% vminDuration = 5; -% vmaxGap = 2; -% ClimTemp = temp; -% ClimTime = time; -% -paramNames = {'Event','Threshold','WindowHalfWidth','smoothPercentileWidth','minDuration',... - 'maxGap','ClimTemp','ClimTime','percentile'}; -defaults = {'MHW',0.9,5,31,5,2,temp,time,'matlab'}; -% -% varargin = reshape(varargin,2,length(varargin)/2); -% -% for i = 1:length(defaults) -% if any(ismember(varargin(1,:),paramNames{i})) -% feval(@()assignin('caller',paramNames{i},varargin{2,ismember(varargin(1,:),paramNames{i})})) -% end -% end - - -[vEvent, vThreshold,vWindowHalfWidth,vsmoothPercentileWidth,vminDuration,vmaxGap,ClimTemp,ClimTime,vpercentile]... - = internal.stats.parseArgs(paramNames, defaults, varargin{:}); - -EventNames = {'MHW','MCS'}; -vEvent = internal.stats.getParamVal(vEvent,EventNames,... - '''Event'''); -PercentileNames={'matlab','python'}; -vpercentile=internal.stats.getParamVal(vpercentile,PercentileNames,... - '''matlab'''); + [vEvent, vThreshold,vWindowHalfWidth,vsmoothPercentileWidth,vminDuration,vmaxGap,ClimTemp,ClimTime,vpercentile]... + = internal.stats.parseArgs(paramNames, defaults, varargin{:}); + + EventNames = {'MHW','MCS'}; + vEvent = internal.stats.getParamVal(vEvent,EventNames,... + '''Event'''); + PercentileNames={'matlab','python'}; + vpercentile=internal.stats.getParamVal(vpercentile,PercentileNames,... + '''matlab'''); +end %% "What if cli_start-window or cli_end+window exceeds the time range of data" @@ -138,7 +151,7 @@ after_date=cli_end+vWindowHalfWidth-time(end); temp_clim=ClimTemp(:,:,ClimTime>=cli_start-vWindowHalfWidth & ClimTime<=cli_end+vWindowHalfWidth); -if ahead_date>0 && after_date>0 +if ahead_date> 0 && after_date>0 temp_clim=cat(3,NaN(size(temp_clim,1),size(temp_clim,2),ahead_date), ... temp_clim,NaN(size(temp_clim,1),size(temp_clim,2),after_date)); elseif ahead_date>0 && after_date<=0 @@ -176,7 +189,9 @@ data_thre=num2cell(temp_clim(:,:,any(ind_fake'>=(ind_fake(fake_doy == i)-vWindowHalfWidth) & ind_fake' <= (ind_fake(fake_doy ==i)+vWindowHalfWidth),2)),3); switch vpercentile case 'matlab' - m90(:,:,i) = quantile(temp_clim(:,:,any(ind_fake'>=(ind_fake(fake_doy == i)-vWindowHalfWidth) & ind_fake' <= (ind_fake(fake_doy ==i)+vWindowHalfWidth),2)),vThreshold,3); + % find out that the quantile function has an issue with singleton dimensions.. + %m90(:,:,i) = quantile(temp_clim(:,:,any(ind_fake'>=(ind_fake(fake_doy == i)-vWindowHalfWidth) & ind_fake' <= (ind_fake(fake_doy ==i)+vWindowHalfWidth),2)),vThreshold,3); + m90(:,:,i) = quantile(squeeze(permute(temp_clim(:,:,any(ind_fake'>=(ind_fake(fake_doy == i)-vWindowHalfWidth) & ind_fake' <= (ind_fake(fake_doy ==i)+vWindowHalfWidth),2)),[3 1 2])),vThreshold);% for singleton dimensions mclim(:,:,i) = mean(temp_clim(:,:,any(ind_fake'>=(ind_fake(fake_doy == i)-vWindowHalfWidth) & ind_fake' <= (ind_fake(fake_doy ==i)+vWindowHalfWidth),2)),3,'omitnan'); case 'python' @@ -207,7 +222,7 @@ date_mhw(:,1)=2000; indextocal = day(datetime(date_mhw),'dayofyear'); -ts=str2double(string(datestr(mhw_start:mhw_end,'YYYYmmdd'))); +ts= mhw_start:mhw_end; mhw_ts=zeros(x_size,y_size,length(ts)); @@ -271,7 +286,8 @@ end mhwstart=NaN(size(potential_event,1),1); - mhwend=NaN(size(potential_event,1),1); + mhwend=NaN(size(potential_event,1),1); + mhwpeak=NaN(size(potential_event,1),1); mduration=NaN(size(potential_event,1),1); mhwint_max=NaN(size(potential_event,1),1); mhwint_mean=NaN(size(potential_event,1),1); @@ -287,8 +303,8 @@ manom=mrow-mcl; mhw_ts(i,j,event_here(1):event_here(2))=manom; - [maxanom,~]=nanmax(squeeze(manom)); - + [maxanom,ind]= max(squeeze(manom),[],'omitnan'); % changed so it doesn't require Stats toobox + mhwpeak(le) = ts(event_here(1))+ind-1; mhwint_max(le)=... maxanom; mhwint_mean(le)=... @@ -301,7 +317,7 @@ mhwend(le)=endtime; mduration(le)=event_here(2)-event_here(1)+1; end - MHW=[MHW;[mhwstart mhwend mduration mhwint_max mhwint_mean mhwint_var mhwint_cum repmat(i,size(mhwstart,1),1) repmat(j,size(mhwstart,1),1)]]; + MHW=[MHW;[mhwstart mhwend mhwpeak mduration mhwint_max mhwint_mean mhwint_var mhwint_cum repmat(i,size(mhwstart,1),1) repmat(j,size(mhwstart,1),1)]]; end end end @@ -365,6 +381,7 @@ mhwstart=NaN(size(potential_event,1),1); mhwend=NaN(size(potential_event,1),1); + mhwpeak=NaN(size(potential_event,1),1); mduration=NaN(size(potential_event,1),1); mhwint_max=NaN(size(potential_event,1),1); mhwint_mean=NaN(size(potential_event,1),1); @@ -380,8 +397,8 @@ manom=mrow-mcl; mhw_ts(i,j,event_here(1):event_here(2))=manom; - [maxanom,~]=nanmin(squeeze(manom)); - + [maxanom,ind]=min(squeeze(manom),[],'omitnan'); + mhwpeak(le) = ts(event_here(1))+ind-1; mhwint_max(le)=... maxanom; mhwint_mean(le)=... @@ -395,7 +412,7 @@ mduration(le)=event_here(2)-event_here(1)+1; end - MHW=[MHW;[mhwstart mhwend mduration mhwint_max mhwint_mean mhwint_var mhwint_cum repmat(i,size(mhwstart,1),1) repmat(j,size(mhwstart,1),1)]]; + MHW=[MHW;[mhwstart mhwend mhwpeak mduration mhwint_max mhwint_mean mhwint_var mhwint_cum repmat(i,size(mhwstart,1),1) repmat(j,size(mhwstart,1),1)]]; end end @@ -409,14 +426,21 @@ end end -MHW=table(MHW(:,1),MHW(:,2),MHW(:,3),MHW(:,4),MHW(:,5),MHW(:,6),MHW(:,7),MHW(:,8),MHW(:,9),... - 'variablenames',{'mhw_onset','mhw_end','mhw_dur','int_max','int_mean','int_var','int_cum','xloc','yloc'}); - +if ~isempty(MHW) + MHW=table(MHW(:,1),MHW(:,2),MHW(:,3),MHW(:,4),MHW(:,5),MHW(:,6),MHW(:,7),MHW(:,8),MHW(:,9),MHW(:,10),... + 'variablenames',{'mhw_onset','mhw_end','mhw_peak','mhw_dur','int_max','int_mean','int_var','int_cum','xloc','yloc'}); +else + disp('No MHWs detected'); +end function p=percentile(data,thre) -if nansum(isnan(data))~=length(data) - data=data(~isnan(data)); +if sum(isnan(squeeze(data)),'omitnan')~=length(data)% added squeeze to work with singleton dimensions + data=data(~isnan(data));% getting rid of NaNs will overestimate the threshold suggest changing it to the median + % Also, if for a reason there is a grid point in the dataset with only + % NaNs, this will create an error pos=1+(length(data)-1)*thre; - p=interp1((1:length(data))',sort(data(:)),pos); + p=interp1((1:length(data))',sort(data(:)),pos); %linear interpolation else p=nan; end + + diff --git a/event_line.m b/event_line.m index 2786e0b..0f35423 100644 --- a/event_line.m +++ b/event_line.m @@ -43,23 +43,40 @@ function event_line(temp,MHW,mclim,m90,loc,data_start,date_start,date_end,vararg % % 'Fontsize' - Default is 16. The size of font in plot. -paramNames = {'Event','Color','Alpha','LineWidth','FontSize'}; -defaults = {'MHW',[1 0.5 0.5],0.5,2,16}; + vEvent = 'MHW'; + vColor = [1 0.5 0.5]; + vAlpha = 0.5; + vLineWidth = 2; + vFontSize = 16; + + paramNames = {'vEvent','vColor','vAlpha','vLineWidth','vFontSize'}; + defaults = {'MHW',[1 0.5 0.5],0.5,2,16}; + + varargin = reshape(varargin,2,length(varargin)/2); -[vEvent, vColor,vAlpha,vLineWidth,vFontSize]... - = internal.stats.parseArgs(paramNames, defaults, varargin{:}); + for i = 1:length(defaults) + if any(ismember(varargin(1,:),paramNames{i})) + feval(@()assignin('caller',paramNames{i},varargin{2,ismember(varargin(1,:),paramNames{i})})) + end + end -EventNames = {'MHW','MCS'}; -vEvent = internal.stats.getParamVal(vEvent,EventNames,... - '''Event'''); +% paramNames = {'Event','Color','Alpha','LineWidth','FontSize'}; +% defaults = {'MHW',[1 0.5 0.5],0.5,2,16}; +% +% [vEvent, vColor,vAlpha,vLineWidth,vFontSize]... +% = internal.stats.parseArgs(paramNames, defaults, varargin{:}); +% +% EventNames = {'MHW','MCS'}; +% vEvent = internal.stats.getParamVal(vEvent,EventNames,... +% '''Event'''); temp_here=squeeze(temp(loc(1),loc(2),:)); MHW=MHW{:,:}; -MHW=MHW(MHW(:,8)==loc(1) & MHW(:,9)==loc(2),:); +MHW=MHW(MHW(:,9)==loc(1) & MHW(:,10)==loc(2),:); period_plot=datenum(data_start,1,1):1:(datenum(data_start,1,1)+size(temp,3)-1); -period_mhw=[datenum(num2str(MHW(:,1)),'yyyymmdd') datenum(num2str(MHW(:,2)),'yyyymmdd')]; +period_mhw=[MHW(:,1) MHW(:,2)]; period_plot_v=datevec(period_plot); period_unique=datevec(datenum(2016,1,1):datenum(2016,12,31)); @@ -79,7 +96,7 @@ function event_line(temp,MHW,mclim,m90,loc,data_start,date_start,date_end,vararg switch vEvent case 'MHW' - for i=1:size(MHW,1) + for i=1:size(MHW,1); MHW_here=period_mhw(i,:); x1=(MHW_here(1):MHW_here(2))'; y1=m90_plot(x1-datenum(data_start,1,1)+1); @@ -96,7 +113,7 @@ function event_line(temp,MHW,mclim,m90,loc,data_start,date_start,date_end,vararg case 'MCS' - for i=1:size(MHW,1) + for i=1:size(MHW,1); MHW_here=period_mhw(i,:); x1=(MHW_here(1):MHW_here(2))'; y1=temp_here(x1-datenum(data_start,1,1)+1); diff --git a/mean_and_trend.m b/mean_and_trend.m index 5225969..8e11de9 100644 --- a/mean_and_trend.m +++ b/mean_and_trend.m @@ -1,4 +1,4 @@ -function [mean_metric,annual_metric,trend_metric,p_metric]=mean_and_trend(MHW,mhw_ts,data_start,varargin) +function [mean_metric,annual_metric,trend_metric,p_metric]=mean_and_trend(MHW,mhw_ts,data_start,varargin); %mean_and_trend - calculating mean states and trends of event metrics % Syntax % [mean_metric,annual_metric,trend_metric,p_metric]=mean_and_trend(MHW,mhw_ts,data_start); @@ -42,27 +42,39 @@ % % p_metric - A 2D numeric matrix (m-by-n) containing p value of % TREND_METRIC. + vMetric = 'Frequency'; + + + paramNames = {'vMetric'}; + defaults = {'Frequency'}; + varargin = reshape(varargin,2,length(varargin)/2); -paramNames = {'Metric'}; -defaults = {'Frequency'}; + for i = 1:length(defaults) + if any(ismember(varargin(1,:),paramNames{i})) + feval(@()assignin('caller',paramNames{i},varargin{2,ismember(varargin(1,:),paramNames{i})})) + end + end -[vMetric]... - = internal.stats.parseArgs(paramNames, defaults, varargin{:}); +% paramNames = {'Metric'}; +% defaults = {'Frequency'}; +% +% [vMetric]... +% = internal.stats.parseArgs(paramNames, defaults, varargin{:}); +% +% MetricNames = {'Frequency','Duration','MaxInt','MeanInt','CumInt','Days'}; +% vMetric = internal.stats.getParamVal(vMetric,MetricNames,... +% '''Metric'''); -MetricNames = {'Frequency','Duration','MaxInt','MeanInt','CumInt','Days'}; -vMetric = internal.stats.getParamVal(vMetric,MetricNames,... - '''Metric'''); - -[x,y,~]=size(mhw_ts); +[x,y,t]=size(mhw_ts); switch vMetric case 'Duration' MHW=MHW{:,:}; - full_mhw_start=datevec(num2str(MHW(:,1)),'yyyymmdd'); - full_mhw_end=datevec(num2str(MHW(:,2)),'yyyymmdd'); - metric=MHW(:,3); + full_mhw_start=datevec(MHW(:,1)); + full_mhw_end=datevec(MHW(:,2)); + metric=MHW(:,4); period_used=datenum(data_start,1,1):(datenum(data_start,1,1)+size(mhw_ts,3)-1); period_used=datevec(period_used); @@ -74,28 +86,28 @@ p_metric=NaN(x,y); annual_metric=NaN(x,y,length(years_mhw)); - loc_full=unique(MHW(:,8:9),'rows'); + loc_full=unique(MHW(:,9:10),'rows'); - for m=1:size(loc_full) + for m=1:size(loc_full); loc_here=loc_full(m,:); - MHW_here=MHW(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2),:); - metric_here=metric(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2)); + MHW_here=MHW(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2),:); + metric_here=metric(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2)); mean_metric(loc_here(1),loc_here(2))=nanmean(metric_here); - for i=1:length(years_mhw) + for i=1:length(years_mhw); year_here=years_mhw(i); - judge_1=(datenum(num2str(MHW_here(:,1)),'yyyymmdd') >= datenum(year_here,1,1)) ... - & (datenum(num2str(MHW_here(:,2)),'yyyymmdd') <= datenum(year_here,12,31)); - judge_2=((full_mhw_start(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2),1)==year_here)) ... + judge_1=(MHW_here(:,1) >= datenum(year_here,1,1)) ... + & (MHW_here(:,2) <= datenum(year_here,12,31)); + judge_2=((full_mhw_start(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2),1)==year_here)) ... & ... - ((datenum(year_here,12,31)-datenum(num2str(MHW_here(:,1)),'yyyymmdd')) >= (datenum(num2str(MHW_here(:,2)),'yyyymmdd')-datenum(year_here+1,1,1))); - judge_3=((full_mhw_end(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2),1)==year_here)) ... + ((datenum(year_here,12,31)-(MHW_here(:,1))) >= (MHW_here(:,2))-datenum(year_here+1,1,1)); + judge_3=((full_mhw_end(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2),1)==year_here)) ... & ... - ((datenum(year_here-1,12,31)-datenum(num2str(MHW_here(:,1)),'yyyymmdd')) <= (datenum(num2str(MHW_here(:,2)),'yyyymmdd')-datenum(year_here,1,1))); + ((datenum(year_here-1,12,31)-(MHW_here(:,1))) <= (MHW_here(:,2))-datenum(year_here,1,1)); metric_judge=metric_here(judge_1 | judge_2 | judge_3); - if ~isempty(metric_judge) + if ~isempty(metric_judge); metric_judge=nanmean(metric_judge); else metric_judge=nan; @@ -115,9 +127,9 @@ case 'MaxInt' MHW=MHW{:,:}; - full_mhw_start=datevec(num2str(MHW(:,1)),'yyyymmdd'); - full_mhw_end=datevec(num2str(MHW(:,2)),'yyyymmdd'); - metric=MHW(:,4); + full_mhw_start=datevec(MHW(:,1)); + full_mhw_end=datevec(MHW(:,2)); + metric=MHW(:,5); period_used=datenum(data_start,1,1):(datenum(data_start,1,1)+size(mhw_ts,3)-1); period_used=datevec(period_used); @@ -129,28 +141,28 @@ p_metric=NaN(x,y); annual_metric=NaN(x,y,length(years_mhw)); - loc_full=unique(MHW(:,8:9),'rows'); + loc_full=unique(MHW(:,9:10),'rows'); - for m=1:size(loc_full) + for m=1:size(loc_full); loc_here=loc_full(m,:); - MHW_here=MHW(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2),:); - metric_here=metric(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2)); + MHW_here=MHW(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2),:); + metric_here=metric(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2)); mean_metric(loc_here(1),loc_here(2))=nanmean(metric_here); - for i=1:length(years_mhw) + for i=1:length(years_mhw); year_here=years_mhw(i); - judge_1=(datenum(num2str(MHW_here(:,1)),'yyyymmdd') >= datenum(year_here,1,1)) ... - & (datenum(num2str(MHW_here(:,2)),'yyyymmdd') <= datenum(year_here,12,31)); - judge_2=((full_mhw_start(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2),1)==year_here)) ... + judge_1=(MHW_here(:,1) >= datenum(year_here,1,1)) ... + & (MHW_here(:,2) <= datenum(year_here,12,31)); + judge_2=((full_mhw_start(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2),1)==year_here)) ... & ... - ((datenum(year_here,12,31)-datenum(num2str(MHW_here(:,1)),'yyyymmdd')) >= (datenum(num2str(MHW_here(:,2)),'yyyymmdd')-datenum(year_here+1,1,1))); - judge_3=((full_mhw_end(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2),1)==year_here)) ... + ((datenum(year_here,12,31)-(MHW_here(:,1))) >= (MHW_here(:,2))-datenum(year_here+1,1,1)); + judge_3=((full_mhw_end(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2),1)==year_here)) ... & ... - ((datenum(year_here-1,12,31)-datenum(num2str(MHW_here(:,1)),'yyyymmdd')) <= (datenum(num2str(MHW_here(:,2)),'yyyymmdd')-datenum(year_here,1,1))); + ((datenum(year_here-1,12,31)-(MHW_here(:,1))) <= (MHW_here(:,2))-datenum(year_here,1,1)); metric_judge=metric_here(judge_1 | judge_2 | judge_3); - if ~isempty(metric_judge) + if ~isempty(metric_judge); metric_judge=nanmean(metric_judge); else metric_judge=nan; @@ -170,9 +182,9 @@ case 'MeanInt' MHW=MHW{:,:}; - full_mhw_start=datevec(num2str(MHW(:,1)),'yyyymmdd'); - full_mhw_end=datevec(num2str(MHW(:,2)),'yyyymmdd'); - metric=MHW(:,5); + full_mhw_start=datevec(MHW(:,1)); + full_mhw_end=datevec(MHW(:,2)); + metric=MHW(:,6); period_used=datenum(data_start,1,1):(datenum(data_start,1,1)+size(mhw_ts,3)-1); period_used=datevec(period_used); @@ -184,28 +196,29 @@ p_metric=NaN(x,y); annual_metric=NaN(x,y,length(years_mhw)); - loc_full=unique(MHW(:,8:9),'rows'); + loc_full=unique(MHW(:,9:10),'rows'); - for m=1:size(loc_full) + for m=1:size(loc_full); loc_here=loc_full(m,:); - MHW_here=MHW(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2),:); - metric_here=metric(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2)); + MHW_here=MHW(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2),:); + metric_here=metric(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2)); mean_metric(loc_here(1),loc_here(2))=nanmean(metric_here); - for i=1:length(years_mhw) + for i=1:length(years_mhw); + year_here=years_mhw(i); - judge_1=(datenum(num2str(MHW_here(:,1)),'yyyymmdd') >= datenum(year_here,1,1)) ... - & (datenum(num2str(MHW_here(:,2)),'yyyymmdd') <= datenum(year_here,12,31)); - judge_2=((full_mhw_start(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2),1)==year_here)) ... + judge_1=(MHW_here(:,1) >= datenum(year_here,1,1)) ... + & (MHW_here(:,2) <= datenum(year_here,12,31)); + judge_2=((full_mhw_start(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2),1)==year_here)) ... & ... - ((datenum(year_here,12,31)-datenum(num2str(MHW_here(:,1)),'yyyymmdd')) >= (datenum(num2str(MHW_here(:,2)),'yyyymmdd')-datenum(year_here+1,1,1))); - judge_3=((full_mhw_end(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2),1)==year_here)) ... + ((datenum(year_here,12,31)-(MHW_here(:,1))) >= (MHW_here(:,2))-datenum(year_here+1,1,1)); + judge_3=((full_mhw_end(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2),1)==year_here)) ... & ... - ((datenum(year_here-1,12,31)-datenum(num2str(MHW_here(:,1)),'yyyymmdd')) <= (datenum(num2str(MHW_here(:,2)),'yyyymmdd')-datenum(year_here,1,1))); + ((datenum(year_here-1,12,31)-(MHW_here(:,1))) <= (MHW_here(:,2))-datenum(year_here,1,1)); metric_judge=metric_here(judge_1 | judge_2 | judge_3); - if ~isempty(metric_judge) + if ~isempty(metric_judge); metric_judge=nanmean(metric_judge); else metric_judge=nan; @@ -226,9 +239,9 @@ case 'CumInt' MHW=MHW{:,:}; - full_mhw_start=datevec(num2str(MHW(:,1)),'yyyymmdd'); - full_mhw_end=datevec(num2str(MHW(:,2)),'yyyymmdd'); - metric=MHW(:,7); + full_mhw_start=datevec(MHW(:,1)); + full_mhw_end=datevec(MHW(:,2)); + metric=MHW(:,8); period_used=datenum(data_start,1,1):(datenum(data_start,1,1)+size(mhw_ts,3)-1); period_used=datevec(period_used); @@ -240,28 +253,28 @@ p_metric=NaN(x,y); annual_metric=NaN(x,y,length(years_mhw)); - loc_full=unique(MHW(:,8:9),'rows'); + loc_full=unique(MHW(:,9:10),'rows'); - for m=1:size(loc_full) + for m=1:size(loc_full); loc_here=loc_full(m,:); - MHW_here=MHW(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2),:); - metric_here=metric(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2)); + MHW_here=MHW(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2),:); + metric_here=metric(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2)); mean_metric(loc_here(1),loc_here(2))=nanmean(metric_here); - for i=1:length(years_mhw) + for i=1:length(years_mhw); year_here=years_mhw(i); - judge_1=(datenum(num2str(MHW_here(:,1)),'yyyymmdd') >= datenum(year_here,1,1)) ... - & (datenum(num2str(MHW_here(:,2)),'yyyymmdd') <= datenum(year_here,12,31)); - judge_2=((full_mhw_start(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2),1)==year_here)) ... + judge_1=(MHW_here(:,1) >= datenum(year_here,1,1)) ... + & (MHW_here(:,2) <= datenum(year_here,12,31)); + judge_2=((full_mhw_start(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2),1)==year_here)) ... & ... - ((datenum(year_here,12,31)-datenum(num2str(MHW_here(:,1)),'yyyymmdd')) >= (datenum(num2str(MHW_here(:,2)),'yyyymmdd')-datenum(year_here+1,1,1))); - judge_3=((full_mhw_end(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2),1)==year_here)) ... + ((datenum(year_here,12,31)-(MHW_here(:,1))) >= (MHW_here(:,2))-datenum(year_here+1,1,1)); + judge_3=((full_mhw_end(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2),1)==year_here)) ... & ... - ((datenum(year_here-1,12,31)-datenum(num2str(MHW_here(:,1)),'yyyymmdd')) <= (datenum(num2str(MHW_here(:,2)),'yyyymmdd')-datenum(year_here,1,1))); + ((datenum(year_here-1,12,31)-(MHW_here(:,1))) <= (MHW_here(:,2))-datenum(year_here,1,1)); metric_judge=metric_here(judge_1 | judge_2 | judge_3); - if ~isempty(metric_judge) + if ~isempty(metric_judge); metric_judge=nanmean(metric_judge); else metric_judge=nan; @@ -298,12 +311,12 @@ p_metric=NaN(x,y); annual_metric=NaN(x,y,length(years_mhw)); - loc_full=unique(MHW(:,8:9),'rows'); + loc_full=unique(MHW(:,9:10),'rows'); - for m=1:size(loc_full) + for m=1:size(loc_full); loc_here=loc_full(m,:); - for i=1:length(years_mhw) + for i=1:length(years_mhw); year_here=years_mhw(i); mhw_here=squeeze(mhw_ts(loc_here(1),loc_here(2),(datenum(year_here,1,1):datenum(year_here,12,31))-datenum(data_start,1,1)+1)); annual_metric(loc_here(1),loc_here(2),i)=nansum(mhw_here~=0 & ~isnan(mhw_here)); @@ -325,8 +338,8 @@ MHW=MHW{:,:}; - full_mhw_start=datevec(num2str(MHW(:,1)),'yyyymmdd'); - full_mhw_end=datevec(num2str(MHW(:,2)),'yyyymmdd'); + full_mhw_start=datevec(MHW(:,1)); + full_mhw_end=datevec(MHW(:,2)); period_used=datenum(data_start,1,1):(datenum(data_start,1,1)+size(mhw_ts,3)-1); period_used=datevec(period_used); @@ -338,24 +351,24 @@ p_metric=NaN(x,y); annual_metric=NaN(x,y,length(years_mhw)); - loc_full=unique(MHW(:,8:9),'rows'); + loc_full=unique(MHW(:,9:10),'rows'); - for m=1:size(loc_full) + for m=1:size(loc_full); loc_here=loc_full(m,:); - MHW_here=MHW(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2),:); + MHW_here=MHW(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2),:); mean_metric(loc_here(1),loc_here(2))=(size(MHW_here,1)./length(years_mhw)); - for i=1:length(years_mhw) + for i=1:length(years_mhw); year_here=years_mhw(i); judge_1=(datenum(num2str(MHW_here(:,1)),'yyyymmdd') >= datenum(year_here,1,1)) ... & (datenum(num2str(MHW_here(:,2)),'yyyymmdd') <= datenum(year_here,12,31)); - judge_2=((full_mhw_start(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2),1)==year_here)) ... + judge_2=((full_mhw_start(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2),1)==year_here)) ... & ... - ((datenum(year_here,12,31)-datenum(num2str(MHW_here(:,1)),'yyyymmdd')) >= (datenum(num2str(MHW_here(:,2)),'yyyymmdd')-datenum(year_here+1,1,1))); - judge_3=((full_mhw_end(MHW(:,8)==loc_here(1) & MHW(:,9)==loc_here(2),1)==year_here)) ... + ((datenum(year_here,12,31)-(MHW_here(:,1))) >= (MHW_here(:,2)-datenum(year_here+1,1,1))); + judge_3=((full_mhw_end(MHW(:,9)==loc_here(1) & MHW(:,10)==loc_here(2),1)==year_here)) ... & ... - ((datenum(year_here-1,12,31)-datenum(num2str(MHW_here(:,1)),'yyyymmdd')) <= (datenum(num2str(MHW_here(:,2)),'yyyymmdd')-datenum(year_here,1,1))); + ((datenum(year_here-1,12,31)-(MHW_here(:,1))) <= (MHW_here(:,2)-datenum(year_here,1,1))); MHW_judge=MHW_here(judge_1(:) | judge_2(:) | judge_3(:),:);