% This code is an example that shows how to calculate a  
% linear fit of target station data to reference station data
% given one year of data for each, and then use that linear fit
% to estimate the target data from the reference data and 
% compare it to the actual target data.

% You will need to modify this code for Problems 2 and 3.

% For Problem 2, you will need to change the year and select the 
% stations assigned to you (see link on class web page).

% For Problem 3, you will need to read in all years from 1997 to 2010 
%  of the reference station data. See the code linked to 
% 'Matlab code for reading daily average wind speeds for 1997-2011'
% in Lecture 11 if you need help to do this. You will then use the
% reference station daily average wind speeds and your linear fit to 
% estimate the daily average wind speeds for the target station,
% and from those, calculate the annual mean wind speed, just 
% as you did for Problem 1 (a).

% You need to specify the year, the target and reference stations, and
% the group that they are in. 
% Note that Zach and David will have to read in data for two groups.

% station IDs:
% group 1: 92 44 52
% group 2: 29 108 50

% construct file name for each file

year = 97
    
yy = mod(year, 100);all

yyS = num2str(yy,'%02i');

% ---select file name for group one or group two---

fn = ['wspd_day_avg_groupone_' yyS '.mat'];

%fn = ['wspd_day_avg_grouptwo_' yyS '.mat'];

load (fn)

disp(' ')
disp(['year = ' yyS])

% wspd_day_avg            365x3

% least squares linear fit 

% P = POLYFIT(X,Y,N) finds the coefficients of a polynomial P(X) of
%     degree N that fits the data Y best in a least-squares sense. P is a
%     row vector of length N+1 containing the polynomial coefficients in
%     descending powers, P(1)*X^N + P(2)*X^(N-1) +...+ P(N)*X + P(N+1).


% ---specify your station pair---

i = 2; % reference station
j = 1; % target station 

% find days for which data is missing for one or the other stations

mday1 = isnan(wspd_day_avg(:,i) );  % missing for reference
mday2 = isnan(wspd_day_avg(:,j) );  % missing for target

mday = mday1 |  mday2;  % missing for reference or target

missing_days = sum(mday)

% x = reference station data

x = wspd_day_avg(~mday,i);
%x = wspd_day_avg(:,i);

% y = target station data

y = wspd_day_avg(~mday,j);
%y = wspd_day_avg(:,j);

% calculate coefficients for linear fit

n = 1; % linear fit

p = polyfit(x,y,n)

% calculate fit to target data from reference data using fit
% yfit = p(1) * x + p(2)

yfit = polyval(p,x);

corrcoef(y,yfit)    % should be same as correlation between y and x
std(y)              % std dev of target data
std(yfit)           % std dev of fit to target data

% check result by plotting time series of target and fit

figure (1)

plot(y)
grid on

hold on

plot(yfit,'-r')
legend('target', 'fit to target')
xlabel('day of year')
ylabel('wind speed')

hold off

% check result by plotting scatter plot of target and fit
figure (2)

scatter(x,y)
grid on
xlabel('reference station wind speed')
ylabel('target station wind speed')

hold on

scatter(x,yfit,'r')

legend('target', 'fit to target')

hold off