Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Detect_edits2 #13

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Detect_edits2 #13

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
387dad2
Detect_edits2
maximemarin Oct 2, 2019
59966ae
Delete detect.m
maximemarin Nov 14, 2019
7b50612
Rename detect01.m to detect.m
maximemarin Nov 14, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
120 changes: 72 additions & 48 deletions detect.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
%
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added mhw_peak for both MHW and MCS algorithms

% 'Event' - Default is 'MHW'.
% - 'MHW' - detecting MHW events.
% - 'MCS' - detecting MCS events.
Expand Down Expand Up @@ -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.
Expand All @@ -99,46 +103,55 @@
% (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''');
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Modified the input variables parsing so that users don't need the Stats toolbox

end

%% "What if cli_start-window or cli_end+window exceeds the time range of data"

ahead_date=time(1)-(cli_start-vWindowHalfWidth);
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
Expand Down Expand Up @@ -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'
Expand Down Expand Up @@ -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;

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because mhw dates are now in datenums, coding can be simplified

mhw_ts=zeros(x_size,y_size,length(ts));

Expand Down Expand Up @@ -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);
Expand All @@ -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)=...
Expand All @@ -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
Expand Down Expand Up @@ -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);
Expand All @@ -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)=...
Expand All @@ -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
Expand All @@ -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),...
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indexing numbers had to be modified due to the addition of a variable.

'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


39 changes: 28 additions & 11 deletions event_line.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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{:});
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Once again, this only works if the user has the Stats toolbox. I've removed it similarly to detect, so that users can use it without the toolbox

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),:);
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Adding peak_date as an output of the detect function, there is one more dimension into the table, to the numbers had to be modified accordingly

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)];
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For some reason, this gave me an error so I modified it


period_plot_v=datevec(period_plot);
period_unique=datevec(datenum(2016,1,1):datenum(2016,12,31));
Expand All @@ -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);
Expand All @@ -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);
Expand Down
Loading