Code covered by the BSD License  

Highlights from
Seis_Pick

image thumbnail

Seis_Pick

by

 

06 Feb 2012 (Updated )

Seis_Pick provides an interactive picking environment for processing seismic data.

[pickmatrix]=seis_pick(arg1,varargin)
function [pickmatrix]=seis_pick(arg1,varargin)
% [pickmatrix]=seis_pick(sac_string_argument)
%       OR
% [pickmatrix]=seis_pick(traces,dt,nc)
% 
% e.g. [pickmatrix]=seis_pick('stem=wey_1.00 numbers=0:7 components=xyz')
%       OR
%      [pickmatrix]=seis_pick(traces,0.00025,3)
%
%
% seis_pick provides an interactive picking environment geared towards
% processing microseismic events (although suitable for any seismic 
% application)
%
% For the most basic analysis, the user can pick P and S wave onsets.
% For polarisation analysis, the user can pick the end of the P wave and 
% S wave windows, and perform hodogram analysis.
% If P and/or S-wave windows are picked, the user can also analyse the 
% frequency content of each phase, as well as pre-event noise
%
% Written by JP Verdon, University of Bristol, 2011, as part of the BUMPS
% (Bristol University Microseismicity ProjectS) consortium
%
% http://www1.gly.bris.ac.uk/BUMPS/
% http://www1.gly.bris.ac.uk/~JamesVerdon/
%
% INPUTS:  
%         There are two input options: reading traces directly from SAC 
%         files, or using workspace variables
%
% READING FROM SAC FILES:
%         a single string argument is specified, listing a number of 
%         parameters used to identify the sac files.
%
%         stem= filename stem for sac events (e.g. event1.)
%         numbers= the numbers appending the event names (e.g. event1.1)
%         components= the component names appending the files (e.g. xyz)
%
%         As an e.g.: A set of sac files with names event1.1.x, event1.1.y,
%         event1.1.z; event1.2.x (y,z); through to event1.8.x(y,z) would be
%         read using:
%
%         sac_string_argument='event1.1.x would be 'stem=event1.
%         numbers=1:8 components=xyz'
%
% USING WORKSPACE VARIABLES:
%         to use pre-assigned workspace variables, 3 values must be input:
%
%         arg1 = a matrix containing the traces. Each column corresponds to
%         a single component
%
%         arg2 = sampling rate (dt)
%
%         arg3 = number of components per station (usually 3 for
%         microseismic data)
%
%      OPTIONAL: pressing 'L' in the picking window will load a prexisting 
%      pick matrix file, which should be an asci file in the same format as
%      the pick matrices output by seis_pick
%
% OUTPUTS:  pickmatrix = a matrix of P and S picks, and P-wave azimuths
% 	           (pickmatrix = [station_number, P-start, S-start, P-end, 
%                               S-end, Azimuth])
%
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%  Acknowledgements:
%   J. Wookey for the msac_read/write subroutines to read SAC format
%   data
%
%   G. Perron and J.B. Ajo-Franklin for their work on seisplot (now 
%   microseis_plot) for plotting seismic traces
%
%   N.I. Fisher for the circular statistics used by this program
%
%   The sponsors of the BUMPS project for providing funding
%
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%
%  This software is distributed under the term of the BSD free software 
%  license.
%
%  Copyright:
%     (c) 2006-2012, James Verdon
%
%  All rights reserved.
%
%   * Redistribution and use in source and binary forms, with or without
%     modification, are permitted provided that the following conditions are
%     met:
%        
%   * Redistributions of source code must retain the above copyright notice,
%     this list of conditions and the following disclaimer.
%        
%   * Redistributions in binary form must reproduce the above copyright
%     notice, this list of conditions and the following disclaimer in the
%     documentation and/or other materials provided with the distribution.
%     
%   * Neither the name of the copyright holder nor the names of its
%     contributors may be used to endorse or promote products derived from
%     this software without specific prior written permission.
%
%
%   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
%   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
%   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
%   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT
%   OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
%   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
%   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
%   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
%   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
%   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
%   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
%--------------------------------------------------------------------------

ihelp=1; % Help option is set to 1 until unselected

