How to find consecutive values above a certain threshold?

Hi.

I'm picking out values from a hourly data set, and I want to pick out values that are above 12 at least three consecutive times. Any tips on how I can do this?

example: A=[0 1 2 5 7 8 13 17 28 11 6 0 2 1 4]

I want to put the three values 13 17 28 into a vector.

Help is greatly appreciated!

- Kristine

Answer 1

[Azzi Abdelmalek]

A=[0 1 2 5 7 8 13 17 28 11 6 0 2 1 40 55 88 47 4  44  ]
idx=A>12;
ii1=strfind([0 idx 0],[0 1]);
ii2=strfind([0 idx 0],[1 0])-1;
ii=(ii2-ii1+1)>=3;
out=arrayfun(@(x,y) A(x:y),ii1(ii),ii2(ii),'un',0);
celldisp(out)

[Kristine]

Thank you!

[Azzi Abdelmalek]

v=[39904 0 9;39904 0.0417 7;39904 0.0833 13;39904 0.1250 14;39904 0.1667 16;39904 0.2083 10  ]
A=v(:,3)'
idx=A>12;
ii1=strfind([0 idx 0],[0 1]);
ii2=strfind([0 idx 0],[1 0])-1;
ii=(ii2-ii1+1)>=3;
out=arrayfun(@(x,y) v(x:y,:),ii1(ii),ii2(ii),'un',0);
celldisp(out)

[Kristine]

Thank you very much!!

[Kristine]

I want to put all the values from out into one matrix. But I can't seem to save the values in my loop. Do you have any suggestions to how I can fix this?

out=1x119 cell

for i=1:length(out)

D=out{1,i}

end

Last questions promise :p

- Kristine

[Image Analyst]

Kristine:

Try this:

allValues = [out{:}]

Then see the FAQ.

[Kristine]

I used allValues = vertcat(out{:}) Thank you again!!

[James]

Can you explain how this solution works?

[Muhammad Shahid Iqbal]

Hello guys, this was very helpful but how about if we also need the index of those consective values which are greater than threshold???

[Image Analyst]

Muhammad: See my solution below. Just ask regionprops for PixelIdxList or PixelList - that will be the indexes.

Answer 2

[Image Analyst]

Kristine:

This pretty easy and straightforward if you have the Image Processing Toolbox to identify stretches where the numbers are above the threshold and measure their lengths. Then just save those stretches of numbers into cells of a cell array.

% Define sample data and a threshold value.
A=[0 1 2 5 7 8 13 17 28 11 6 0 2 1 4 91 49 37 79 9 100 101 3]
threshold = 12;
% Find logical vector where A > threshold
binaryVector = A > 12
% Label each region with a label - an "ID" number.
[labeledVector, numRegions] = bwlabel(binaryVector)
% Measure lengths of each region and the indexes
measurements = regionprops(labeledVector, A, 'Area', 'PixelValues');
% Find regions where the area (length) are 3 or greater and
% put the values into a cell of a cell array
for k = 1 : numRegions
  if measurements(k).Area >= 3
    % Area (length) is 3 or greater, so store the values.
    ca{k} = measurements(k).PixelValues;
  end
end
% Display the regions that meet the criteria:
celldisp(ca)

In the command window, this is what you'll see:

A =
     0     1     2     5     7     8    13    17    28    11     6     0     2     1     4    91    49    37    79     9   100   101     3
binaryVector =
     0     0     0     0     0     0     1     1     1     0     0     0     0     0     0     1     1     1     1     0     1     1     0
labeledVector =
     0     0     0     0     0     0     1     1     1     0     0     0     0     0     0     2     2     2     2     0     3     3     0
numRegions =
     3
ca{1} =
    13    17    28
ca{2} =
    91    49    37    79

[Kristine]

Thank you!

[Ahmad Bayhaqi]

Hi @Image Analyst. It's good to see your script. Thank you. I also get the same condition here. But how to loop regionprops in matrix?

So I have a=76x80,60266 and b=76x80x60266. ; 76=lon, 78=lat and 60266 is time

