Last Change

2011-10-30 20:13

Introduction

The Computer Program in Seismology location program elocate is a simple program that can easily be run in a batch mode. The program requires only a velocity model file, VEL.MOD and the data file elocate.dat. The program is described in Chapter 4 of the distributed PROGRAMS.330/DOC/GSAC.pdf/cps330g.pdf.

There are many references to the "leave one out" jackknife procedure for estimating means and variances. A recent reference is

G. A. Prieto, D. J. Thomas, F. L. Vernon, P. M. Shearer and R. L. Parker, (2007). Confidence intervals for earthquake source parameters, Geophys. J. Int 168, 1227-1234. doi:10.1111/j.1365-246X.2006.03257.x

The essence of the the procedure is to systematically leave one observation out of a data set, and then to run the inversion procedure. Then the inverted results are examined. Define K as the original number of observations. The result of each inversion will be a Vj for j=1,K. The mean value is 

     Vbar = (1/K) SUM Vj
and the variance is 
     Var = (K-1)/K  SUM [Vj - Vbar]2

Installation and execution

To run this procedure, you must have installed the Computer Program in Seismology location program elocate. You must also have access to the programs awk/gawk and the C compiler gcc.

Download the following:

or t.tgz and then gunzip -c t.tgz | tar xvf - to get the DOIT and elocate.dat files.

Make the shell script executable: 

	chmod +x DOIT

Execute the shell script: 
	DOIT
or
	./DOIT

Output

All processing results are placed in a sub-directory RESULT: 

     ls RESULT
     elocate.sum.1	elocate.sum.13	elocate.sum.3  elocate.sum.7  elocate.txt.10  elocate.txt.14  elocate.txt.4  elocate.txt.8
     elocate.sum.10	elocate.sum.14	elocate.sum.4  elocate.sum.8  elocate.txt.11  elocate.txt.15  elocate.txt.5  elocate.txt.9
     elocate.sum.11	elocate.sum.15	elocate.sum.5  elocate.sum.9  elocate.txt.12  elocate.txt.2   elocate.txt.6
     elocate.sum.12	elocate.sum.2	elocate.sum.6  elocate.txt.1  elocate.txt.13  elocate.txt.3   elocate.txt.7

The files starting with elocate.txt are the complete output of the elocate program run in batch mode.
The files starting with elocate.sum are the simple summary of the location: 
     cat RESULT/elocate.sum.10
          0.035      37.9515    -77.9603      5.69 20111026223658.022       1319668618.022
The entries are the RMS error, latitude, longitude, depth, origin time as YYYYMMDDHHmmSS.MSEC and the epoch time in seconds after January 1, 1970.

The listed output of the DOIT script is: 

    37.9501 0.0040 -77.9569 0.0059 5.6960 0.4665 1319668617.987 0.057 20111026223657.987
where the columns are 
    Latitude Dlat   Longitude Dlon Depth Ddep Epoch_time Dtime and  origin time in human form.

The location estimate obtained using all data is in the elocate.txt output: 

 Error Ellipse  X=   0.3721 km  Y= 0.6225 km  Theta = 130.7784 deg

 RMS Error        :               0.036              sec
 Travel_Time_Table:          CUS     
 Latitude         :             37.9500 +-    0.0048 N         0.5303 km
 Longitude        :            -77.9570 +-    0.0057 E         0.4946 km
 Depth            :                5.69 +-      0.88 km
 Epoch Time       :      1319668617.988 +-      0.11 sec
 Event Time       :  20111026223657.988 +-      0.11 sec
 Event (OCAL)     :  2011 10 26 22 36 57 988
 HYPO71 Quality   :                  BB
 Gap              :                 102              deg

The jackknife estimated standard error are similar to those for latitude and longitude, and smaller for depth and origin time for this data set.

Description of the shell script

#!/bin/sh


######
#   preserve the original elocate.dat
######
cp elocate.dat elocate.dat.save

#####
#     create a work subdirectory and populate it
#     if it already exists start anew!
#####
rm -fr TEMP
mkdir TEMP
cp VEL.MOD TEMP

#####
#    create a result directory to preserve the 
#    different run output
#####
rm -fr RESULT
mkdir RESULT

#####
#    initialize to get the VEL.MOD file
#####
elocate -VELMOD

#####
#    define the velocity model and starting depth
#####
MODEL=4
DEPTH=10

#####
#    determine the number of lines in elocate.dat
#####
NLINE=`wc elocate.dat.save | awk '{print $1}' `

