Detecting storms from wave height data

I'm using the Image Analyst code to detect sea storms from wave height data using a threshold value as suggested by https://it.mathworks.com/matlabcentral/answers/2119581-detection-of-storms-from-precipitation-data#answer_1459176.

I fixed a threshold for wave height and a minimum storm duration.

Now I would like to modify the code to include small holes in the subsamples (for example 2 or more missing data) or few values below threshold. I would like to prevent storm splitting (see attached figure).

Thanks a lot to Image Anayst for his support.

load ('H.mat'); %3 hour data
threshold=1 %
% Find time periods with H >= threshold
stormPeriods = bwconncomp(H >= threshold);
props = regionprops(stormPeriods, H, 'Area', 'MeanIntensity','MaxIntensity',"SubarrayIdx");
values = [props.Area];
props = props(values*3 > 12 ); % storms with duration >12 h
A_cell = (struct2cell(props));

[Image Analyst]

"include small holes in the subsamples (for example 2 or more missing data" <== So you want all NaN values to be considered as storms no matter how long the run of NaN's is?

"few values below threshold" <== like for example, what? 5 values below should be considered part of the storm on either side of that run of low values? 10 values? I guess we can just set a variable and you cann set it to whatever you want.

In your plot above, how many storms do you want there to be and where do they start and stop?

Answer 1

[Star Strider]

I am not certain what you want.

Try this —

load ('H.mat'); %3 hour data
whos('-file','H.mat')

  Name          Size             Bytes  Class     Attributes

  H         67200x1             537600  double              

threshold=1 %

threshold = 1

t = linspace(0, numel(H)-1, numel(H));                         % Supply Missing Time Vector

Storms = H >= threshold;

Stormsa = [Storms; false];
StormStart = strfind(Stormsa(:).', [0 1])+1;
StormEnd = strfind(Stormsa(:).', [1 0]);
StormDur = StormEnd - StormStart

StormDur = 1×2174

0 6 13 5 9 0 3 12 17 1 13 11 1 0 1 0 1 19 3 12 6 2 0 2 5 0 3 11 5 2

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

StormDur2 = StormDur(StormDur > 1)

StormDur2 = 1×1258

6 13 5 9 3 12 17 13 11 19 3 12 6 2 2 5 3 11 5 2 11 5 8 3 5 10 15 3 19 9

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

b = fitdist(StormDur2(:), 'exponential')

b = 
  ExponentialDistribution

  Exponential distribution
    mu = 7.67409   [7.26707, 8.11646]

StormDurStats = [min(StormDur) max(StormDur) mean(StormDur) median(StormDur) std(StormDur) mode(StormDur)]

StormDurStats = 1×6

0 44.0000 4.5892 2.0000 6.0523 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
histfit(StormDur2, 100, 'exponential')
grid

StormIdx = [StormStart(StormDur>1); StormEnd(StormDur>1)].';
StormSplitThreshold = mean(StormDur)

StormSplitThreshold = 4.5892

% StormIdx = StormIdx(1:end-1,:)
for k = 1:size(StormIdx,1)-1
    % DD = (StormIdx(k+1,1) - StormIdx(k,2))
    if (StormIdx(k+1,1) - StormIdx(k,2)) <= StormSplitThreshold
        StormIdx(k+1,1) = StormIdx(k,2);
    end
end

StormIdx

StormIdx = 1258×2

15 21 21 37 74 79 148 157 175 178 223 235 279 296 320 333 343 354 492 511

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[ts,Hs] = stairs(t, H);
format shortG

figure
stairs(t, H)
hold on
% patch([ts; flip(ts)], [zeros(size(Hs)); flip(Hs)], 'r', FaceAlpha=0.3, EdgeColor='r')
for k = 1:size(StormIdx,1)
    idx = ts >= t(StormIdx(k,1)) & ts <= t(StormIdx(k,2));
    [findidx1,findidx2] = bounds(find(idx));
    StormTimes(k,:) = [findidx1,findidx2,ts(findidx1),ts(findidx2)];
    % EndStormTimes = [k StormTimes(end,:)]
    patch([ts(idx); flip(ts(idx))], [zeros(size(Hs(idx))); flip(Hs(idx))], 'r', FaceAlpha=0.5, EdgeColor='none', EdgeAlpha=0)
    % plot(t(StormIdx(k,1) : StormIdx(k,2)), H(StormIdx(k,1) : StormIdx(k,2)), 'r.')
    AUC(k,:) = [ts(findidx1) ts(findidx2) trapz(ts(idx(1:2:end)), Hs(idx(1:2:end)))];
end

Results = array2table(AUC, 'VariableNames',{'Start Time','End Time','Area'})

Results = 1258x3 table
    Start Time    End Time    Area
    __________    ________    ____

        14           20       1.7 
        20           36       8.1 
        73           78       1.7 
       147          156       4.7 
       174          177       0.9 
       222          234       2.2 
       278          295         6 
       319          332       2.7 
       342          353       5.9 
       491          510       2.9 
       547          550       0.2 
       619          631       1.4 
       674          680       1.5 
       680          685       2.5 
       718          720       1.1 
       766          771       2.5 

% StormTimes

hold off
grid
xlim([0 1E+3])
yline(threshold)
xlabel('Time (Units)')
ylabel('Height')

.

[Maria Francesca]

Thank you very much, this suggestion fits perfectly my case.

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.