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));
1 个评论
Image Analyst
2024-11-14
"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?
采纳的回答
I am not certain what you want.
Try this —
load ('H.mat'); %3 hour data
whos('-file','H.mat')
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
StormDur2 = StormDur(StormDur > 1)
b = fitdist(StormDur2(:), 'exponential')
StormDurStats = [min(StormDur) max(StormDur) mean(StormDur) median(StormDur) std(StormDur) mode(StormDur)]
figure
histfit(StormDur2, 100, 'exponential')
grid
StormIdx = [StormStart(StormDur>1); StormEnd(StormDur>1)].';
StormSplitThreshold = mean(StormDur)
% 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
[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'})
% StormTimes
hold off
grid
xlim([0 1E+3])
yline(threshold)
xlabel('Time (Units)')
ylabel('Height')
.
2 个评论
Maria Francesca
2024-11-19
Star Strider
2024-11-19
My pleasure!
If my Answer helped you solve your problem, please Accept it!
.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Axis Labels 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
