Selecting a rain event from data series

Hi everyone. I hope you can help me. I have a time series of data: rain data (h in mm) recorded every 10 min in a day (column vector of 144 data). From the data I must select an event preceded and followed by 6 hours (36 data) without rain (h=0). How do I select the data range?

[Dyuman Joshi]

A scalar data value recorded every 10 min for a day would give 6*24 == 144 data points per day. (6 data point an hour * 24 hours in a day)

How many data points do you collect each hour? If 1, then how did you get 1440 elements? or did you collect data for 10 days?

"From the data I must select an event preceded and followed by 6 hours (36 data) without rain (h=0). How do I select the data range?"

Using conditional operations and logical indexing - Find Array Elements That Meet a Condition

[Me]

Sorry, the zero is a typo. In one hour I have 6 logged data, so in one day 144 data. How do I select the data?

Answer 1

[Alexander]

I think this should work now:

clear;
load rain2.txt; rain = rain2; % Only to reuse the code above
NoOfDry = 1; 
rHmm = rain(:,2);rHmm = rHmm(:)'; % To force a row vector
tH = rain(:,1)/6; % Only to have a timeline in hours, doesn't needed here
figure(1);plot(rHmm);grid;
dryTime = zeros(1,36); % 6h dry, 36*10 minutes
yDry = strfind(rHmm, dryTime); % Now we are looking for the indices dry >= 6h
if(isempty(yDry))
    fprintf('In your data set are NO period(s) with ge. 6h no rain\n');
    NoOfDry = 0;
    return;
end
yDryDiff = diff(yDry); % 1 means > 6h dry
iNoOffDry = find(yDryDiff > 1); % Here are the positions of dry >= 6h
NoOfDry = NoOfDry+length(iNoOffDry); % This is the number of occurences dry >= 6h
fprintf('In your data set are %i period(s) with ge. 6h no rain\n',NoOfDry);
commandwindow

[Me]

A thousand thanks! You solved my problem!

I'm looking for >= 6h dry, and I didn't know how to work on indices.

Thanks again, good job!

[Alexander]

Cheers

[Alexander]

I think it's a good idea to accept this answer if it fulfills your needs. ;-)

Answer 2

[Alexander]

You can use strfind:

y = strfind(data, [0 0 0 0 0 0])

[Me]

Thanks for your quick reply! "strfind" is not good for my problem.

I am attaching a txt file to better explain the problem. In the txt file you will find the indices in the first column, the rainfall height in the second. The data that I have to select are between the index 37 and 105, these data are preceded by 6 hours without rain (at least 36 data in which h=0) and are followed by 6 hours without rain (at least 36 data in which h = 0). How can I do?

Answer 3

[Alexander]

I have to contradict, strfind is perfect if you want to avoid loops. Have a look at my rapid programmed script:

load rain.txt
tH = rain(:,1)/6; tH = tH(:)';% I prefere hours
rHmm = rain(:,2);rHmm = rHmm(:)'; % only for convenience
tHobs = tH(36:108); % I modified the observation time: 06:00 - 18:00
rHmmObs = rHmm(36:108);
figure(1);plot(tH,rHmm);grid;

dryTime = zeros(1,36); % 6h dry
figure(2);plot(tHobs,rHmmObs);grid;

yDry = strfind(rHmmObs, dryTime)

yDry =

     []

% Result: yDry = [], which means that no 6h w/o rain between 6-18:00
% Let's reduce the dry time to one hour:
dryTime = zeros(1,6);
yDry = strfind(rHmmObs, dryTime);
disp(yDry)

    15    16    17    18    19    20    21    22    23    24    25    26    51    60    61    62    63    64

% There are a lot of consecutive numbers in, which means that there are dry
% periods which are longer than one hour. Let's test these:
yDryDiff = diff(yDry);
disp(yDryDiff)

     1     1     1     1     1     1     1     1     1     1     1    25     9     1     1     1     1

% This means, you have one period of (6+11)*10 minutes w/o rain
%             you have 2 periods with exact one hour  w/o rain
%             you have one period of (6+4)*10 minutes w/o rain
% To check item 1:
dryTime = zeros(1,17);
yDry = strfind(rHmmObs, dryTime);
disp(yDry)

    15

[Alexander]

Thanks

[Dyuman Joshi]

Alexander, your code still does not find the indices which OP wants to obtain.

[Alexander]

Please explain.

strfind([10, 20, 300, 23570],0)

ans = []

Answeres correct.

strfind([10, 20, 300, 23570],[0 0])

Also.

Do you have an example for my understanding?

[Dyuman Joshi]

My bad, I misunderstood something.

[Me]

Thanks a lot for the replies. Unfortunately, I can't solve my problem, probably because I don't know how to use matlab very well. Could you give me some further help?

In both suggested codes the indexes of the data I need are imposed, but in reality I don't know the position of the data I need, and I have to automate the problem because my data set will consist of data recorded every 10 min for at least 10 years.I have verified that for more than two days the codes do not work.I am attaching a new txt file. In this case "observing" the data, I single out 3 groups of data preceded and followed by 6 dry hours (h=0). I need to automatically know where these 3 data groups are.

[Alexander]

I think I don't understand your problem. Are you looking for exact 6h dry (=0)? Or >= 6h dry? Just some code:

clear;
load rain2.txt; rain = rain2; % Only to reuse the code above

Error using load
Unable to find file or directory 'rain2.txt'.

rHmm = rain(:,2);rHmm = rHmm(:)'; % To force a row vector
tH = rain(:,1)/6; % Only to have a timeline in hours, doesn't needed here
figure(1);plot(rHmm);grid;
dryTime = zeros(1,36); % 6h dry, 36*10 minutes
yDry = strfind(rHmm, dryTime); % Now we are looking for the indices dry >= 6h
yDryDiff = diff(yDry); % 1 means > 6h dry
iNoOffDry = find(yDryDiff > 1); % Here are the positions of dry >= 6h
NoOffDry = length(iNoOffDry); % This is the number of occurences dry >= 6h
fprintf('In your data set are %i period(s) with ge. 6h no rain\n',NoOffDry)
commandwindow

If you want to get the exact time or index you have just work off the indices backward.

Just a very small piece of advice: If you really want to use the large advantages of Matlab, you should get a little bit more familiar with it. ;-)

[Alexander]

There is an exception I've forgotten: the first matching pattern. But I'm currently not at my PC. Hence, I will update my code tomorrow. Sorry.