Detecting storms from wave height data

Question

Maria Francesca bruno 2024-11-14，9:37

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2166269-detecting-storms-from-wave-height-data

评论： Star Strider 2024-11-19，13:46

H.mat

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));

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Image Analyst 2024-11-14，16:36

"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 2024-11-14，14:39

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2166269-detecting-storms-from-wave-height-data#answer_1545373

在 MATLAB Online 中打开

H.mat

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>

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>

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

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

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>

[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')

.

2 个评论
显示无隐藏无

Maria Francesca bruno 2024-11-19，13:42

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

Star Strider 2024-11-19，13:46

My pleasure!

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

.

请先登录，再进行评论。

Detecting storms from wave height data

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

2 个评论
显示无隐藏无

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Detecting storms from wave height data

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

2 个评论 显示 无隐藏 无

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

2 个评论
显示无隐藏无