% Find out if input is SAC files or workspace variables
if ischar(arg1) == 1
    isac=1;
    stringarg=arg1;
elseif isfloat(arg1) ==1
    isac=0;
else
    error('Something is wrong with the input arguments')
end

if isac == 1
    % Read SAC events:
    [tr,~,nr,nc]=read_sac_event(stringarg);

    if nc < 1 || nc > 4  % Check that the number of components is ok
        error('number of components must be between 1 - 3')
    end

    dt=tr(1).delta;  % Sampling
    traces=[tr.trace];   % Construct traces

elseif isac == 0
    % read from workspace variables
    % Check that the correct number of arguments is available
    if length(varargin) < 1
        error('dt must be specified in input')
    end
    traces=arg1;
    dt=varargin{1};
    if length(varargin) > 1
        nc=varargin{2};
    else 
        nc = 3; % Default to 3 component data if nc not set
        disp('nc not specified - defaulting to 3 components');
    end
    ntraces=size(traces,2);
    % Check that nc seems reasonable (this is no guarantee)
    if mod(ntraces,nc) ~= 0
        warning('SEIS_PICK:Input_error','nc does not appear to be correct: mod(ntraces,nc) ~= 0')
    end
    nr=ntraces/nc;
end

npts=size(traces,1);  % Number of samples in each trace
tracetime=(npts-1)*dt;  % Length of traces    
ny=0.5/dt;    % Nyquist frequency    
    
% Initialise relevant parameters
bugshift=0;       % Shift to deal with plotting bug
ihod=0;           % Tag to enter hodogram analysis on single trace
ihodall=0;        % Tag to enter hodogram analysis on all traces
ispec=0;          % Tag to enter spectral analysis
ishod(1:nr)=0;    % Tag for S-wave polarisation  
iphod(1:nr)=0;    % Tag for P-wave polarisation
isspec(1:nr)=0;   % Tag for S-wave spectral analysis  
ipspec(1:nr)=0;   % Tag for P-wave spectral analysis
inspec(1:nr)=0;   % Tag for spectral analysis of pre-event noise
ifil(1:nr)=0;     % Tag for any filtering
iflip=0;           % Tag to flip trace order

psang(1:nr)=NaN; % Angle between P and S wave polarisations (should ideally be orthogonall, 90 degrees)
az(1:nr)=NaN;    % Event Azimuth (from P-wave polarisations)
azs(1:nr)=NaN;   % S-wave azimuths
tpb(1:nr)=NaN;   % P-wave arrival times
tsb(1:nr)=NaN;   % S-wave arrival times
tpe(1:nr)=NaN;   % End of P-phase window 
tse(1:nr)=NaN;   % End of S-phase window   

pickmatrix=NaN(nr,6); % Matrix containing picks
pickmatrix(:,1)=1:nr;
uftrace=NaN(size(traces)); % Backup of unfiltered traces

% Read picks from SAC files
if isac == 1
    for i=1:nr
        tpb(i)=tr(nc*(i-1)+1).t0; 
        if tpb(i) < 0
            tpb(i) = NaN;
        end
            tpe(i)=tr(nc*(i-1)+1).t2 ;
        if tpe(i) < 0
            tpe(i) = NaN;
        end
            tsb(i)=tr(nc*(i-1)+1).t1 ;
        if tsb(i) < 0
            tsb(i) = NaN;
        end
            tse(i)=tr(nc*(i-1)+1).t3 ;
        if tse(i) < 0
            tse(i) = NaN;
        end
        pickmatrix(i,:)=[i tpb(i) tsb(i) tpe(i) tse(i) az(i)];
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Begin trace plotting screen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lloc=(3:2:(nr*2+1))/2;   % Dividing lines on plot
cuts=[0 tracetime];      % X limits
zoom=[0.5 lloc(nr)];     % Y limits 
scrsz = get(0,'ScreenSize');   % Initialise figure
figure('Position',[50 50 0.4*scrsz(3) 0.85*scrsz(4)])