#####
#    now loop over the lines
#####
COUNT=0
while [ ${COUNT} -lt ${NLINE} ]
do
        head -${COUNT} elocate.dat.save > TEMP/elocate.dat
        COUNTM1=`echo ${COUNT} ${NLINE} | awk '{print $2 -$1 -1}' `
	tail -${COUNTM1} elocate.dat.save >> TEMP/elocate.dat
	COUNT=` echo ${COUNT} | awk '{print $1 + 1}' `
#####
#       use subshell to do the work
#####
	(cp VEL.MOD TEMP; cd TEMP; elocate -M ${MODEL} -D ${DEPTH} -BATCH > elocate.txt)
#####
#       preserve the results with a unique name
#####
	mv TEMP/elocate.txt RESULT/elocate.txt.${COUNT}
	mv TEMP/elocate.sum RESULT/elocate.sum.${COUNT}


done

#####
#    Use Jackknife procedure to get location and one sigma
#####

#####
#    Reference: Equations (1) - (4)
#     G. A. Prieto, D. J. Thomas, F. L. Vernon, P. M. Shearer 
#     and R. L. Parker, (2007).  Confidence intervals for earthquake
#     source parameters, Geophys. J. Int 168, 1227-1234. 
#      doi:10.1111/j.1365-246X.2006.03257.x


#####
#    scan the elocate.sum.* files
#    to perform the leave one out
#
#    the entries in elocate.sum are
# RMS  LAT LON DEP OT OT_CSS
# 0.007 37.9478 -77.9593 5.64 20111026223658.017 1319668618.017

#####
#    make lists of the variables
#####
rm -f list.lat list.lon list.dep list.otcss
cat RESULT/elocate.sum.* | awk '{print $2}' > list.lat
cat RESULT/elocate.sum.* | awk '{print $3}' > list.lon
cat RESULT/elocate.sum.* | awk '{print $4}' > list.dep
cat RESULT/elocate.sum.* | awk '{print $6}' > list.otcss

#####
#    now get the jackknife estimates for each parameter.
#    this ignores any possible covariance
#####

for P in LAT LON DEP OT
do
case ${P} in
	LAT) IND=1 ;FILE=list.lat ;FMTSD="%9.4f" ;;
	LON) IND=2 ;FILE=list.lon ;FMTSD="%9.4f" ;;
	DEP) IND=3 ;FILE=list.dep ;FMTSD="%9.4f" ;;
	OT)  IND=4 ;FILE=list.otcss ;FMTSD="%20.3f" ;;
esac
#####
#     M is the number of entries in the file
#####
M=`wc ${FILE} | awk '{print $1}' `
#####
#     create and awk scrip to compute the mean
#####
cat > awkscmn << EOF
BEGIN {SUM=0.0}
{SUM+=\$1}
END { printf "${FMTSD}\n", SUM/NR} 
EOF
MEAN=`cat ${FILE} | awk -f awkscmn  `
#####
#    save the mean in an array
#####
MVAL[${IND}]="${MEAN}"

#####
#     now computed the standard error from the square root of the variance
#####
cat > awkscsd << EOF
BEGIN {SUM=0.0}
{
SUM+= (\$1 -(${MEAN}) )*(\$1 -(${MEAN}) )
}
END { printf "${FMTSD}\n", sqrt( ( NR - 1 ) * SUM / NR ) }
EOF
SD=`cat ${FILE} | awk -f awkscsd`
SDVAL[${IND}]="${SD}"

done


rm -f awkscmn awkscsd list.dep list.lat list.lon list.otcss

#####
#     at this point compile a C program uner LINUX to convert
#      epoch time to human time
#####

#####
#     now compile a C program to convert the time in seconds from
#     January 1, 1970 00:00:00.000 to UTC
#####
cat > c.c << EOF
#include <stdio.h>
#include <time.h>

double TIME=${MVAL[4]} ;
time_t TIMES;
struct tm tm;
char ostr[30];

main()
{
	TIMES=(time_t)TIME;
	gmtime_r( &TIMES, &tm);
	strftime(ostr,sizeof(ostr),"%Y%m%d%H%M%S", &tm);
	printf("%s.%3.3d\n",ostr,(int)( 1000*(TIME - (double)TIMES )+0.49));
}
EOF
gcc c.c
TIMEHUMAN=`a.out`

######
#   now format the nice summary
#####
echo ${MVAL[1]} ${SDVAL[1]}  ${MVAL[2]} ${SDVAL[2]}  ${MVAL[3]} ${SDVAL[3]}  ${MVAL[4]} ${SDVAL[4]} ${TIMEHUMAN}

rm -f a.out c.c