%Lecture 7 Stepwise scripts
%
%This script uses loops to do repetitive tasks.  In this case it is taking
%the Yellowstone Earthquake information .txt file and plotting loops and
%nested loops, plotting earthquake hypocentrer locations in 3D and making
%scaled images if hypocenter locations with text to show the depth.  This
%script performs these funations at a stepsize of 0.25, 0.05 and 0.1.
%
%By Wendy Bohon, Modified from Lecture 7 of Computers in Geology glg 598.
%Last modified 2/18/2009

%loads in the yellowstone earthquake txt file and assigns axis properties
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); 

%makes figure 1 w/ x-axis as longitude and y-axis as latitude. clf clears
%the figure
figure(1)
clf
plot(x_locations,y_locations,'k.')
xlabel('longitude')
ylabel('latitude')
axis equal
hold on


%makes a nested loop with stepsize 0.25.  Puts blue circles across the map
%at distances designated by the stepsize
stepsize = 0.25;
long_range=min(x_locations):stepsize:max(x_locations);
lat_range = min(y_locations):stepsize:max(y_locations);
for long = long_range
    for lat = lat_range
        plot(long,lat,'bo')
    end
end

%Works within the nested loop to apply depth (as text) to the earthquake
%hypocenters.  The program searches for the depth values within each step 
%(+/- half the step size from the lat long location of the estimate). Then
%converts the mean and some text to a string (sprintf), and plots that
%amount using text.
stepsize = 0.25;
long_range=min(x_locations):stepsize:max(x_locations);
lat_range = min(y_locations):stepsize:max(y_locations);
for long = long_range
    for lat = lat_range
       tf = x_locations <= long+stepsize./2 & x_locations >= long-stepsize./2 & y_locations <= lat+stepsize./2 & y_locations >= lat-stepsize./2;
       locations = find(tf);
       mean_depth = mean(depth(locations));
       s = sprintf('%2.1f',mean_depth);
       text(long, lat, s);
    end
end

%plots earthquake hypocenters in 3 dimensions.
%The grid on command is used to display the x,y,z grid in the background to
%give some perspective to the plot. 
figure(2)
plot3(x_locations,y_locations,-depth,'r.')
grid on
title('Yellowstone Earthquake Hypocenters')
xlabel('longitude')
ylabel('latitude')
zlabel('depth')

%makes a scaled images with a stepsize of 0.25 (makes a gridded map with
%color).  It gather these values into a matrix that mimics our grid.
stepsize = 0.25;
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)

        % Notice that we have to index "long_range" by n and "lat_range" by m
        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

%This plots a scaled image figure (gridded color).
%Use Ydir to plot the y-axis so that larger values are
%above smaller values. The first argument in set is gca, 
%which tells the set command:"I want to set a property in the current
%axis."
figure(3)
clf
imagesc(long_range,lat_range,mean_depth)
xlabel('longitude')
ylabel('latitude')
set(gca,'Ydir','normal')
colorbar
hold on

%this plots the text of the hypocenter depth on the scaled image
plot(x_locations,y_locations,'k.')
for n = 1:length(long_range)
    for m = 1:length(lat_range)
        s = sprintf('%2.1f',mean_depth(m,n));
        text(long_range(n), lat_range(m), s); 
    end
end

%This script is the same as the script above but with a different stepsize.

%loads in the yellowstone earthquake txt file and assigns axis properties.
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);

%makes figure 1 w/ x-axis as longitude and y-axis as latitude.
figure(4)
clf
plot(x_locations,y_locations,'k.')
xlabel('longitude')
ylabel('latitude')
axis equal
hold on

%makes a nested loop with stepsize 0.05.  Puts blue circles across the map
stepsize = 0.05;
long_range=min(x_locations):stepsize:max(x_locations);
lat_range = min(y_locations):stepsize:max(y_locations);
for long = long_range
    for lat = lat_range
        plot(long,lat,'bo')
    end
end

%works within the nested loop to apply depth (text) to the earthquake hypocenters
stepsize = 0.05;
long_range=min(x_locations):stepsize:max(x_locations);
lat_range = min(y_locations):stepsize:max(y_locations);
for long = long_range
    for lat = lat_range
       tf = x_locations <= long+stepsize./2 & x_locations >= long-stepsize./2 & y_locations <= lat+stepsize./2 & y_locations >= lat-stepsize./2;
       locations = find(tf);
       mean_depth = mean(depth(locations));
       s = sprintf('%2.1f',mean_depth);
       text(long, lat, s);
    end
end

%plots earthquake hypocenters in 3 dimensions
figure(5)
plot3(x_locations,y_locations,-depth,'r.')
grid on
title('Yellowstone Earthquake Hypocenters')
xlabel('longitude')
ylabel('latitude')
zlabel('depth')

%makes a scaled images with a stepsize of 0.05 (makes a gridded map with color)
stepsize = 0.05;
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)

        % Notice that we have to index "long_range" by n and "lat_range" by m
        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

%This plots a scaled image figure (gridded color)
figure(6)
clf
imagesc(long_range,lat_range,mean_depth)
xlabel('longitude')
ylabel('latitude')
set(gca,'Ydir','normal')
colorbar
hold on

%this plots the text of the hypocenter depth on the scaled image
plot(x_locations,y_locations,'k.')
for n = 1:length(long_range)
    for m = 1:length(lat_range)
        s = sprintf('%2.1f',mean_depth(m,n));
        text(long_range(n), lat_range(m), s); 
    end
end

%This script is the same as the ones above but with a different stepsize

%loads in the yellowstone earthquake txt file and assigns axis properties
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);

%makes figure 1 w/ x-axis as longitude and y-axis as latitude.
figure(7)
clf
plot(x_locations,y_locations,'k.')
xlabel('longitude')
ylabel('latitude')
axis equal
hold on


%makes a nested loop with stepsize 0.1.  Puts blue circles across the map
stepsize = 0.1;
long_range=min(x_locations):stepsize:max(x_locations);
lat_range = min(y_locations):stepsize:max(y_locations);
for long = long_range
    for lat = lat_range
        plot(long,lat,'bo')
    end
end

%works within the nested loop to apply depth (text) to the earthquake hypocenters
stepsize = 0.1;
long_range=min(x_locations):stepsize:max(x_locations);
lat_range = min(y_locations):stepsize:max(y_locations);
for long = long_range
    for lat = lat_range
       tf = x_locations <= long+stepsize./2 & x_locations >= long-stepsize./2 & y_locations <= lat+stepsize./2 & y_locations >= lat-stepsize./2;
       locations = find(tf);
       mean_depth = mean(depth(locations));
       s = sprintf('%2.1f',mean_depth);
       text(long, lat, s);
    end
end

%plots earthquake hypocenters in 3 dimensions
figure(8)
plot3(x_locations,y_locations,-depth,'r.')
grid on
title('Yellowstone Earthquake Hypocenters')
xlabel('longitude')
ylabel('latitude')
zlabel('depth')

%makes a scaled images with a stepsize of 0.1 (makes a gridded map with color)
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)

        % Notice that we have to index "long_range" by n and "lat_range" by m
        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

%This plots a scaled image figure (gridded color)
figure(9)
clf
imagesc(long_range,lat_range,mean_depth)
xlabel('longitude')
ylabel('latitude')
set(gca,'Ydir','normal')
colorbar
hold on

%this plots the text of the hypocenter depth on the scaled image
plot(x_locations,y_locations,'k.')
for n = 1:length(long_range)
    for m = 1:length(lat_range)
        s = sprintf('%2.1f',mean_depth(m,n));
        text(long_range(n), lat_range(m), s); 
    end
end