iend = 0;
while iend < 1  
    clf
    
    % Plot a box with instructions
    if ihelp == 1
        subplot('Position',[0.8 0.05 0.15 0.9])
        box on
        set(gca,'Xtick',[],'Ytick',[])
        
        text(0.05,0.97,'PICKING CONTROLS:')
        text(0.05,0.95,'p  : Pick P wave start')
        text(0.05,0.91,'P  : Delete P wave start')
        text(0.05,0.89,'s  : Pick S wave start')
        text(0.05,0.87,'S  : Delete S wave start')
        text(0.05,0.85,'o  : Pick P wave end')
        text(0.05,0.83,'O  : Delete P wave end')
        text(0.05,0.81,'a  : Pick S wave end')
        text(0.05,0.79,'A  : Delete S wave end')
        text(0.05,0.77,'w  : Delete ALL picks on one trace')
        text(0.05,0.75,'W  : Delete ALL picks on ALL traces')
        
        text(0.05,0.7,'FILTERING CONTROLS:')
        text(0.05,0.68,'f  : Predictive filter')
        text(0.05,0.66,'b  : Bandpass filter ')
        text(0.05,0.64,'F  : Remove all filtering')
        
        text(0.05,0.59,'ANALYSIS OPTIONS:')
        text(0.05,0.57,'h  : Enter hodogram analysis')
        text(0.05,0.55,'q  : Enter frequency analysis')
        
        text(0.05,0.50,'PLOTTING CONTROLS:')
        text(0.05,0.48,'c  : Cut all traces')
        text(0.05,0.46,'C  : Uncut to view all of the trace')
        text(0.05,0.44,'z  : Zoom in on selected traces')
        text(0.05,0.42,'Z  : Unzoom to see all traces')
        text(0.05,0.40,'V  : Flip trace order') 
        
        text(0.05,0.35,'FILE OPTIONS:')
        text(0.05,0.33,'X  : Exit picking')
        text(0.05,0.31,'E  : Save image of traces as .eps')
        text(0.05,0.29,'L  : Load a pick matrix from file')
        text(0.05,0.27,'m  : Save picks back into sac files')
        text(0.05,0.25,'M  : Save the pick matrix to file')
        text(0.05,0.23,'J  : Toggle help panel')

    end
    
    % Plot the seismic traces
    if ihelp == 1
        subplot('Position',[0.05 0.05 0.7 0.9])
    end
    hold on
    
    % Seisplot_energy function to plot traces and energy envelopes
    microseis_plot(traces,nc,0,tracetime,1,nr,dt,2,pickmatrix);
    
    % Set plot limits
    xlim([cuts(1) cuts(2)]) 
    for j=1:nr
        plot([0,tracetime],[lloc(j),lloc(j)],'k')
    end
    ylim(zoom)
    
    % Write P-wave azimuths on screen
    xlen=(cuts(2)-cuts(1));
    for j=1:nr
        if isfinite(az(j)) == 1
            text(cuts(1)+0.85*xlen,lloc(j)-0.25,['P-wave azimuth = ',num2str(floor(az(j))),'^{\circ}'])
        end
    end
    % Compute average azimuth
    az_for_mean=az(~isnan(az));
    if ~isempty(az_for_mean)
        [azmean,ul,ll] = circ_mean(az_for_mean'.*pi/180);
        azmean=azmean*180/pi;
    	ul=ul*180/pi;
        ll=ll*180/pi;
        text(cuts(1)+0.85*xlen,zoom(2)+0.30,['Average azimuth = ',...
            num2str(floor(azmean)),'^{\circ}, Range = ',num2str(floor(ll)),...
            ' - ',num2str(floor(ul)),'^{\circ}']);
    end
    if iflip == 1
        axis ij
    end
    
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ginput function used to pick traces and perform the appropriate operation 
% based on which button pressed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    [x,y,button]=ginput(1);
    
    % Get the trace pick number
    for j=nr:-1:1
        if y < lloc(j) + bugshift; % + bugshift needed sometimes when picks appear on wrong traces - why????
            nt=j;
        end
    end
    
    % Switch to decide what the button has asked for
    switch button
    
    
    
        %%%%%%%%%% Picks %%%%%%%%%%%%%
        case 112   % 'p' - P wave begins
            tpb(nt)=x;
            
        case 80   % 'P' - delete P wave begins
            tpb(nt)=NaN;
            
        case 111   % 'o' - P wave ends
            tpe(nt)=x;
            
        case 79   % 'O' - delete P wave ends
            tpe(nt)=NaN;
            
        case 115   % 's' - S wave starts
            tsb(nt)=x;
            
        case 83   % 'S' - delete S wave starts
            tsb(nt)=NaN;
            
        case 97    % 'a' - S wave ends
            tse(nt)=x;
            
        case 65    % 'A' - delete S wave ends
            tse(nt)=NaN;
            
        case 119   % 'w' - clear picks for one trace
            tse(nt)=NaN;tpe(nt)=NaN;tpb(nt)=NaN;tsb(nt)=NaN;
            az(nt)=NaN;azs(nt)=NaN;iphod(nt)=0;ishod(nt)=0;
            isspec(nt)=0;ipspec(nt)=0;inspec(nt)=0;
        
        case 87    % 'W' - clear all picks
            tse(1:nr)=NaN;tpe(1:nr)=NaN;tpb(1:nr)=NaN;tsb(1:nr)=NaN;
            az(1:nr)=NaN;azs(1:nr)=NaN;iphod(1:nr)=0;ishod(1:nr)=0;
            isspec(1:nr)=0;ipspec(1:nr)=0;inspec(1:nr)=0;
            
            
            
        %%%%%%%%%% Filtering %%%%%%%%%%%%%
        
        case 102   % 'f' - use predictive filter
            if ifil(nt) == 0
                if isfinite(tpb(nt)) 
                    for j=1:nc
                        uftrace(:,nc*nt+j-nc)=traces(:,nc*nt+j-nc);
                        tmptrace=traces(:,nc*nt+j-nc);
                        noise=(tpb(nt)/dt)-100;
                        filtrace=predictive_filter(tmptrace,noise,dt, [2 2 1 .5]);
                        traces(:,nc*nt+j-nc)=filtrace;
                        clear filtrace tmptrace
                        ifil(nt)=1; 
                    end
                else
                    warning('SEIS_PICK:Input_error','P wave arrival must be picked to define pre-event noise')
                end
            else
                warning('SEIS_PICK:Input_error','Trace has already been filtered')
            end
            
        case 98   % 'b' - bandpass filter
            if ifil(nt) == 0
                freq = input('High and Low pass frequencies: ');
                [b,a] = butter(2,freq(1)/ny,'high');
                for j=1:nc
                    uftrace(:,nc*nt+j-nc)=traces(:,nc*nt+j-nc);
                    filtrace = filter(b,a,traces(:,nc*nt+j-nc));
                    traces(:,nc*nt+j-nc)=filtrace;
                    clear filtrace
                end
                [b,a] = butter(2,freq(2)/ny,'low');
                for j=1:nc
                    filtrace = filter(b,a,traces(:,nc*nt+j-nc));
                    traces(:,nc*nt+j-nc)=filtrace;
                    clear filtrace
                end
                ifil(nt)=1;
            else
                warning('SEIS_PICK:Input_error','Trace has already been filtered')
            end
            
            
        case 70   % 'F' - undo filtering to restore to original data
            if ifil(nt) == 1
                for j=1:nc
                    traces(:,nc*nt+j-nc)=uftrace(:,nc*nt+j-nc);
                end
                ifil(nt)=0;
            end
                  
                  
                  
        %%%%%%%%%% Zooming and re-plotting %%%%%%%%%%%%%
            
        case 99    % 'c' - cut traces
            cuts(1)=x;
            [x,~,button]=ginput(1);
            if button == 99
                cuts(2)=x;
            end
            if cuts(1) > cuts(2)
                cuts=fliplr(cuts);
            end
            
        case 67    % 'C' - uncut trace
            cuts=[0 tracetime];
            
        case 122    % 'z' - zoom on traces
            iz(1) = nt;
            [~,y,button]=ginput(1);
            if button == 122
                for j=nr:-1:1
                    if y < lloc(j) + bugshift; 
                        nt=j;
                    end
                end
            end
            iz(2)=nt;
            if iz(2) < iz(1)
                iz=fliplr(iz);
            end
            zoom=[lloc(iz(1))-1 lloc(iz(2))];
            
        case 90    % 'Z' - unzoom
            zoom=[0.5 lloc(nr)];
            
        case 86    % 'V' - flip trace order
            if iflip == 1
                iflip=0;
            elseif iflip == 0
                iflip=1;
            end


            
        %%%%%%%%%% Other %%%%%%%%%%%%%
            
        case 66    % 'B' - shift to deal with plotting bug
            bugshift=input('Shift for plotting bug: ');
            
        case 104   % 'h' - enter hodogram analysis
            ihod=nt;
            
        case 72    % 'H' - hodogram analysis on all picked traces
            ihodall=1;
            
        case 113   % 'q' - enter spectral analysis
        	ispec=nt;
            
        case 76    % 'L' - load a pick matrix from file
            matrix_fname = input('Load pick matrix file name : ','s');
            data=load(matrix_fname);
            if size(data,1) ~= nr
                warning('SEIS_PICK:Input_error','Size of pick file does not correspond to number of stations');
            end
            tpb=data(:,2)';
            tsb=data(:,3)';
            tpe=data(:,4)';
            tse=data(:,5)';
            az=data(:,6)';
            
            
        case 69    % 'E' - print image in eps format
            set(gcf, 'PaperPositionMode','auto')
            traces_fname=input('Traces plot file name: ','s');%'traces.eps';
            print('-depsc2',traces_fname);
            
        case 77    % 'M' - save pick matrix to ascii file
            matrix_fname = input('Output pick matrix file name : ','s');
            save(matrix_fname,'pickmatrix','-ascii');

        case 74     % 'H' - save pick matrix to ascii file
            if ihelp==0
                ihelp=1;
            elseif ihelp == 1
                ihelp=0;
            end
            
        case 109    % 'm' - save picks back into SAC files
            if isac == 1
                disp('Saving picks to SAC files')
                for k=1:nr
                    if isfinite(tpb(k))
                        for j=1:nc
                            tr(nc*(k-1)+j).t0=tpb(k);
                        end
                    else
                        for j=1:nc
                            tr(nc*(k-1)+j).t0=-12345;
                        end
                    end
                    if isfinite(tsb(k))
                        for j=1:nc
                            tr(nc*(k-1)+j).t1=tsb(k);
                        end
                    else
                        for j=1:nc
                            tr(nc*(k-1)+j).t1=-12345;
                        end
                    end
                    if isfinite(tpe(k))
                        for j=1:nc
                            tr(nc*(k-1)+j).t2=tpe(k);
                        end
                    else
                        for j=1:nc
                            tr(nc*(k-1)+j).t2=-12345;
                        end
                    end
                    if isfinite(tse(k))
                        for j=1:nc
                            tr(nc*(k-1)+j).t3=tse(k);
                        end
                    else
                        for j=1:nc
                            tr(nc*(k-1)+j).t3=-12345;
                        end
                    end
                end
                write_sac_event(stringarg,tr);
            else
                warning('SEIS_PICK:Input_error','SAC files have not been read, and so cannot be saved');
            end
        
        %%%%%%%%%% Ending %%%%%%%%%%%%%
            
        case 88   % 'X' - end picking
        	close
            break
            
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Other things to be done
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Go to hodogram analysis for single trace
    if nc == 1 && ihod > 0
        warning('SEIS_PICK:Input_error','Cannot do polarisations with 1-component data')
    else
        if ihod > 0 
            nt=ihod;
            if isfinite(tpb(nt)) && isfinite(tpe(nt)) 
                iphod(nt)=1;
                ptrace=traces(floor(tpb(nt)/dt):floor(tpe(nt)/dt),(nt-1)*nc+1:(nt-1)*nc+nc);
                pfact=max(sqrt(sum(ptrace.^2,2)));
                ptrace_norm=ptrace/pfact;
                [coeff,score] = princomp(ptrace_norm);
                dirVect = coeff(:,1);
                Np=dirVect;
                az(nt) = atand(dirVect(1)/dirVect(2));
            % Check to ensure az is between 0 - 180 degrees
                iazcheck=0;
                while iazcheck < 1
                    if az(nt) < 0
                        az(nt) = az(nt)+360;
                    elseif az(nt) > 360 
                        az(nt) = az(nt)-360;
                    else
                        break
                    end
                end 
                
                meanX = mean(ptrace_norm,1);
                t = [min(score(:,1))-.2, max(score(:,1))+.2];
                endptsp = [meanX + t(1)*dirVect'; meanX + t(2)*dirVect'];
            else
                iphod(nt)=0;
            end
            if isfinite(tsb(nt)) && isfinite(tse(nt))
                ishod(nt)=1;
                strace=traces(floor(tsb(nt)/dt):floor(tse(nt)/dt),(nt-1)*nc+1:(nt-1)*nc+nc);
                sfact=max(sqrt(sum(strace.^2,2)));
                strace_norm=strace/sfact;
                [coeff,score] = princomp(strace_norm);
                dirVect = coeff(:,1);
                Ns=dirVect;
                azs(nt) = atand(dirVect(1)/dirVect(2));
                meanX = mean(strace_norm,1);
                t = [min(score(:,1))-.2, max(score(:,1))+.2];
                endptss = [meanX + t(1)*dirVect'; meanX + t(2)*dirVect'];
            else 
                ishod(nt)=0;
            end    
            
            % Hodogram plots
            clf
            if nc == 3
                subplot(2,3,1)
            elseif nc ==2
                subplot(1,2,1)
            end
            hold on
            box on
            set(gca,'DataAspectRatio',[1 1 1],'XLim',[-1 1],'YLim',[-1 1])
            xlabel('X');ylabel('Y');title('X - Y')    
            if iphod(nt)==1
                plot(ptrace_norm(:,1),ptrace_norm(:,2),'r')
                plot(endptsp(:,1),endptsp(:,2),'k--');
            end
            if ishod(nt)==1
                plot(strace_norm(:,1),strace_norm(:,2),'g')
                plot(endptss(:,1),endptss(:,2),'b--');
            end
            
            if nc == 3
                subplot(2,3,2)
                hold on
                box on
                set(gca,'DataAspectRatio',[1 1 1],'XLim',[-1 1],'YLim',[-1 1])
                xlabel('X');ylabel('Z');title('X - Vertical')  
                if iphod(nt)==1
                    plot(ptrace_norm(:,1),ptrace_norm(:,3),'r')
                    plot(endptsp(:,1),endptsp(:,3),'k--');
                end
                if ishod(nt)==1
                    plot(strace_norm(:,1),strace_norm(:,3),'g')
                    plot(endptss(:,1),endptss(:,3),'b--');
                end
        
                subplot(2,3,3)
                hold on
                box on
                set(gca,'DataAspectRatio',[1 1 1],'XLim',[-1 1],'YLim',[-1 1])
                xlabel('Y');ylabel('Z');title('Y - Vertical')   
                if iphod(nt)==1
                    plot(ptrace_norm(:,2),ptrace_norm(:,3),'r')
                    plot(endptsp(:,2),endptsp(:,3),'k--');
                end
                if ishod(nt)==1
                    plot(strace_norm(:,2),strace_norm(:,3),'g')
                    plot(endptss(:,2),endptss(:,3),'b--');
                end
            end
            
            if nc == 3
                subplot(2,3,[4 6])
            elseif nc ==2
                subplot(1,2,2)
            end
            hold on
            box on
            set(gca,'XTick',[],'YTick',[])
            text(0.02,0.1,'Press x or H to exit')
            text(0.02,0.3,'Press s to save as hodogram-n.eps')

            if iphod(nt) == 1
                text(0.02,0.9,['P-wave polarisation azimuth = ',num2str(floor(az(nt))),'^{\circ}'])
            else
                text(0.02,0.9,'P-wave window not selected')
            end
            if ishod(nt) == 1
                text(0.02,0.7,['S-wave polarisation azimuth = ',num2str(floor(azs(nt))),'^{\circ}'])
            else
                text(0.02,0.7,'S-wave window not selected')
            end

            if iphod(nt) == 1 && ishod(nt) == 1
                psang(nt)=abs(acosd(dot(Np,Ns)));
                text(0.02,0.5,['Angle between P and S = ',num2str(floor(psang(nt))),' ^{\circ}'])
            end
        
            % Buttons to exit or save hodograms
            iexithodogram=0;
            while iexithodogram < 1
                [~,~,button]=ginput(1);
                if button == 120 || button == 72 % 'x' or 'H' - close figure and end
                    break
                end
                if button == 115 % 's' - save hodogram figure as .eps
                    % Replotting to get rid of help panel
                    
                    figure('Position',[50 50 0.1*nc*scrsz(3) 0.3*scrsz(4)],'Visible','off')
                    set(gcf, 'PaperPositionMode','auto')
                    if nc == 3
                        subplot(1,3,1)
                    end
                    hold on
                    box on
                    set(gca,'DataAspectRatio',[1 1 1],'XLim',[-1 1],'YLim',[-1 1])
                    xlabel('X');ylabel('Y');   
                    if iphod(nt)==1
                        plot(ptrace_norm(:,1),ptrace_norm(:,2),'r')
                        plot(endptsp(:,1),endptsp(:,2),'k--');
                    end
                    if ishod(nt)==1
                        plot(strace_norm(:,1),strace_norm(:,2),'g')
                        plot(endptss(:,1),endptss(:,2),'b--');
                    end
                    
                    if nc == 3
                        subplot(1,3,2)
                        hold on
                        box on
                        set(gca,'DataAspectRatio',[1 1 1],'XLim',[-1 1],'YLim',[-1 1])
                        xlabel('X');ylabel('Z');  
                        if iphod(nt)==1
                            plot(ptrace_norm(:,1),ptrace_norm(:,3),'r')
                            plot(endptsp(:,1),endptsp(:,3),'k--');
                        end
                        if ishod(nt)==1
                            plot(strace_norm(:,1),strace_norm(:,3),'g')
                            plot(endptss(:,1),endptss(:,3),'b--');
                        end
        
                        subplot(1,3,3)
                        hold on
                        box on
                        set(gca,'DataAspectRatio',[1 1 1],'XLim',[-1 1],'YLim',[-1 1])
                        xlabel('Y');ylabel('Z');   
                        if iphod(nt)==1
                            plot(ptrace_norm(:,2),ptrace_norm(:,3),'r')
                            plot(endptsp(:,2),endptsp(:,3),'k--');
                        end
                        if ishod(nt)==1
                            plot(strace_norm(:,2),strace_norm(:,3),'g')
                            plot(endptss(:,2),endptss(:,3),'b--');
                        end
                    end
                    hodogram_fname=['hodogram_',num2str(nt),'.eps'];
                    print('-depsc2',hodogram_fname);
                    close
                
                end
            end
            ihod=0;
            clear ptrace strace ptrace_norm strace_norm pfact sfact ...
                coeff score dirVect meanX t endptss endptsp Np Ns
        end
    end
    
    % Hodogram analysis on all traces
    if ihodall == 1;
        for nt=1:nr
            if isfinite(tpb(nt)) && isfinite(tpe(nt)) 
                ptrace=traces(floor(tpb(nt)/dt):floor(tpe(nt)/dt),(nt-1)*nc+1:(nt-1)*nc+nc);
                pfact=max(sqrt(sum(ptrace.^2,2)));
                ptrace_norm=ptrace/pfact;
                [coeff,~] = princomp(ptrace_norm);
                dirVect = coeff(:,1);
                az(nt) = atand(dirVect(1)/dirVect(2));
                % Check to ensure az is between 0 - 360 degrees
                iazcheck=0;
                while iazcheck < 1
                    if az(nt) < 0
                        az(nt) = az(nt)+360;
                    elseif az(nt) > 360 
                        az(nt) = az(nt)-360;
                    else
                        break
                    end
                end 
            end
        end
        ihodall=0;
        clear ptrace pfact ptrace_norm coeff score dirvect
    end
    
    % Go to spectral analysis
    if ispec> 0
        nt=ispec;    
        
       % Full trace
       % [spec1,f]=plotfft(traces(:,1),dt);
       % [spec2,~]=plotfft(traces(:,2),dt);
       % [spec3,~]=plotfft(traces(:,3),dt);
        
        % Pre-event noise
        if isfinite(tpb(nt))
            inspec(nt)=1;
            noisetrace=zeros(size(traces,1),nc);
            noisetrace(1:floor(tpb(nt)/dt),:)=...
                traces(1:floor(tpb(nt)/dt),(nt-1)*nc+1:(nt-1)*nc+nc);
            for j=1:nc
                [specn(:,j),f]=plotfft(noisetrace(:,j),dt);
            end
        else
            inspec(nt)=0;
        end
        
        % P-wave window
        if isfinite(tpb(nt)) && isfinite(tpe(nt)) 
            ipspec(nt)=1;  
            ptrace=zeros(size(traces,1),nc);
            ptrace(floor(tpb(nt)/dt):floor(tpe(nt)/dt),:)=...
                traces(floor(tpb(nt)/dt):floor(tpe(nt)/dt),(nt-1)*nc+1:(nt-1)*nc+nc);
            for j=1:nc
                [specp(:,j),~]=plotfft(ptrace(:,j),dt);
            end
        else
            ipspec(nt)=0;
        end
        
        % S-wave window
        if isfinite(tsb(nt)) && isfinite(tse(nt)) 
            isspec(nt)=1;
            strace=zeros(size(traces,1),nc);
            strace(floor(tsb(nt)/dt):floor(tse(nt)/dt),:)=...
                traces(floor(tsb(nt)/dt):floor(tse(nt)/dt),(nt-1)*nc+1:(nt-1)*nc+nc);
            for j=1:nc
                [specs(:,j),~]=plotfft(strace(:,j),dt);
            end
        else
            isspec(nt)=0;
        end
        
        
        clf
        for k=1:nc
            subplot(2,nc,k)
            box on
            if inspec(nt) == 1;
                loglog(f,specn(:,k),'b');
                hold on
            end
            if ipspec(nt) == 1;
                loglog(f,specp(:,k),'r');
                hold on
            end
            if isspec(nt) == 1;
                loglog(f,specs(:,k),'g');
                hold on
            end
            title(['Component ',num2str(k)])
            xlabel('Frequency (Hz)')
            ylabel('Power')
        end
        subplot(2,nc,[nc+1 2*nc])
        hold on
        box on
        set(gca,'XTick',[],'YTick',[])
        text(0.02,0.1,'Press x or Q to exit')
        text(0.02,0.3,'Press s to save as spectral-n.eps')
        text(0.02,0.9,'Blue = pre-event noise')
        text(0.02,0.7,'Red = P-wave')
        text(0.02,0.5,'Green = S-wave')
        
        % Buttons to exit or save hodograms
        iexitspectral=0;
        while iexitspectral < 1
        	[~,~,button]=ginput(1);
            if button == 120 || button == 81% 'x' - close figure and end
                break
            end
            if button == 115 % 's' - save spectral figure as .eps
                % Replotting to get rid of help panel
                figure('Position',[50 50 0.2*nc*scrsz(3) 0.3*scrsz(4)],'Visible','off')
                set(gcf, 'PaperPositionMode','auto')
                for k=1:nc
                    subplot(1,nc,k)
                    box on
                    if inspec(nt) == 1;
                        loglog(f,specn(:,k),'b');
                        hold on
                    end
                    if ipspec(nt) == 1;
                        loglog(f,specp(:,k),'r');
                        hold on
                    end
                    if isspec(nt) == 1;
                        loglog(f,specs(:,k),'g');
                        hold on
                    end
                    title(['Component ',num2str(k)])
                    xlabel('Frequency (Hz)')
                    ylabel('Power')
                end
                spectral_fname=['spectral_',num2str(nt),'.eps'];
                print('-depsc2',spectral_fname);
                close
            end
        end

        ispec=0;
        clear ptrace strace ntrace specn specp specs f
    end
    
    
    % Last thing to do - update pick matrix
    for j=1:nr
        pickmatrix(j,:)=[j tpb(j) tsb(j) tpe(j) tse(j) az(j)];
    end
        
end    
   

    

    
    
    
    
    
    
    

Contact us