% crustal_depth
% 
% Script to plot crustal thickness and wave speed ratios for a the
% southwest United States.  The script reads in a data file containing data
% for the entire world, takes the subset of that data in the southwest
% U.S., and plots it in  several different ways.
%
% Author: Robert Wagner
% Date: 2nd March 2009
% Last Modified: 12th March 2009



clear

% Import the data into a set of cell arrays of strings, one per column of 
% data.  Then, convert the needed columns into numeric data.

% Columns, in order:
%Network Code, Station Code, Station Latitude, Station Longitude, Station Elevation,
%Crustal Thickness, Std. Dev. (thickness), Vp/Vs, Std. Dev. (Vp/Vs),
%Assumed Vp, Vs, Poissons Ratio, Number of Earthquakes Used, Residual
%Complexity 
ImportCSV('ears_overview.txt')
stationLats = convert_string_array(Column03);
stationLons = convert_string_array(Column04);
stationEles = convert_string_array(Column05);
crustThickness = convert_string_array(Column06);
vpVs = convert_string_array(Column08);

% subset the data by -120 < longitude < -105 && 31 < latitude < 42
subsetIndices = find(stationLats >= 31 ...
                    & stationLats <= 42 ...
                    & stationLons >= -120 ...
                    & stationLons <= -105);
% Reduce the data set to just the data points in the desired area.
stationLats = stationLats(subsetIndices);
stationLons = stationLons(subsetIndices);
stationEles = stationEles(subsetIndices);
crustThickness = crustThickness(subsetIndices);
vpVs = vpVs(subsetIndices);

% Variables for determining the grid spacing.
gridSize = 0.5; % Degrees
stepNumber = 20; % Blocks to divide area into

% Plot of the data point locations.
set_up_figure(1, 'Locations of siesmic stations in western U.S.', '°East', '°North');
plot_data(stationLons, stationLats);
print('assign7_thick_raw.png','-dpng')

% Triangle mesh of the data points for crustal thickness.
set_up_figure(2, 'Triangles for linear interpolation of thickness', '°East', '°North');
triangle_mesh(stationLons, stationLats, crustThickness)
zlabel('Depth');
print('assign7_thick_tri.png','-dpng')

% Suppress further duplicate point warnings, after the triangle-mesh plot
% throws one.  Let it be known that this data set contins duplicate points.
warning('off', 'MATLAB:griddata:DuplicateDataPoints');
warning('off', 'MATLAB:delaunayn:DuplicateDataPoints');

% Contour map of thickness. Gridded by dividing area into stepNumber
% blocks.  Includes markers for the grid points.
set_up_figure(3, 'Contour map of crustal thickness, with grid markers', '°East', '°North');
grid_contour(stationLons, stationLats, crustThickness, stepNumber, stepNumber, 3, 1)
print('assign7_thick_contour_grid.png','-dpng')

% As above, but without the grid markers.
set_up_figure(4, 'Contour map of crustal thickness', '°East', '°North');
basic_contour(stationLons, stationLats, crustThickness, gridSize, gridSize, 3, 0)
print('assign7_thick_contour.png','-dpng')

% As two above, but for the p-wave/s-wave velocity ratio, not thickness.
set_up_figure(5, 'Contour map of Vp/Vs, with grid markers', '°East', '°North');
grid_contour(stationLons, stationLats, vpVs, gridSize, gridSize, 0.05, 0)
print('assign7_vpvs_contour_grid.png','-dpng')

% As above, but without the grid markers.
set_up_figure(6, 'Contour map of Vp/Vs', '°East', '°North');
basic_contour(stationLons, stationLats, vpVs, gridSize, gridSize, 0.05, 0)
print('assign7_vpvs_contour.png','-dpng')

% Scaled image of the crustal thickness.
set_up_figure(7, 'Color plot of crustal thickness', '°East', '°North');
scaled_image(stationLons, stationLats, crustThickness, gridSize, gridSize, 0)
print('assign7_thick_colormap.png','-dpng')

% Scaled image of the p-wave/s-wave velocity ratio.
set_up_figure(8, 'Color plot of Vp/Vs', '°East', '°North');
scaled_image(stationLons, stationLats, vpVs, gridSize, gridSize, 0)
print('assign7_vpvs_colormap.png','-dpng')

% 3D plot of crustal thickness.
set_up_figure(9, 'Surface map of crustal thickness', '°East', '°North');
plot_surface(stationLons, stationLats, crustThickness, gridSize, gridSize, 0)
zlabel('Depth');
view(-20,62);
print('assign7_thick_3d.png','-dpng')

% 3D plot of p-wave/s-wave velocity ratio.
set_up_figure(10, 'Surface map of Vp/Vs', '°East', '°North');
plot_surface(stationLons, stationLats, vpVs, gridSize, gridSize, 0)
zlabel('Depth');
view(-20,62);
print('assign7_vpvs_3d.png','-dpng')

% Turn duplicate point wrnings back on.
warning('on', 'MATLAB:griddata:DuplicateDataPoints');
warning('on', 'MATLAB:delaunayn:DuplicateDataPoints');