Error using trapz to calculate area of each event under the peaks

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.

[Walter Roberson]

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.

[Kasih Ditaningtyas Sari Pratiwi]

Sorry, I've just updated the copy of my data and an error message "Inputs must be scalar". Thanks Walter.

Answer 1

[Walter Roberson]

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

[Kasih Ditaningtyas Sari Pratiwi]

Thank you for your reply Walter. So I better to combine eventStarts and eventEnds as one logical mask? Do you know how to do it other than find?

I still did not get about calculating the area separately, so do you mean I have to filter each region so that I can calculate it with trapz?

[Walter Roberson]

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

[Kasih Ditaningtyas Sari Pratiwi]

Thanks Walter! :D