%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));