%% Frederik J Simons - Lab03 - fjsimonsl03
%% Basic data entry, and plotting
% Let us begin by downloading the *.mat file that I made using save, using the function websave, into a local file with the same name. And then let us use load to bring the variable alldata into our workspace. Notice how I handle the incoming struct using getfield. At the end of the two lines, I have a variable names GPS that has many, many rows, and seventeen columns with GPS locations and metadata. Note: these data were collected by students in previous years, and I spent some time (not enough!) editing them for outliers and nonsensical values.

websave('controlpoints.mat','http://geoweb.princeton.edu/people/simons/PIX/controlpoints.mat');
GPS=getfield(load('controlpoints.mat','alldata'),'alldata');

% You'll notice that load is much, much faster than webread. It ALWAYS pays to organize your data (and I do know what the values mean, and the columns in which certain data are being kept, and so on). Bottom line is: I can move along.

cpl=GPS(:,2);
easting=GPS(:,4);
northing=GPS(:,5);
elevation=GPS(:,6);
zyear=GPS(:,13);
numsats=GPS(:,16);
weather=GPS(:,17);

% Make a little plot, to get a quick idea of where we stand.

plot(easting,northing,'o')

%
%% In-class exercise 1: Remove the few outliers that are really obvious
% Go to last week's LiveScript, and figure out a clean and simple way to do this!

% Here is where you find the outliers
condisjun=
% Here is how you count them
nout=
% Here is how you reassign the variables after applying the logical condition
cpl=
easting=
northing= 
elevation= 
zyear= 
numsats= 
weather= 
% Update the plots
plot(easting,northing,'o')

%
%% In-class exercise 2: Split the data into control-point categories
% Go to last week's LiveScript, and figure out a clean and simple way to do this! Use a for loop. Reassign the variables to a cell array. Please note: there is no control point 3 in this set! Also, the tags 5 refer to a line.

% Here is where you do that, with a for loop, into a cell
for index=[1 2 4 5]
    % Here is the condition
    condisjun=
    % Here is how you assign them to the entries of a cell array.
    eastcpl{index}=easting(
    northcpl{index}= 
    elevcpl{ 
    zyrcpl{ 
    numscpl{ 
    weathcpl{ 
end

%
%% In-class exercise 3: Plot the data per control point type
% Go to last week's LiveScript, and figure out a clean and simple way to do this!

% Here is where you do that
% Assign colors and symbols
cols={'k','y','b','r','c'};
symb={'o','v','^','s','d'};
% Use a for loop
for index=[1 2 4 5]
    % Inside of the loop, supply the plot command
     plot( ,...
        'MarkerFaceColor',cols{index},'MarkerEdgeColor',cols{index},...
        'Marker',symb{index},'LineStyle','none')
    hold on
end 
hold off

%
%% Make a histogram of the data for control point 2
% First, figure out the counts, using hist. Next, produce a bar graph. That's altogether a histogram, although I recommend avoiding that specific function. Use subplot to get multiple panels.

% Produce the counts
% This is the index of the control point whose histogram you make
index=2;
% This is the number of bins that you will be using
nbins=35;
% The first output of "hist" is the value, the second is the count
% The first index (=1) into the cells means "easting"
% The second index (=index) into the cells means "cpl"
[ycpl{1,index},xcpl{1,index}]=hist(eastcpl{index},nbins);
% The first index (=2) into the cells means "northing"
% The second index (=index) into the cells means "counts"
[ycpl{2,index},xcpl{2,index}]=hist(northcpl{index},nbins);
% Make the histogram, which is a bar graph of the counts
% Hence, this is the histogram of the eastings in questions
subplot(1,2,1)
% The last entry into the function "bar" is the width of the bars
b(1)=bar(xcpl{1,index},ycpl{1,index},1);
xl(1)=xlabel(sprintf('easting of control point %i',index));
% Hence, this is the histogram of the eastings in questions
subplot(1,2,2)
b(2)=bar(xcpl{2,index},ycpl{2,index},1);
xl(2)=xlabel(sprintf('northing of control point %i',index));

%
%% In-class exercise 4: Prettify the histograms
% Label, title, colors, reasonable limits, etc. Label the axes. Report the mean and the stdev, using sprintf in the title.
%
%% Determine a regression line for control line 5
% Begin by figuring out the value of the correlation coefficient, using corr. Determine the coefficients of the regression line using polyfit. Overlay the regression line using polyval.

clf
% Keep it simple, rename the variables "x" and "y"
x=eastcpl{5};
y=northcpl{5};
% Plot the values using the simplified variable names
plot(x,y,'o')
% The correlation coefficient
rho=corr(x,y);
% The coefficients of the regression line (a straight line!)
N=2;
[P,S,MU]=polyfit(x,y,N);
% The evaluated regression line... N points suffice!
xvals=linspace(min(x),max(x),2);
yvals=polyval(P,xvals,[],MU);
% Plot it right on top of there
hold on
plot(xvals,yvals,'k-')
hold off

%
%% In-class exercise 5: Remove outliers, redetermine the regression line

% Freshen up your variables, in case you're doing it wrong
x=eastcpl{5};
y=northcpl{5};
% Here is where you find some of the the outliers
condisjun= 
x=x(condisjun);
y=y(condisjun);
plot(x,y,'o')
% The correlation coefficient
rho=corr(x,y);
% The coefficients of the regression line (a straight line!)
N=2;
[P,S,MU]=polyfit(x,y,N);
% The evaluated regression line... N points suffice!
xvals=linspace(min(x),max(x),2);
yvals=polyval(P,xvals,[],MU);
% Plot it right on top of there
hold on
plot(xvals,yvals,'k-')
hold off
axis tight

%
%% In-class exercise 6: Prettify the regression plot
% Label, title, axis, further outliers, and so on. Label the axes. Report the correlation coefficient using sprintf in the title.
%
%% Your job
% Tell us a story about this data set, one that involves making at least one histogram, and one cross-plot, including a simple regression analysis and reporting on the correlation coefficient.