Calculate duration of Peaks using findpeaks function

Dear reader,

I am currently working on an analysis of storm events. I am interested in finding the duration of storms defined by windspeeds that exceed a particular threshold:

Currently, I have used the findpeaks() function to identify the location (timestamps) and maximum intensity (peak value) of the windspeeds. I have done this as follows:

[PKS,LOCS] = findpeaks(Windspeed,time,'MinPeakHeight',th)

Here, windspeed is a vector containing hourly windspeed measurements, time is vector with corresponding datenum values for the measurements of the windspeed. Finally, th is the thresholds for when windspeeds are considered a storm events.

Now, I also want to know the duration of the storm events. That is, the time between when the windspeed first exceeds th until it agains drops below th.

I have tried using 'WidthReference', however, this only gives the options 'halfprom' and 'halfheight' and both don't give the values I am looking for.

Can anyone help me with this. Many thanks in advance!

Answer 1

[Star Strider]

Try something like this —

t = linspace(0, 24*30-1, 2*30);                                 % Simulate Time
wv = 2 + sin(2*pi*t/30) .* cos(2*pi*t/24);                      % Simulate Data

th = 2.4;                                                       % Set Threshold
idx = find(diff(sign(wv - th)));                                % Approximate Indices Where 'wv' Crosses 'th'
for k = 1:numel(idx)
    idxrng = max([1,idx(k)-1]) : min([numel(t),idx(k)+1]);      % Index Range For Interpolation
    t_th(k,:) = interp1(wv(idxrng), t(idxrng), th);             % More Precise Time Estimates For 'th' Crossings
end
% t_th
t_th = t_th(1:2*floor(numel(t_th)/2));                          % Delete Incomplete Events
t_thm = reshape(t_th, 2, []);                                   % Matrix Of Paired Events
EventDuration = diff(t_thm)

EventDuration = 1×4

45.9247 34.8986 33.6824 45.7171

figure
plot(t, wv)
hold on
plot(t, ones(size(t))*th, ':k')
plot(t_thm, ones(size(t_thm))*th, '|-r')
hold off
grid
text(mean(t_thm), ones(size(mean(t_thm)))*th, compose('Duration = %.2f\\rightarrow', EventDuration), 'Horiz','right', 'Vert','bottom', 'Rotation',-60)

It may be necessary for you to adapt this to your data. It may work as written, however it should be straightforward to change it if required.

.

[Michiel Smit]

Thanks for the suggestion. Took me a bit of time to get it working, but in the end works as intented.

[Star Strider]

As always, my pleasure!

(Having a sample of the data to work with always helps when coding the answer.)

[Michiel Smit]

I noticed a mistake in the script above. If wv(idxrng) has succesive values that are the same (very rare but it may happen), interp1() will not work. I have solved it by adding the following lines to the for-loop for k:

% Find the threshold crossings for time and wind velocity
cross_wv = wv(idxrng);
cross_t = time(idxrng);

% Make sure all points are unique for interp1 to function
[~,id_unique] = unique(cross_wv);

% Now Approximate crossing time by running interp1()
t_th(kk,:) = interp1(cross_wv(id_unique), cross_t(id_unique), th);

Then proceed as normally as shown in the example script above.

It does not give the most accurate results, but for my work, accurate to the nearest 4 hours is more than sufficient.

[Star Strider]

The way I solve that is usually to interpolate to a higher resolution (more data points).

In a timetable, that can be done with the retime function, using the 'pchip' method (my favourite, unless the data are truly linear), otherwise using interp1 (again with 'pchip'). This has the disadvantage of creating data where none previously existed, however here, where the storm durations are the desired result and the time tolerance is forgiving of such inaccuracies, it likely would be appropriate and straightforward.

The problem with unique is that it elimiinates some of the data. That likely creates more problems than it solves.

Answer 2

[David Hill]

f=num2str(Windspeed>=th);
f=f(f~=' ');
[startIndex,endIndex]=regexp(f,'1*');
Duration=diff(datenum([time(startIndex),time(endIndex)])');