Calculate the duration of each storm event

Hi all! I am trying to do some statistical anaylsis of precipitation and I would like to implement a code which do the following: I have hourly data of precipitation and what I would like to do is to calculate the duration of each storm. If there are 24 hours with no rainfall (value 0), the storm ends and I would like the code to start calculating the duration of the following storm. The resulting output should be a vector in which there are listed the durations of the different storms. Does anybody know how I could solve it? Thanks a lot Camilla

Answer 1

[Guillaume]

You haven't specified the criteria for the beginning of a storm, so I can only give you a generic solution. In any case, what you ultimately want to do is find the run length of a series of true / false (0 / 1). The first thing is thus to convert your precipitation data into 0/1, possibly like this:

precipitation = [0 0 0.2 0.4 1 1.3 1.1 0.2 0 0 0 0.3 2 2.1 1.9 0.2 0.1 0 0 0]; %made up data
isstorm = precipitation > 0.1; %for example

From there, a combination of diff and find is all you need to do to find the run length of the storms.

%ensure that we start and end with no storm
%so the first transition is always from sunny to storm 
%and last transition is always from storm to sunny
isstorm = [0 isstorm 0]; 

%find transitions between sunny and storm and back
beforetransition = find(diff(isstorm));

%elements 1:2:end are starts of storms
%elements 2:2:end are ends o storms
%duration is the difference
%because we made sure to start and end sunny, there's always an even number of transitions
stormduration = beforetransition(2:2:end) - beforetransition(1:2:end);

[Guillaume]

Camilla's comment copied here: Criteria of the storm is just precipitation >0 mm. My precipitation data are in mm s they range from 0 to around 60 mm...I could divide everything by 100 and try to code your solution! Thank you

No need to divide anything, just change the isstorm criteria to:

isstorm = precipitation > 0;

[Guillaume]

Camilla's comment copied here: But in your code I didn't understand how can I take into consideration the fact that a new storm has to be considered only if 24 hours of 0 mm of rain has occurred before the next one...so between one storm and the other there has to be at least 24 hr of 0 values

Sorry, I completely missed that. The simplest way to solve this is to remove the sunny periods that are too short. It's easy to calculate the duration of the sunny periods:

%duration of sunny periods
sunnyduration = beforetransition(3:2:end) - beforetransition(2:2:end-2);
%find the indices of sunnyduration that are too short:
sunnytooshort = find(sunnyduration < 24); 
%these indices * 2 are the starts that need to be removed
%and * 2 + 1 the ends that need to be removed
beforetransition([2*sunnytooshort, 2*sunnytooshort+1]) = [];

You can then calculate the storm duration as before:

stormduration = beforetransition(2:2:end) - beforetransition(1:2:end);

[Guillaume]

To sum up, this is the code

 precipitation = [
   rand(1, 12), zeros(1, 4), rand(1, 4), ... %1st storm duration = 20
   zeros(1, 26), ...
   rand(1, 4), zeros(1, 5), rand(1, 6), zeros(1, 23), rand(1, 12), ... %2nd storm duration = 50
   zeros(1, 24), ...
   rand(1, 14)]; %last storm duration = 14

isstorm = [0 precipitation > 0 0];
transitions = find(diff(isstorm));
tooshort = find(transitions(3:2:end) - transitions(2:2:end-2) < 24);
transitions([2*tooshort 2*tooshort+1]) = [];
stormduration = transitions(2:2:end) - transitions(1:2:end)

[Camilla Santicoli]

It works! Thank you very much for your help Guillaume, I was trying to solve it with index and diff but it wasn't working! Thanks again, very kind!