Error using trapz to calculate area of each event under the peaks
2 次查看(过去 30 天)
显示 更早的评论
Kasih Ditaningtyas Sari Pratiwi
2017-11-6
评论: Kasih Ditaningtyas Sari Pratiwi
2017-11-7
Hi. I need your help. I try to calculate area of each peaks event, but I always got an error.
%%Peak flow analysis 2 (Figure 3)
flowtable = finalCSVnew(:,[1,7:8]); % create table containing DateAndTime, Durchflusslm, and SummeaktuellerTagm data
peakflowEvent = flowtable{:,2} > 3 ; % determine the threshold of flow(m3/h) for peakflowEvent
% use false as logical vector to determine transition. With function diff,
% transitions from false (0) to true (1) will be 1 and transitions from true
% to false will be -1. This will be 1 at the start of a dry period and -1 after the end
peakTransitions = diff([false; peakflowEvent; false]);
eventStarts = find(peakTransitions == 1);
eventEnds = find(peakTransitions == -1) -1;
% define the peak flow of each event through the flow data (peakflow) and
% the time when peak flow is happened (peakflowtime)
[peakflow, peakflowlocrel] = arrayfun(@(s, e) max(flowtable.Durchflusslm(s:e)), eventStarts, eventEnds);
peakflowlocabs = peakflowlocrel + eventStarts - 1;
peakflowtime=flowtable.DateAndTime(peakflowlocabs);
% create result table containing start and end time for peak flow event, the duration
% between start and end time, and peak flow
peakflowanalysis2 = table(flowtable.DateAndTime(eventStarts), flowtable.DateAndTime(eventEnds), ...
flowtable.DateAndTime(eventEnds) - flowtable.DateAndTime(eventStarts), ...
peakflow, peakflowtime, ...
'VariableNames', {'Start', 'End', 'Duration','PeakFlow','PeakFlowTime'});
volume=trapz(flowtable.Durchflusslm(flowtable.DateAndTime(eventStarts):flowtable.DateAndTime(eventEnds)));
It got an error on the last function to determine volume with trapz "Input must be scalars". Do you know how to properly use trapz for calculating area of each peak event? I have been trying since morning to solve this, but I always got an error. I would really appreciate for your help.
2 个评论
Walter Roberson
2017-11-6
You did not post a copy of the error message, and you did not post a copy of your data for us to test with.
采纳的回答
Walter Roberson
2017-11-6
In the line
volume=trapz(flowtable.Durchflusslm(flowtable.DateAndTime(eventStarts):flowtable.DateAndTime(eventEnds)));
you have
flowtable.DateAndTime(eventStarts):flowtable.DateAndTime(eventEnds)
eventStarts and eventEnds are vectors, so flowtable.DateAndTime(eventStarts) and flowtable.DateAndTime(eventEnds) are vectors. It is not possible to use vectors on either side of a ":" operator.
You also have the problem that by default the ":" operator uses an increment of 1. In the context of datetime objects, an increment of 1 corresponds to one day.
What I would suggest would be to not do find() in creating eventStarts and eventEnds -- to just construct a single logical mask of the values you want included and then use that logical mask to do the extraction of Durchflusslm
Note, though, that the setup you were trying to use, and logical masks, would share the properties that the non-adjacent regions would be pushed beside each other and the trapz would be calculated as if they had always been contiguous. I would suggest to you that it would make more sense to calculate the area of each region separately. The difference in outputs of these two cases would correspond to 1/2 the sum of each of the boundary values. (Interior points in trapazoid rule are weighted twice as high as the exterior points.)
3 个评论
Walter Roberson
2017-11-7
Suppose you have
1111___1111
then if you trapz each group separately, you would get 1/2*1+1+1+1/2*1 = area 3 for the first group, and again 3 for the second group, for a total of 6. But if you squish them together by removing the unwanted part in the middle, to get 11111111 then you would get 1/2*1+1+1+1+1+1+1+1/2*1 = 7. Therefor, squishing the groups together overestimates the area.
If you want a logical mask for which areas are in peak flow, then it would simply be your peakflowEvent mask.
trapz(flowtable.Durchflusslm(peakflowEvent))
would have the effect of squishing the areas together and doing the trapz like your code currently tries to calculate. But I am suggesting that you should instead be looping over your segments and calculating the areas separately:
nEvents = length(eventStarts);
volumes = zeros(1, nEvents);
for K = 1 : nEvents
volumes(K) = trapz(flowtable.Durchflusslm(eventStarts(K):eventEnds(K));
end
volume = sum(volumes);
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!