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.

[Maria Francesca bruno]

Wonderful solution. It works perfectly.

How can I modify the code to include any small holes in the subsample (for example 2 missing data) ? I wold like to prevent storm splitting.

Thank you

[Image Analyst]

If it's missing, I assume the value in the vector would be nan. I'm also going to assume that there is data on each side of the nan runs (so it's not zero and the nan run is inside a storm). So then you can use isnan to find indexes where the nans live. Now you have to decide what value you want those nans to take. You could set them to 0.0001 - some small non-zero value :

nanLocations = isnan(rain);
rain(nanLocations) = 0.0001;

OR you could use interp1 to set those values to interpolated values across the nan run. Can you do it? Untested code to get you started.

nanLocations = isnan(rain);
if ~isempty(nanLocations)
    rain2 = rain; % Make a copy
    xFull = 1 : numel(rain2);
    xNoNan = xFull(~nanLocations); % Indexes of valid locations only.
    rain2 = rain2(xNoNan); % Extract only the valid values, not nans.
    rain_corrected = interp1(rain2, xFull); % Or something like that - check documentation!
end

If not, let me know how you want to handle the missing values.

Now that the nan/missing locations are filled in with valid values, the storm(s) will not be split into separate storms.

[Maria Francesca bruno]

Thank you very much, for missing values the idea is good.

For my reasearch I'm looking for a procedure that is able to handle both NANs and values below a theshold.

For example i want to extract events with values greater than a fixed threshold including also few values that are smaller than the threshold (if present).

I have used bwconncomp for a timeseries named H

stormPeriods = bwconncomp(H >= threshold);

props = regionprops(stormPeriods, H, 'Area', 'MeanIntensity','MaxIntensity',"SubarrayIdx");

But the selected storms are often splitted by NaNs and small clusters of data less than the threshold. How can I fix it? Possibly with a single procedure.

[Image Analyst]

I thought I answered that but maybe not since it was just off the top of my head because you're not submitting your data so I can't test/run the code. Like I asked, please start a new question with your data.

If you have any more questions, then attach your data and code to read it in with the paperclip icon after you read this:

TUTORIAL: How to ask a question (on Answers) and get a fast answer