Matlab Routines

You can copy the text and paste it into a MatLab .m file. Then save the file, giving it the name of the corresponding function, such as 'sergeiinterpolate.m'

function  [out] = sergeiinterpolate (inarr)
% This function interpolates a 1D array by substututing the NaN data points by preceding values
d=size(inarr);
prev=inarr(1);
arr=inarr;
for i=2:d
    if isnan(inarr(i))==1
        arr(i)=prev;
    else
        prev=inarr(i);
    end
end   
out = arr;

% end of function


Sequences of commands to process the  Lake Temperatures data (see also functions below)

% Convert the "date" columns in the raw data into a single "julian day" column 
jd0=juliandate(2002,7,4);    % set starting date, look up the correct date in your data file!
jd=juliandate(data(:,1),data(:,2),data(:,3)) - jd0;
temp=data(:,4);   %temp = original temperature data
%%

%Remove NaNs
inarr(:,1)=jd;
inarr(:,2)=temp;
Torig=sergeiRemoveNaNs(inarr,2);     % The output is Torig (jd,temp)= original data without NaNs
%%

% Merge multiple years into one
Tnew(:,1)=mod(Torig(:,1),365);
Tnew(:,2)=Torig(:,2);
Tnewsorted=sortrows(Tnew,1);
Torig=Tnewsorted;
Tbackup=Tnewsorted;


Routines for the analysis of the lake temperature data (copy them into a .m file and save in Matlab directory)

function  [out] = sergeiRemoveNaNs (inarr,dim)
% This function scans a two-column array along the dimension dim and removes the rows
% that contain NaNs
d=max(size(inarr));
k=1;
for i=1:d
    if isnan(inarr(i,dim))==0
        arr(k,:)=inarr(i,:);
        k=k+1;
    end
end   
out = arr;

% end of function




function  [out] = sergeiRunning (inarr,dim,ws)
% This function scans a two-column array (julian day, temperature) along the dimension dim (1 or 2) and
% computes the running mean and running variance. The averaging window ws should be specified
% in terms of time (e.g., number of days).  
% For each averaging interval, the function records the number of data
% points that fall within that interval.
% The "julian day" column in the input array should contain values between
% 0 and 365.
% The output is an array with 365 rows and columns that correspond to:
% 1) julian day 2) running mean 3) number of points found within the
% averaging window 4) running variance 5) "error on the mean" 

d=max(size(inarr));
% make the data cyclic by extending the sereies for another year 
extended=cat(1,inarr,inarr);
jd=extended(:,1);
jd(d+1:2*d,1)=jd(1:d)+365;
T=extended(:,2);

%cycle over the days of the year
i=1;
for m=1:365
    %find the closest matching days of the year in the data array
    while jd(i)>m
        i=i-1;
    end
    while jd(i)<m
        i=i+1;
    end
    %calculate running mean by summing up data within the averaging window
    k=0;sum=0;count=0;
    while jd(i+k)<=jd(i)+ws
        sum=sum+T(i+k);
        count=count+1;
        k=k+1;
    end
    time(m)=jd(i);
    rmean(m)=sum/count;
    npoints(m)=count;
end

% now use the running mean to calculate the running variance
i=1;
for m=1:365
    while jd(i)>m
        i=i-1;
    end
    while jd(i)<m
        i=i+1;
    end

    k=0;sum2=0;count=0;
    while jd(i+k)<=jd(i)+ws
        sum2=sum2+(T(i+k)-rmean(m))^2;
        k=k+1;
    end
   
    if npoints(m)<=1
        rstd(m)=NaN;
    else
        rstd(m)=sqrt(sum2/(npoints(m)-1));
    end

end  

% center the running values to the middle of the averaging window
time(1:365)=mod(time(1:365)+ws/2,365);

% form the output array
out(:,1)=time(1:365);
out (:,2) = rmean(1:365);
out(:,3)=npoints(1:365);
out(:,4)=rstd(1:365);
out(:,5)=rstd(1:365)/sqrt(npoints(1:365));

% end of function