%% Frederik J Simons - Lab02 - fjsimonsl02
%% Basic data entry, and plotting
% Let us begin by reading in the Comma-Separated-Value data table, using the function webread:

GPS=webread('http://geoweb.princeton.edu/people/simons/PIX/controlpoints.csv');

% What did it actually read in? Take a look at what the function fieldnames returns:

fn=fieldnames(GPS);

% The strings containing the table headers are inside an array of the class cell, e.g.:

fn{4}
fn{5}

% So now we know how to address the entries in the data table: they are available via explicit addressing as:

easting=GPS.East;

% Alternatively, we can pull them out of the data table via an implicit addressing, as:

northing=GPS.(fn{5});
elevation=GPS.(fn{6});

% We can now proceed to making a very simple plot, of the variable "easting" versus the variable "northing".

plot(easting,northing,'o')
xlabel('easting')
ylabel('northing')

%
%% Selection of the data we really want to plot
% Somewhat messy! We are not using all of the information in the table. For starters, the second column of the table appears to have information about which type of marker every GPS location (easting, northing, elevation) refers to. Let's pull out that information, again using the explicit form of addressing:

cpl=GPS.CP_L_;

% Let us work with just the first four categories of control points, and plot those. I'm going to do this by reassigning the variables "easting", "northing" and "elevation" to only contain those for which the "cpl" variable takes the values 1, 2, 3, or 4, and then replot!

easting=easting(cpl<5);
northing=northing(cpl<5);
elevation=elevation(cpl<5);
plot(easting,northing,'o')
xlabel('easting')
ylabel('northing')

% Don't forget that now you have reassigned everything, you also need to shrink the variable "cpl"! You must do this last, since you're using it up until this point.

cpl=cpl(cpl<5);

%
%% Obvious outlier removal, and reporting thereof
% Excellent. We now appear to have collected just the four control points. But what do we see? Massive outliers. Unless somebody's been moving the beacons, there's definitely two points out there that don't belong. In this case, it's easy enough to see which points we should excise. They are the two only points with northings between 4.4635 and 4.4664. We'll use the function find to track down those points, to count how many there are (our eyes say: two), and then we will remove them from our data. First, report on how many there are:

condisjun=northing<4.4664e6 & northing>4.46635e6;
nout=sum(condisjun)

% The variable condisjun is of the class logical. It gets a 1 when the condition is met, a 0 when it isn't. So the sum of the ones is the number of hits, namely, two. Now we can negate the condition (converting all the 0s to 1s and the 1s to 0s), and we'll have the logical array that will KEEP the entries into the vector that we really want to keep. These, we reassign - for the third time - to the variables "easting", "northing" and "elevation".

easting=easting(~condisjun);
northing=northing(~condisjun);
elevation=elevation(~condisjun);
cpl=cpl(~condisjun);

% And then we plot again!

plot(easting,northing,'o')
xlabel('easting')
ylabel('northing')

% Never forget that if you remove any outliers, you MUST tell. Best to put this information in a title. We will make a little string using the function sprintf, and feed the result to the title command.

tstring=sprintf('Campus Survey, %i data after removing %i outliers',...
    length(easting),nout);
title(tstring)

%
%% Distinguishing different data by colored symbols
% We have neglected to visualize one vital piece of information. The four clusters refer to four different control points. We'd definitely like to distinguish them, via colors AND symbols (there are color blind people out there!) Below, I show you how to do such a thing, making a bare-bones plot.
% I find it advisable to pick colors and symbols first! There will be four in total, and we collect them in a cell array. Don't pick red and green together, if you can avoid it. The list of possible choices is under help plot.

cols={'k','y','b','r'};
symb={'o','v','^','s'};

% Next, using logical operators, we'll identify the four different categories. Time for a for loop!

for index=1:4
    plot(easting(cpl==index),northing(cpl==index),...
        'MarkerFaceColor',cols{index},'MarkerEdgeColor',cols{index},...
        'Marker',symb{index},'LineStyle','none')
    hold on
end
hold off

% And it wouldn't be complete without a legend. Now you need to know what those things actually were.

legend('Guyot','Wright','Cannon','Statue')

% Don't forget to put the labels and the title on! And slap some grid lines on it.

xlabel('easting')
ylabel('northing')
title(tstring)
grid on

%
%% Proper referencing of the data
% The axes of the plot remain hard to read. Who has time to decipher all these numbers, and multiply by the powers of 10, etc? We better pick a reference point from which "easting" and "northing" will simply be distances east and north, in meters.
% Let the data determine what's a good reference! How about the lower-left corner of the data span, i.e. the minimum northing and the minimum easting.

mineast=min(easting);
minnorth=min(northing);

% We better reassign the variables.

easting=easting-mineast;
northing=northing-minnorth;

% And we owe it to the gentle reader to tell us what the reference point is. And what the units are. How about we put that into the axis labels themselves. This time we will collect "axis handles" to the x and the y labels, so we can change their font sizes, later on. Same for the title.

for index=1:4
    plot(easting(cpl==index),northing(cpl==index),...
        'MarkerFaceColor',cols{index},'MarkerEdgeColor',cols{index},...
        'Marker',symb{index},'LineStyle','none')
    hold on
end
hold off
xl=xlabel(sprintf('easting (m) from %i',mineast));
yl=ylabel(sprintf('northing (m) from %i',minnorth));
tl=title(tstring); 
grid on

