-
Notifications
You must be signed in to change notification settings - Fork 14
/
detectc.m
425 lines (356 loc) · 19.1 KB
/
detectc.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
function [MHW,mclim,m90,mhw_ts,category_ts]=detectc(temp,time,cli_start,cli_end,mhw_start,mhw_end,varargin)
%detect - Detecting spatial MHW/MCS
% Syntax
%
% [MHW]=detect(temp,time,cli_start,cli_end,mhw_start,mhw_end)
% [MHW,mclim,m90,mhw_ts]=detect(temp,time,cli_start,cli_end,mhw_start,mhw_end);
% [MHW,mclim,m90,mhw_ts]=detect(temp,time,cli_start,cli_end,mhw_start,mhw_end,'Event','MCS','Threshold',0.1);
%
% Description
%
% [MHW]=detect(temp,time,cli_start,cli_end,mhw_start,mhw_end) returns
% all detected MHW events for the m-by-n-by-t matrix TEMP starting in the
% year DATA_START. m, n and t separately indicate two spatial dimensions
% (m and n) and one temporal dimension (t). Climatologies used to
% determine events are calculated based on TEMP from CLI_START to
% CLI_END. MHW is a table where each row corresponds to a particular
% event during MHW_START to MHW_END and each column indicates a metric.
%
% [MHW,mclim,m90,mhw_ts]=detect(temp,time,cli_start,cli_end,mhw_start,mhw_end)
% also return the spatial climatology MCLIM (m-by-n-by-366) and threshold
% M90 (m-by-n-by-366) to calculate MHW events and resultant MHW time series
% MHW_TS (m-by-n-by-t). In the condition that there is no missing value in
% data TEMP, NaN in all outputs indicates lands and 0 in MHW_TS indicates
% the corresponding day in that grid is not in a MHW event.
%
% [MHW,mclim,m90,mhw_ts]=detect(temp,time,cli_start,cli_end,mhw_start,mhw_end,'Event','MCS','Threshold',0.1)
% returns MCS events based on 10th percentile threshold.
%
% Input Arguments
%
% temp - 3D daily temperature to detect MHW/MCS events, specified as a
% m-by-n-by-t matrix. m and n separately indicate two spatial dimensions
% and t indicates temporal dimension.
%
% time - datenum(start_year,start_month,start_day):datenum(end_year,
% end_month,end_day)
%
% cli_start - A numeric value in format of datennum(yyyy,mm,dd), indicating the start date for the period
% across which the spatial climatology and threshold are calculated.
%
% cli_end - A numeric value in format of datennum(yyyy,mm,dd) indicating the end year for the period across
% which the spatial climatology and threshold are calculated.
%
% data_start - A numeric value in format of datennum(yyyy,mm,dd) indicating the start year of your input
% data TEMP.
%
% mhw_start - A numeric value in format of datennum(yyyy,mm,dd) indicating the start year for the period
% across which MHW/MCS events are detected.
%
% 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.
%
% 'Event' - Default is 'MHW'.
% - 'MHW' - detecting MHW events.
% - 'MCS' - detecting MCS events.
%
% 'Threshold' - Default is 0.9. Threshold percentile to detect MHW/MCS
% events.
%
% 'windowHalfWidth' - Default is 5. Width of sliding window to calculate
% spatial climatology and threshold.
%
% 'smoothPercentileWidth' - Default is 31. Width of moving mean window to smooth spatial
% climatology and threshold.
%
% 'minDuration' - Default is 5. Minimum duration to accept a detection of MHW/MCS
% event.
%
% 'maxGap' - Default is 2. Maximum gap accepting joining of MHW events.
%
% 'ClimTemp' - Default is TEMP. The data used to calculate climatology
% and thresholds.
%
% 'ClimTime' - A vector of datenum() corresponding to ClimTemp.
%
% 'percentile' - Default is 'matlab'. Indicating the way to calculate
% percentile, using either 'matlab way' or 'python way'.
%
% Output Arguments
%
% MHW - A table containing all detected MHW/MCS events where each row
% 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_dur - duration of each event.
% - int_max - maximum intensity of each event.
% - int_mean - mean intensity of each event.
% - int_var - variance of intensity during each event.
% - int_cum - cumulative intensity across each event.
% - xloc - location of each event in x-dimension of TEMP.
% - yloc - location of each event in y-dimension of TEMP.
%
% mclim - A 3D matrix (m-by-n-by-366) containing climatologies.
%
% m90 - A 3D matrix (m-by-n-by-366) containing thresholds.
%
% mhw_ts - A 3D matrix
% (m-by-n-by-(datenum(MHW_end,1,1)-datenum(MHW_start)+1)) containing
% spatial intensity of MHW/MCS in each day.
%
% category_ts - A 3D matrix
% (m-by-n-by-(datenum(MHW_end,1,1)-datenum(MHW_start)+1)) containing
% category of MHW/MCS in each day.
% 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''');
%% "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
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
temp_clim=cat(3,NaN(size(temp_clim,1),size(temp_clim,2),ahead_date), ...
temp_clim);
elseif ahead_date<=0 && after_date>0
temp_clim=cat(3, ...
temp_clim,NaN(size(temp_clim,1),size(temp_clim,2),after_date));
else
end
temp_mhw=temp(:,:,time>=mhw_start & time<=mhw_end);
%% Calculating climatology and thresholds
date_true=datevec(cli_start-vWindowHalfWidth:cli_end+vWindowHalfWidth);
date_true=date_true(:,1:3);
date_false = date_true;
date_false(:,1) = 2012;
fake_doy = day(datetime(date_false),'dayofyear');
ind = 1:length(date_false);
mclim=NaN(size(temp,1),size(temp,2),366);
m90=NaN(size(temp,1),size(temp,2),366);
for i=1:366
if i == 60
else
ind_fake=ind;
ind_fake(fake_doy==i & ~ismember(datenum(date_true),cli_start:cli_end))=nan;
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);
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'
m90(:,:,i) = cellfun(@percentile, data_thre,repmat({vThreshold},size(temp,1),size(temp,2)));
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');
end
end
end
% Dealing with Feb29
m90(:,:,60) = mean(m90(:,:,[59 61]),3,'omitnan');
mclim(:,:,60) = mean(mclim(:,:,[59 61]),3,'omitnan');
% Does running averages of threshold and clim..
m90long=smoothdata(cat(3,m90,m90,m90),3,'movmean',vsmoothPercentileWidth);
m90=m90long(:,:,367:367+365);
mclimlong=smoothdata(cat(3,mclim,mclim,mclim),3,'movmean',vsmoothPercentileWidth);
mclim=mclimlong(:,:,367:367+365);
[x_size,y_size]=deal(size(m90,1),size(m90,2));
%mbigadd=temp_mhw;
date_mhw=datevec(mhw_start:mhw_end);
date_mhw(:,1)=2000;
indextocal = day(datetime(date_mhw),'dayofyear');
ts=str2double(string(datestr(mhw_start:mhw_end,'YYYYmmdd')));
mhw_ts=NaN(x_size,y_size,length(ts));
category_ts=NaN(x_size,y_size,length(ts));
MHW=cell(x_size,y_size);
%% Detecting MHW/MCS in each grid
switch vEvent
case 'MHW'
for i=1:x_size
for j=1:y_size
mhw_ts(i,j,isnan(squeeze(temp_mhw(i,j,:))))=nan;
if sum(isnan(squeeze(temp_mhw(i,j,:))))~=size(temp_mhw,3)
maysum=zeros(size(temp_mhw,3),1);
maysum(squeeze(temp_mhw(i,j,:))>=squeeze(m90(i,j,indextocal)))=1;
trigger=0;
potential_event=[];
for n=1:size(maysum,1)
if trigger==0 && maysum(n)==1
start_here=n;
trigger=1;
elseif trigger==1 && maysum(n)==0
end_here=n-1;
trigger=0;
potential_event=[potential_event;[start_here end_here]];
elseif n==size(maysum,1) && trigger==1 && maysum(n)==1
trigger=0;
end_here=n;
potential_event=[potential_event;[start_here end_here]];
end
end
if ~isempty(potential_event)
potential_event=potential_event((potential_event(:,2)-potential_event(:,1)+1)>=vminDuration,:);
if ~isempty(potential_event)
gaps=(potential_event(2:end,1) - potential_event(1:(end-1),2) - 1);
while min(gaps)<=vmaxGap
% potential_event(find(gaps<=vmaxGap),2)=potential_event(find(gaps<=vmaxGap)+1,2);
% potential_event(find(gaps<=vmaxGap)+1,:)=[];
% gaps=(potential_event(2:end,1) - potential_event(1:(end-1),2) - 1);
potential_event(find(gaps<=vmaxGap),2)=potential_event(find(gaps<=vmaxGap)+1,2);
loc_should_del=(find(gaps<=vmaxGap)+1);
loc_should_del=loc_should_del(~ismember(loc_should_del,find(gaps<=vmaxGap)));
potential_event(loc_should_del,:)=[];
gaps=(potential_event(2:end,1) - potential_event(1:(end-1),2) - 1);
end
mhwstart=NaN(size(potential_event,1),1);
mhwend=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);
mhwint_var=NaN(size(potential_event,1),1);
mhwint_cum=NaN(size(potential_event,1),1);
mhwcategory=NaN(size(potential_event,1),1);
for le=1:size(potential_event,1)
event_here=potential_event(le,:);
endtime=ts(event_here(2));
starttime=ts(event_here(1));
mcl=squeeze(mclim(i,j,indextocal(event_here(1):event_here(2))));
mth=squeeze(m90(i,j,indextocal(event_here(1):event_here(2))));
mrow=squeeze(temp_mhw(i,j,event_here(1):event_here(2)));
manom=mrow-mcl;
mca=mth-mcl;
mhw_ts(i,j,event_here(1):event_here(2))=manom;
[maxanom,~]=nanmax(squeeze(manom));
mhwint_max(le)=...
maxanom;
mhwint_mean(le)=...
mean(manom);
mhwint_var(le)=...
std(manom);
mhwint_cum(le)=...
sum(manom);
mhwstart(le)=starttime;
mhwend(le)=endtime;
mduration(le)=event_here(2)-event_here(1)+1;
mhwcategory(le)=nanmax(floor(manom./mca));
category_ts(i,j,event_here(1):event_here(2))=nanmax(floor(manom./mca));
end
mhwcategory(mhwcategory>4)=4;
MHW{i,j}=[mhwstart mhwend mduration mhwint_max mhwint_mean mhwint_var mhwint_cum repmat(i,size(mhwstart,1),1) repmat(j,size(mhwstart,1),1) mhwcategory];
end
end
end
end
end
case 'MCS'
for i=1:x_size
for j=1:y_size
mhw_ts(i,j,isnan(squeeze(temp_mhw(i,j,:))))=nan;
if sum(isnan(squeeze(temp_mhw(i,j,:))))~=size(temp_mhw,3)
maysum=zeros(size(temp_mhw,3),1);
maysum(squeeze(temp_mhw(i,j,:))<=squeeze(m90(i,j,indextocal)))=1;
trigger=0;
potential_event=[];
for n=1:size(maysum,1)
if trigger==0 && maysum(n)==1
start_here=n;
trigger=1;
elseif trigger==1 && maysum(n)==0
end_here=n-1;
trigger=0;
potential_event=[potential_event;[start_here end_here]];
elseif n==size(maysum,1) && trigger==1 && maysum(n)==1
trigger=0;
end_here=n;
potential_event=[potential_event;[start_here end_here]];
end
end
if ~isempty(potential_event)
potential_event=potential_event((potential_event(:,2)-potential_event(:,1)+1)>=vminDuration,:);
if ~isempty(potential_event)
gaps=(potential_event(2:end,1) - potential_event(1:(end-1),2) - 1);
while min(gaps)<=vmaxGap
% potential_event(find(gaps<=vmaxGap),2)=potential_event(find(gaps<=vmaxGap)+1,2);
% potential_event(find(gaps<=vmaxGap)+1,:)=[];
% gaps=(potential_event(2:end,1) - potential_event(1:(end-1),2) - 1);
potential_event(find(gaps<=vmaxGap),2)=potential_event(find(gaps<=vmaxGap)+1,2);
loc_should_del=(find(gaps<=vmaxGap)+1);
loc_should_del=loc_should_del(~ismember(loc_should_del,find(gaps<=vmaxGap)));
potential_event(loc_should_del,:)=[];
gaps=(potential_event(2:end,1) - potential_event(1:(end-1),2) - 1);
end
mhwstart=NaN(size(potential_event,1),1);
mhwend=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);
mhwint_var=NaN(size(potential_event,1),1);
mhwint_cum=NaN(size(potential_event,1),1);
mhwcategory=NaN(size(potential_event,1),1);
for le=1:size(potential_event,1)
event_here=potential_event(le,:);
endtime=ts(event_here(2));
starttime=ts(event_here(1));
mcl=squeeze(mclim(i,j,indextocal(event_here(1):event_here(2))));
mth=squeeze(m90(i,j,indextocal(event_here(1):event_here(2))));
mrow=squeeze(temp_mhw(i,j,event_here(1):event_here(2)));
manom=mrow-mcl;
mca=mth-mcl;
mhw_ts(i,j,event_here(1):event_here(2))=manom;
[maxanom,~]=nanmin(squeeze(manom));
mhwint_max(le)=...
maxanom;
mhwint_mean(le)=...
mean(manom);
mhwint_var(le)=...
std(manom);
mhwint_cum(le)=...
sum(manom);
mhwstart(le)=starttime;
mhwend(le)=endtime;
mduration(le)=event_here(2)-event_here(1)+1;
mhwcategory(le)=nanmax(floor(manom./mca));
category_ts(i,j,event_here(1):event_here(2))=nanmax(floor(manom./mca));
end
mhwcategory(mhwcategory>4)=4;
MHW{i,j}=[mhwstart mhwend mduration mhwint_max mhwint_mean mhwint_var mhwint_cum repmat(i,size(mhwstart,1),1) repmat(j,size(mhwstart,1),1) mhwcategory];
end
end
end
end
end
end
category_ts(category_ts>4)=4;
function p=percentile(data,thre)
if nansum(isnan(data))~=length(data)
data=data(~isnan(data));
pos=1+(length(data)-1)*thre;
p=interp1((1:length(data))',sort(data(:)),pos);
else
p=nan;
end