Detection of storms from precipitation data
6 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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.
2 commentaires
the cyclist
le 16 Mai 2024
Can you upload the data (or a representative sample)? You can use the paper clip icon in the INSERT section of the toolbar.
Réponses (2)
the cyclist
le 16 Mai 2024
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)
2 commentaires
the cyclist
le 17 Mai 2024
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.
Image Analyst
le 17 Mai 2024
Modifié(e) : Image Analyst
le 17 Mai 2024
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')
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)
% 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)
% 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))
0 commentaires
Voir également
Catégories
En savoir plus sur Data Type Conversion dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!