Detection of storms from precipitation data

I have 5 min intervals of precipitation data that I am trying to sum into separate storm events.

6 hours of no precipitation (or 72 zero values) signify the end of the storm event.

I'm trying to create a code that

a) Give me the array of matrices (diifferent names e.g M1,M2 etc) of the storm events which lies between 6 hours of no precipitaion.

b) Sums the precip values until the 6 hours of no rainfall and then starts summing the following storm event. The output will be an array with individual storm event depths of the time series.

Kindly help.

[the cyclist]

Can you upload the data (or a representative sample)? You can use the paper clip icon in the INSERT section of the toolbar.

[Joydeb Saha]

I am giving you an idea on it. It is a structure file(S). It has time in one folder.. upto see below

'01-Jan-2000 00:00:00'

'01-Jan-2000 00:05:00'

'01-Jan-2000 00:10:00'

'01-Jan-2000 00:15:00'

'01-Jan-2000 00:20:00'

.....

31-Dec-2021 23:55:00

And 'vals'

Which is precipitation values.

Both are in different folders.

I have extracted those, and sending you the vals folder. Time folder is larger than limit that I can send here.

Answer 1

[the cyclist]

I think this does what you want. It is easy to make minor indexing mistakes, so you should double-check:

% Load data
load("VE_0234_5min_vals.mat")

% Find all the sequences of at least 72 consecutive zeros. However, this
% will overcount. (For example, a string of exactly 73 zeros will be counted twice.)
zeroStartRedundant = strfind(l',zeros(1,72))';

% Find the redundant sequences (because they are one ahead of the prior sequence)
redundantSequenceIndex = [true; diff(zeroStartRedundant)==1];

% Remove the redundant ones
zeroStart = [1; zeroStartRedundant(not(redundantSequenceIndex))];

% Count the number of storms
numberStorms = numel(zeroStart);

% Preallocate memory for the storm total
stormTotal = zeros(numberStorms,1);

% For each storm, sum the rainfall
for ns = 1:numberStorms-1
    stormTotal(ns) = sum(l(zeroStart(ns):(zeroStart(ns+1)-1)));
end

% Display a histogram of the storm totals
figure
histogram(stormTotal)

[Joydeb Saha]

Dear the Cyclist,

This is OK, but some modification is needed for the further development.

I have reduced the data size, extracted only Oct, 2000 data values here in the new attachment.

1) I need to have the storm matrix separately, that is.. say from row 125 to row 196 you have values geater than zero (numeric values), before that, you have 72 or more than zero values and after that you have 72 or more than 72 zero values. so we need to extract out that matrix separately, there may be many such matrices detecting these storms.

2) Whether this storm is erosive or not, to find it out.. we need to check sum of these storms, and make another folder. If the sum of these storms is greater than 12.7 then it is erosive.

3) There are further calculation also. That we can discuss after the development of the code up to that.

Hope my description is clear. Can you kindly help me on this part ..

Thank you so much in advance.

[the cyclist]

I didn't understand #1 or #2.

The way you phrase #3, you are assuming that I am willing to take on a long development of code with you (or for you). You should not assume that. Instead, you should post individual questions for each part of the task that you do not understand how to do yourself. Then you can put it all together into your project.

Answer 2

[Image Analyst]

If you have the Image Processing Toolbox, this is really trivial because it has function specially built for this kind of thing. All you need to do is call bwareaopen to find periods of where there are 72 or more time points with zero rain. Then simply call regionprops to get the summed rain fall in the non-dry regions. It's only 2 or 3 lines of code.

dryPeriods = bwareaopen(rain == 0, 72);

props = regionprops(~dryPeriods, rain, 'Area', 'MeanIntensity');

amountOfRainPerStorm = [props.Area] .* [props.MeanIntensity];

You'll have a vector, amountOfRainPerStorm, that for each element gives you the total amount of rain summed over the rainy period (i.e. between 72+ runs of dry periods). Here is a full demo.

s = load('VE_0234_5min_vals.mat')

s = struct with fields:
    l: [2314369x1 double]

rain = s.l;
% Find time periods with no rain for at least 6 time periods
dryPeriods = bwareaopen(rain == 0, 72);
% OPTIONAL: Find the number of dry periods.
[labeledDryPeriods, numDryPeriods] = bwlabel(dryPeriods);
fprintf('I found %d dry periods.\n', numDryPeriods)

I found 2271 dry periods.

% The other time periods, ~dryPeriods, has periods with rain
% or else periods shorter than 6 periods between rain times.
% OPTIONAL: Find the number of rainy periods.
[labeledRainyPeriods, numRainyPeriods] = bwlabel(~dryPeriods);
fprintf('I found %d rainy periods.\n', numRainyPeriods)

I found 2270 rainy periods.

% We want the sum of each of those rainy periods.
rainyPeriods = ~dryPeriods;
props = regionprops(rainyPeriods, rain, 'Area', 'MeanIntensity');
amountOfRainPerStorm = [props.Area] .* [props.MeanIntensity];
% ALL DONE!

%==========================================================================
% Now display results nicely and informatively.
% Show rain per storm in plot (though it will be subsampled to fit all 6794
% points on your screen.
subplot(2, 1, 1);
plot(amountOfRainPerStorm, 'b.');
grid on;
title('Rainfall per rainy period')
xlabel('Storm Number');
ylabel('Summed rainfall');
% Show histogram of the rain amounts.
subplot(2, 1, 2);
histogram(amountOfRainPerStorm);
grid on;
title('Histogram of rainfalls per storm')
xlabel('Summed rainfall');
ylabel('Number of storms with this total rainfall')

fprintf('The longest rainy period had %f of rain.\n', max(amountOfRainPerStorm))

The longest rainy period had 158.600000 of rain.