#!/bin/sh
#####
# define the location of the deconvolved, QCd waveforms
# of both the large event, DIRBIG and the smaller event used
# for the empirical Green function, e.g.,
#DIRBIG=/home/rbh/PROGRAMS.310t/MOMENT_TENSOR/MECH.TEL/20070815234057/DAT.TEL
#DIRSMALL=/home/rbh/PROGRAMS.310t/MOMENT_TENSOR/MECH.TEL/20070818025235/DAT.TEL
#####
DIRBIG=path_to_big_event_waveform_directory
DIRSMALL=path_to_small_event_waveform_directory
#####
# ONLY CHANGE BELOW HERE TO FINE TUNE
#####
#####
# define the CUT window with respect to P
#####
CUTL=-40
CUTH=200
#####
# define the frequency band for the pre-filter
#####
FL=0.01
FH=0.20
#####
# define the DELTA for the traces used in the decon
#####
INTDELTA=0.25
#####
# modified 10 AUG 2009 from the script used for Wells and L'Aquila
# Previously:
# ls $DIRBIG $DIRSMALL | sort | uniq > uniq.list
# this will not work since early IRIS data did not use location codes
# So
# 1. Assume thaat only one location code is used
#` 2. build up a list of unique STATION COMPONENT codes
#####
#####
# initialize by cleaning up
#####
rm -f stacmp
rm -f *.decon
if [ -d DECONDIR ]
then
rm -fr DECONDIR
fi
mkdir DECONDIR
for i in ${DIRBIG}/*[ZRT] ${DIRSMALL}/*[ZRT]
do
saclhdr -KSTNM -KCMPNM -NL $i >> stacmp
done
cat stacmp | sort | uniq > uniqstacmp
#####
# do the deconvolution
#####
#####
# select the ALPHA parameter
# The approximate relation between the ALPHA and Gaussian pulse duration is
# ALP Duration Corner Frequency
# 0.5 8 s ~ 0.2 Hz
# 1.0 4 s ~ 0.4 Hz
# 2.5 1.6s ~ 1 Hz
#
# The ALP parameter will affect the resolution of the final fit
# A value of 2.5 may be too extreme for real data
#
# The choice of the ALPHA parameter affects the resolution of the
# resulting empirical source time function
#
# The interative decon program computes a measure of fit, we use the
# LEVEL parameters to guide the acceptance of the results
#####
for ALP in 0.5 1.0
do
case ${ALP} in
0.5) LEVEL=70. ;;
1.0) LEVEL=60. ;;
2.5) LEVEL=50. ;;
esac
#####
# select the component
# later use a case statement to set the cut, e.g., a vs t1
#####
for COMP in HZ
do
grep ${COMP} uniqstacmp > unuq.list
rm -fr BIG.DIR SML.DIR
mkdir BIG.DIR
mkdir SML.DIR
while read KSTNM KCMPNM
do
if [ -f ${DIRBIG}/${KSTNM}*${KCMPNM} ]
then
if [ -f ${DIRSMALL}/${KSTNM}*${KCMPNM} ]
then
cp ${DIRBIG}/${KSTNM}*${KCMPNM} BIG.DIR/${KSTNM}${KCMPNM}
cp ${DIRSMALL}/${KSTNM}*${KCMPNM} SML.DIR/${KSTNM}${KCMPNM}
gsac << EOF
cuterr fillz
cut a ${CUTL} a ${CUTH}
r BIG.DIR/${KSTNM}${KCMPNM} SML.DIR/${KSTNM}${KCMPNM}
hp c ${FL} n 3
lp c ${FH} n 3
interpolate delta ${INTDELTA}
w N D
q
EOF
saciterd -FN N -FD D -2 -ALP ${ALP} -POS -D 20 -N 100
USER5=`saclhdr -USER5 decon.out`
#####
# only use decon that filt 80%
#####
cat > awkprog << EOF
{
ANS="YES"; if(\$1 < ${LEVEL}) ANS="NO" ; print ANS
}
EOF
ANS=`echo $USER5 | awk -f awkprog`
if [ ${ANS} = "YES" ]
then
mv decon.out ${KSTNM}${KCMPNM}.${ALP}.decon
fi
fi
fi
done < unuq.list
#####
# make the map
#####
gsac << EOF
rh BIG.DIR/* SML.DIR/*
map kstnm on r on
quit
EOF
sh map.sh
gm convert -trim map.eps ${COMP}.${ALP}.map.png
done
for i in Z
do
case $i in
Z) KOLOR=2 ;;
T) KOLOR=4 ;;
esac
gsac << EOF
bg plt
r *${i}.${ALP}.decon
prs shd pos color ${KOLOR} tl -20 200 amp 0.3 vl 0 360 az sa 0 ann sta
q
EOF
if [ -f PRS001.PLT ]
then
plotgif -C16 -F7 -W10 -K < PRS001.PLT > ${i}decon.gif
fi
done
mv Zdecon.gif Zdecon.${ALP}.gif
#mv Tdecon.gif Tdecon.${ALP}.gif
mv *.decon DECONDIR
done
#####
# cleanup
#####
rm -fr BIG.DIR
rm -fr SML.DIR
rm -fr N D numerator denominator observed predicted stacmp PRS001.??? awkprog
rm -fr uniqstacmp unuq.list
rm -fr decon.out
|
This is the most important part since this this tells where
to find the files. Each directory has waveform files that end in 'Z'
'R' or 'T' The initial part of this script finds
station-component sets that are in both directories since we need a
pair for the deconvolution
This says that we will consider the time window 40 seconds before and
200 seconds after the predicted P-wave first arrival
If we filter both sets of waveforms identically, the deconvolution
should yield the same result. However, the high-pass removes some low
frequency noise.
To run the deconvolution, both traces must have the same sample rate.
In addition resampling make the process faster since we do not wish to
fit high frequencies. This value may be adequate when the deconvolved
function is much longer than a few seconds. When the function is
shorter, this must be smaller
This result of this loop is a list of the unique station-component
pairs that are common to both waveform directories
Do the deconvolution for these two values of the filter parameter ALP
This is set up for the P-wave. First because we use the vertical
component and because we cut with respect to the P-wave first arrival
time. If you wish to use SH, then the loop should start as
for COMP in HT
and the gsac command cut should be with 't0' instead of 'a'
These two 'if' statements ensure that the station-component is in both
directories, then
copy them to a local work directory
cut both waveforms with respect to p
prefilter both waveforms
resample to a smaller dt
save as N for numerator (big event) and D for denominator (small event)
Do the interative deconvolution and examine the goodness of fit
parameter. To see how 'saciterd' works
saciterd -h
If the fit matches the desired level, save the result with a name that
reflects station, component, ALP
All processing is complete for this value of ALP
Use the 'gsac' command to make a GMT script to plot the stations and
ray paths used for the deconvolution
Use the GraphicsMagick package to convert the map.eps to a PNG
graphics file. If you have ImageMagick installed, the command is
convert -trim map.eps ${COMP}.${ALP}.map.png
All processing is complete
Using 'gsac' plot some aximuthal record sections of the deconvolved
traces.
Convert these directly to a 'gif' image instead of an eps, and then a
PNG
Clean up
|