%
%% Data-driven beautification of your plot
% I'm finding these axes a little, well, tight. Let's force the axis limits to be "opened up", how about, by 10% of the data span. First, by setting the axes to be really, really tight.

axis tight

% Then, collecting what those axis limits are.

xli=xlim;
yli=ylim;

% And then adding a percentage of the data range of space on either side.

xperc=5; xspace=round([xli(2)-xli(1)]*xperc/100);
yperc=5; yspace=round([yli(2)-yli(1)]*xperc/100);
xlim(xli+[-1 1]*xspace);
ylim(yli+[-1 1]*yspace);

% Seems a little roundabout, but it's very robust. You've let the data themselves decide where the reference point is, and what a reasonable amount of space is for the box around them. This will always work. Personally, I don't like inward-facing ticks. Let's change that. The variable "gca" refers to the handle of the current axis.

set(gca,'TickDir','out','TickLength',[0.02 0.025]/2)

% And then, for the final embellishment, change the fontsizes to something readable, and well-proportioned. You can deal your values to all the variables on the left-hand side.

[xl.FontSize,yl.FontSize,tl.FontSize]=deal(15);

% Last thing, I promise. I do not much like how the automatic placement of the title makes it hug the box so much. Move it! We did just make a reasonable spacing variable, "yspace", so we reuse it. If I divide by two, the title won't run off the page. (If I want it any higher, I'd have to shrink the plot, but that is for another day).

tl.Position=tl.Position+[0 yspace/2 0];

%
%%
%
%% Your job... which we've all but done for you!
% Make a NEW LiveScript that plots the location of your rocks, color-coded by sedimentary, igneous or metamorphic, properly titled and with comments for every line that you use in the way that we've built the examples above. Save that as yournetidl02 and submit to BlackBoard. Here's some thoughts to start.
% First, pick some colors and symbols. Collect what will be fed to the legend and title.

% Prep work for labeling and annotation
cols={'r','b','g'};
symb={'o','v','s'};
legs={'sedimentary','igneous','metamorphic'};
tstring={'fjsimonsl02 Campus Rocks by rock type'};

% Assemble your data set. Rather than reading from any presaved filed, you could hand-type them below, for which again I suggest creating a cell.

% Sedimentary rock locations, my random examples (you supply your own!)
lox{1}=[529300     4466206;
        528978     4466250;
        528951     4466467];
% Igneous rock locations, my random examples (you supply your own!)
lox{2}=[528971     4466440;
        528976     4466440];
% Metamorphic rock locations, my random examples (you supply your own!)
lox{3}=[528980     4466248;
        529083     4466537;
        529075     4466541;
        528951     4466461];

% Figure out the referencing points! You'll need the function cat to collect all of the locations together, then use min and max to find the corners.

% Determine the minimum easting and northings across all rock types
mineastnorth=min(cat(1,lox{:}));
mineast=mineastnorth(1);
minnorth=mineastnorth(2);

% Now you can complete the prep work by making axis label strings ahead of time.

% Further prep work for labeling and annotation
xstring=sprintf('easting (m) from %i',mineast);
ystring=sprintf('northing (m) from %i',minnorth);

% And then, finally, you can do what the code ACTUALLY does, which is plot! Instead of reassigning the location to new, referenced variables, you can perform the subtraction inside of the plot command.

% For loop with the plot commands for each of the rock types
clf
for index=1:length(lox)
    % lox{index} contains a rock type of index=1,2,3...
    % The FIRST column is their eastings, the SECOND column their northings
    % We reference to a sensible origin by subtracting the reference values inside of the plot command
    plot(lox{index}(:,1)-mineast,lox{index}(:,2)-minnorth,...
        'MarkerFaceColor',cols{index},'MarkerEdgeColor',cols{index},...
        'Marker',symb{index},'LineStyle','none')
    hold on
end 
hold off 
grid on

% Inspect your plot. Looking about right? Then the labels, legend, and title are alrady there, and just need to be added.

% Use the prep-work label, legend and title strings to make the annotations all at the very end
tl=title(tstring);
xl=xlabel(xstring);
yl=ylabel(ystring);
lg=legend(legs);

% Pleased? Now open up the axes just a tad, as we have done in the example.

% Final cosmetics
axis tight
% Old axis ranges, made by MATLAB
xli=xlim; yli=ylim;
% Percentage by which you are going to open them up a bit
xperc=5; yperc=5; 
% New axis ranges, made by YOU
nxli=xli+[-1 1]*round([xli(2)-xli(1)]*xperc/100);
nyli=yli+[-1 1]*round([yli(2)-yli(1)]*xperc/100);
% Apply those new axis ranges
xlim(nxli); ylim(nyli);
% Clean up
set(gca,'TickDir','out','TickLength',[0.02 0.025]/2)
[xl.FontSize,yl.FontSize,tl.FontSize]=deal(15);
tl.Position=tl.Position+[0 yspace/2 0];

% Want to get massively fancy? Load the database of campus buildings, and add those to the plot

% Get the necessary files of the webpage
websave('CampusBuildings.zip','http://geoweb.princeton.edu/people/simons/PIX/CampusBuildings.zip');
% Now, leave MATLAB, make sure that the locally stored CampusBuildings.zip gets unzipped!
!unzip -o CampusBuildings.zip
% Then read it in
S=shaperead('CampusBuildings/buildings_reprojected.shp');
% And add it to the plot... but remember you need to subtract the reference!
hold on
plot([S.X]-mineast,[S.Y]-minnorth, 'k-') 
hold off