%% Exercise 4: Location Map of ODP 659 Benthic Oxygen Isotope Record

%% SECTION 6.3 ETOPO1 DATA SET
% Next ETOPO1 data set
%NCOLS  1201
%NROWS  1201
%XLLCENTER  -30.00000
%YLLCENTER  10.00000
%CELLSIZE 0.01666667
%NODATA_VALUE  -32768
clear
ETOPO2 = load('grid01.asc');

%%
ETOPO2 = flipud(ETOPO2);
ETOPO2(find(ETOPO2 == -32768)) = NaN;

%%
[LON,LAT] = meshgrid(-30:1/60:-10,10:1/60:30);

%%
% Simple graph from MATLAB Recipes - 3rd Edition
surf(LON,LAT,ETOPO2)
shading interp
view(0,90)
colorbar

hold on
plot(-21.0262,18.0772,...
    'Marker','o','MarkerSize',10,...
    'MarkerFaceColor','r','MarkerEdgeColor','k')
text(-21.0262-5,18.0772,'ODP 659',...
    'FontSize',20,'Color','r')

%%
% We can create a version of the graph compatible with the graphics
% in the previous chapters, i.e., x-Axis 10 cm wide, line thickness of axis 
% 0.6 pt, line thickness of curve 0.5 pt, axis labels 8 pt sans serif,
% title 10 pt sans serif. NOTE: set current axis property VIEW to [0 90].
figure1 = figure('Color',[1 1 1]);
axes1 = axes('Units','Centimeters',...
    'Box','on',...
    'LineWidth',0.6,...
    'FontSize',8);
surf1 = surf(LON,LAT,ETOPO2,...
    'Parent',axes1,...
    'FaceColor','interp',...
    'EdgeColor','none');
title1 = title('ETOPO2 Data Set',...
    'FontSize',10);
xlabel1 = xlabel('Longitude',...
    'FontSize',8);
ylabel1 = ylabel('Latitude',...
    'FontSize',8);
set(gca,'View',[0 90])

hold on
% 8°04.63'N, 21°01.57'W
plot3(-21.0262,18.0772,interp2(LON,LAT,ETOPO2,-21.0262,18.0772)+100,...
    'Marker','o','MarkerSize',10,...
    'MarkerFaceColor','r','MarkerEdgeColor','k')
text(-21.0262-4.5,18.0772,'ODP 659',...
    'FontSize',20,'Color','r')

%%
% Define v
v = [-6000 -4000 -2000 -1000 -500 0 500 1000 2000 4000 6000];

%%
% We can create a version of the graph compatible with the graphics
% in the previous chapters, i.e., x-Axis 10 cm wide, line thickness of axis 
% 0.6 pt, line thickness of curve 0.5 pt, axis labels 8 pt sans serif,
% title 10 pt sans serif. NOTE: set current figure property COLORMAP to
% DEMCMAP which is a colormap contained in the Mapping Toolbox appropriate
% to terrain elevation data.
figure1 = figure('Color',[1 1 1]);
axes1 = axes('Units','Centimeters',...
    'Box','on',...
    'LineWidth',0.6,...
    'FontSize',8);
contourf1 = contourf(LON,LAT,ETOPO2,v);
title1 = title('ETOPO2 Data Set',...
    'FontSize',10);
xlabel1 = xlabel('Longitude',...
    'FontSize',8);
ylabel1 = ylabel('Latitude',...
    'FontSize',8);
set(gcf,'Colormap',demcmap(ETOPO2))

hold on
plot(-21.0262,18.0772,...
    'Marker','o','MarkerSize',10,...
    'MarkerFaceColor','r','MarkerEdgeColor','k')
text(-21.0262-4.5,18.0772,'ODP 659',...
    'FontSize',20,'Color','r')

%%
% Change view of a surface plot
figure1 = figure('Color',[1 1 1]);
axes1 = axes('Visible','off',...
    'Units','centimeters',...
    'FontSize',8);
hold(axes1,'all');
surf1 = surf(LON,LAT,ETOPO2,...
    'FaceColor','interp',...
    'EdgeColor','none');
set(gca,'View',[20 60])
set(gcf,'Colormap',demcmap(ETOPO2))

hold on
plot(-21.0262,18.0772,...
    'Marker','o','MarkerSize',10,...
    'MarkerFaceColor','r','MarkerEdgeColor','k')
text(-21.0262-3,18.0772+1,'ODP 659',...
    'FontSize',20,'Color','r')

%%
% Combination with 3D contours using v
figure1 = figure('Color',[1 1 1]);
axes1 = axes('Visible','off',...
    'Units','centimeters',...
    'FontSize',8);
hold(axes1,'all');
surf1 = surf(LON,LAT,ETOPO2,...
    'FaceColor','interp',...
    'EdgeColor','none');
contour31 = contour3(LON,LAT,ETOPO2,v,'k');
set(gca,'View',[20 60])
set(gcf,'Colormap',demcmap(ETOPO2))

hold on
plot(-21.0262,18.0772,...
    'Marker','o','MarkerSize',10,...
    'MarkerFaceColor','r','MarkerEdgeColor','k')
text(-21.0262-3,18.0772+1,'ODP 659',...
    'FontSize',20,'Color','r')

%%
% Create light
figure1 = figure('Color',[1 1 1]);
axes1 = axes('Visible','off',...
    'Units','centimeters',...
    'FontSize',8);
hold(axes1,'all');
surf1 = surf(LON,LAT,ETOPO2,...
    'SpecularExponent',20,...
    'FaceLighting','phong',...
    'FaceColor','interp',...
    'EdgeColor','none');
light1 = light('Parent',axes1,...
    'Style','local',...
    'Position',[145 70 155801]);
set(gca,'View',[20 60])

hold on
plot(-21.0262,18.0772,...
    'Marker','o','MarkerSize',10,...
    'MarkerFaceColor','r','MarkerEdgeColor','k')
text(-21.0262-3,18.0772+1,'ODP 659',...
    'FontSize',20,'Color','r')