% Laboratory in Oceanography - Data and Methods
% SMAST, UMass Dartmouth
%
% For going through simple examples of mapping toolbox
%
% Written by Miles A. Sundermeyer, 4/30/09

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% create an empty map axis
figure(1)
clf

worldmap world

% Import low-resolution world coastlines
whos -file coast.mat

% Load and plot the coastlines on the world map:
load coast
plotm(lat, long)

sum(isnan(lat))			% Note 'nan' separators

%%
% Reveal the underlying axis
axis on
set(gca,'XTickLabelMode','auto')
set(gca,'YTickLabelMode','auto')
set(gca,'xtick',[-10:10])
set(gca,'ytick',[-10:10])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% Create a new map axes for plotting data over Europe
h = worldmap('Europe')

% Inspect the map axes properties
mstruct = getm(h);
mstruct.mapprojection

% Add data to the map
geoshow('landareas.shp', 'FaceColor', [0.15 0.5 0.15])
geoshow('worldlakes.shp', 'FaceColor', 'cyan')
geoshow('worldrivers.shp', 'Color', 'blue')
geoshow('worldcities.shp', 'Marker', '.',...
                           'MarkerEdgeColor', 'red')
% Put a label on the map
labelLat = 35;
labelLon = 14;
textm(labelLat, labelLon, 'Mediterranean Sea')

%%
% Reveal the underlying axis
axis on
set(gca,'XTickLabelMode','auto')
set(gca,'YTickLabelMode','auto')
set(gca,'xtick',[-10:10])
set(gca,'ytick',[-10:10])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run through Boston aerial photo example
%%

% start a Map Viewer session by typing
mapview

% for ease of importing Mapping Toolbox demo data, set working directory to:
% (can also get here by navigating to 'toolbox/map/mapdemos')
cd(fullfile(matlabroot,'toolbox','map','mapdemos'))

% Select "Import From File" from File menu to open GeoTIFF file boston.tif

% Pointer coordinates are shown in the readout at the lower left, 
% e.g., as 774,114.36 feet easting (X), 2,955,685.56 feet northing (Y), 
% in Massachusetts State Plane coordinates.

% Import a vector data layer
boston_roads = shaperead('boston_roads.shp');

% Look at the information contained in this
boston_roads(1)

% 'unitsratio' function gives conversion factors between units of length.
%surveyFeetPerMeter = unitsratio('survey feet','meter');
%for k = 1:numel(boston_roads)
%    boston_roads(k).X = surveyFeetPerMeter * boston_roads(k).X;
%    boston_roads(k).Y = surveyFeetPerMeter * boston_roads(k).Y;
%end 

% Now map data that is already in the workspace using menu:
% "Import From Workspace > Vector Data >  Geographic Data Structure  
% from the File menu
% Specify boston_roads as the data to import from the workspace.

% Can designate any layer to be the active layer
% Doing so allows you to query attributes of the layer you select

% See the attributes for a vector layer using the Info tool

% Choose an attribute for the Datatip tool to inspect: Layers > boston_roads > Set Label Attribute

% Use the Datatip tool to identify attribute of any road displayed.

% Color roads according to their CLASS attribute, which takes on the values 1:6.
roadcolors = makesymbolspec('Line', ...
{'CLASS',1,'Color',[1 1 1]}, {'CLASS',2,'Color',[1 1 0]}, ...
{'CLASS',3,'Color',[0 1 0]}, {'CLASS',4,'Color',[0 1 1]}, ...
{'CLASS',5,'Color',[1 0 1]}, {'CLASS',6,'Color',[0 0 1]})
% ... then select menu: Layers > boston_roads > Set Symbol Spec.

% Add another layer, a set of points that identify 13 Boston landmarks.
boston_placenames = shaperead('boston_placenames.shp');
surveyFeetPerMeter = unitsratio('survey feet','meter');
%for k = 1:numel(boston_placenames)
%    boston_placenames(k).X = ...
%        surveyFeetPerMeter * boston_placenames(k).X;
%    boston_placenames(k).Y = ...
%        surveyFeetPerMeter * boston_placenames(k).Y;
%end
% Choose File > Import From Workspace > Vector Data > Geographic Data Structure

% To make markers more visually prominent, create symbolspec to represent 
% them as red filled circles:
places = makesymbolspec('Point',{'Default','Marker','o', ...
	'MarkerEdgeColor','r','MarkerFaceColor','r'})
% To activate this symbolspec: Layers > boston_placenames > Set Symbol Spec. 
% In the Layer Symbols dialog, highlight 'places'.

% Save select region as image file: File > Save As Raster Map > Selected Area
% Two files are created:
%   An image file (file.tif, file.png, or file.jpg)
%   A worldfile that georeferences the image (file.tfw, file.pgw, or file.jgw)

% Reload this as a georeferenced image:
[image cmap] = imread('mapping_MIT.tif');
R = worldfileread('mapping_MIT.tfw');
figure(1)
clf
mapshow(image, cmap, R);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
cd('C:\Documents and Settings\msundermeyer\Desktop\LabInOceanography\data+examples/Allens_Pond_mapping/hd25k_buzz')

mapview

% Import a vector data layer
hd25k_buzz_poly_MAS = shaperead('hd25k_buzz_poly_MAS.shp');