I want to take values 'a' that are above 'b' threshold with the 5 consecutive days. However, I want to do it for my all regions. What I have done is I did several loops to follow your script.

first loop : to get labeled vector and numregions for all regions

second loop: to get the measurement value using regionprops.

binaryvector=a>b;
for i=1:1:length(lon)
    for j=1:1:length(lat)
     [labeledvector(i,j,:), numregions(i,j,:)]=bwlabeln(binaryvector(i,j,:)); 
 end
end

for i=1:length(lon)
    for j=1:length(lat)
        measurements(i,j,:)=squeeze(regionprops(labeledvector(i,j,:),sst(i,j,:),'Area','PixelValues'));
    end
end

first loop is success. However, I got stuck for using regionprops in loop?

Do you have any idea? and also how to do loop to get ca?

Thank you

[Image Analyst]

I believe it would go something like this:

mask = a > b;
for i=1:length(lon)
    for j = 1:length(lat)
        thisSlice = mask(i, j, :);
        [labeledMatrix, numRegions] = bwlabel(thisSlice);
        props = regionprops(labeledMatrix, 'Area', 'PixelValues')
        % Now do something with props....
    end
end

[Ahmad Bayhaqi]

Thank you @Image Analyst but, in my data, the threshold in every grid is also different, so the output would be different props in every location.

The code that you showed just like producing the one grid.

so, I tried this

mask = a > b;

for i=1:length(lon)
    for j = 1:length(lat)
        thisSlice(i,j,:) = mask(i, j, :);
        [labeledMatrix(i,j,:), numRegions(i,j,:)] = bwlabeln(thisSlice(i,j,:));
        props(i,j,:) = regionprops(labeledMatrix(i,j,:),a(i,j,:),'Area','PixelValues')
    end
end

but, it ends up with the error in the props. It always said about dimension error.

Do you have any idea?

Thank you

Answer 3

[Jan]

With FEX: RunLength:

A        = [0 1 2 5 7 8 13 17 28 11 6 0 2 1 4];
[B, N]   = RunLength(A > 12);
B(N < 3) = false;
mask     = RunLength(B, N);
Result   = A(mask);

[Kristine]

I tried running what you wrote above, but it gave me the error: Undefined function 'RunLength' for input arguments of type 'logical'. Could it be that I don't have that function or because I need to change the vector from double to something else? Thank you for helping!!

[Image Analyst]

You'd need to download that "RunLength" function from the File Exchange using the link he gave you.

[Kristine]

Oh I didn't see the link. Thank you Image Analyst!!

Answer 4

[Stephen23]

Here is a conceptually very simple method (I changed the third value to 20 as well, to provide a sequence of values > 12 but shorter then three):

>> N = 12;
>> A = [0,1,20,5,7,8,13,17,28,11,6,0,2,1,4];
>> idx = [true,A(1:end-1)>N] & A>N & [A(2:end)>N,true];
>> idx = [false,idx(1:end-1)] | idx | [idx(2:end),false];
>> A(idx)
ans =
    13    17    28

[Kristine]

Thank you!

[Kristine]

Okay I'm probably being annoying right now, but how do I do it if my time series looks like this:

A = 504223x3 double:

39904 0 9

39904 0.0417 7

39904 0.0833 13

39904 0.1250 14

39904 0.1667 16

39904 0.2083 10

...

[Azzi Abdelmalek]

Kristine, What is your question?

[Kristine]

Same question as above, but I also need to include the date and hour for my values above 12. Does that make sense?

Answer 5

[Lane Foulks]

A=[0 1 2 5 7 8 13 17 28 11 6 0 2 1 4 14 18 0 2 16 15 18 13 0 ] % added a passing and failing block above threshold
Adiff = diff(A>12);
ind_start = find(Adiff==1);
ind_stop = find(Adiff==-1);
block_length = ind_stop-ind_start; % list of consecutive section lengths
blocks_ind = find(block_length>2)% list of blocks above min length
for ii = 1:numel(blocks_ind) % loops through each block
    A(ind_start(blocks_ind(ii))+1:ind_stop(blocks_ind(ii)))
end