%mean_depth_pcolor
%
%    This script plots the mean depth of earthquakes in a grid over
%    Yellowstone.
%
%    Kevin C. Eagar
%    April 27, 2009
%    Last Updated: 04/27/2009
clear
clc

% Load the text file and assign columns variables
%--------------------------------------------------------------------------
load('ystone_eqs_sort.txt','-ascii');
year = ystone_eqs_sort(:,1);
x_locations = ystone_eqs_sort(:,4);
y_locations = ystone_eqs_sort(:,3);
depth = ystone_eqs_sort(:,5);
magnitudes = ystone_eqs_sort(:,6);

% Set up the range of latitudes and longitudes
%--------------------------------------------------------------------------
stepsize = 0.1;
long_range=min(x_locations):stepsize:max(x_locations);
lat_range = min(y_locations):stepsize:max(y_locations);

% Here we are looping over the range of longitudes
%--------------------------------------------------------------------------
for n = 1:length(long_range)
    
    % Here we are looping over the range of latitudes
    %----------------------------------------------------------------------
    for m = 1:length(lat_range)
        tf = x_locations <= long_range(n)+stepsize./2 & x_locations >= long_range(n)-stepsize./2 & y_locations <= lat_range(m)+stepsize./2 & y_locations >= lat_range(m)-stepsize./2;
        locations = find(tf);

        % Build the "mean_depth" matrix with rows = latitudes, columns =
        % longitudes
        %------------------------------------------------------------------
        mean_depth(m,n) = mean(depth(locations));
    end
end

% Here is a little processing to get the right vertices
%--------------------------------------------------------------------------
long_range_vert = long_range - (stepsize./2);
lat_range_vert = lat_range - (stepsize./2);
long_range_vert = [long_range_vert (long_range_vert(end)+stepsize)];
lat_range_vert = [lat_range_vert (lat_range_vert(end)+stepsize)];
mean_depth(end+1,:) = NaN .* ones(1,size(mean_depth,2));
mean_depth(:,end+1) = NaN .* ones(size(mean_depth,1),1);

% Plot the figure
%--------------------------------------------------------------------------
figure(4)
clf
winsize = get(gcf,'Position');
pcolor(long_range_vert,lat_range_vert,mean_depth)
set(gca,'Ydir','normal')
set(gca,'Units','Normalized','Position',[0 0 1 1])
% hold on
% plot(x_locations,y_locations,'k.')
axis off
print('pcolor2.jpg','-djpeg')

% Print the East, West, North, and South constraints
%--------------------------------------------------------------------------
fprintf(1,'North: %.3f\n',max(lat_range_vert));
fprintf(1,'South: %.3f\n',min(lat_range_vert));
fprintf(1,'East: %.3f\n',max(long_range_vert));
fprintf(1,'West: %.3f\n',min(long_range_vert));