newLeadData_sacformat.tar

LeadDataLoc.sh0000755002211700040140000001270614052773553012567 0ustar  zjiazwzhan#!/bin/sh

SAC_DISPLAY_COPYRIGHT=0
basepath=`pwd`

evefile="eveinfo.lst"
cd datatel
saclst nzyear nzjday nzhour nzmin nzsec nzmsec evla evlo evdp mag f *.z | head -1 | gawk '{print $2,$3,$4,$5,$6,$7,$8,$9,$10,$11;}'  > ../${evefile}
cd -

# --> set directories
if [ -d RESP ]; then
rm RESP/*
echo "found RESP/"
else
mkdir RESP
fi
if [ -d data ]; then
echo "found Vel/";
else
mkdir data;
fi
if [ -d Vel ]; then
echo "found Vel/";
else
mkdir Vel;
fi

# --> rdseed the seed file, modify the sac header
NZYEAR=`cat ${evefile} | gawk '{print $1;}'`;
NZJDAY=`cat ${evefile} | gawk '{print $2;}'`;
NZHOUR=`cat ${evefile} | gawk '{print $3;}'`;
NZMIN=`cat ${evefile} | gawk '{print $4;}'`;
NZSEC=`cat ${evefile} | gawk '{print $5;}'`
NZMSEC=`cat ${evefile} | gawk '{print $6;}'`
EVLA=`cat ${evefile} | gawk '{print $7;}'`
EVLO=`cat ${evefile} | gawk '{print $8;}'`
EVDP=`cat ${evefile} | gawk '{print $9;}'`
MAG=`cat ${evefile} | gawk '{print $10;}'`

################     remove the instrumental response     ###############################
find . -name "SACPZ*" > 1PZs
cat 1PZs | gawk '{split($1,mm,"[_/.]"); if (mm[7]=="--") mm[7]="";print "mv "$1" RESP/"mm[5]"."mm[6]"."mm[7]"."mm[8]}'  | sh
sed -i '6,8s/Z//g' RESP/*
ls -1 *.BH? > BH1
saclst o b f *.BH? > otimelist.dat
cat otimelist.dat | gawk '{print "r "$1;print "ch o 0 b "$3-$2; print "wh"} END{print "q"}' | sac
cat BH1 |gawk '{
        print "r",$1;
        print "rmean";
        print "rtr";
        print "taper";
        print "transfer from polezero subtype RESP/"$1,"to vel freq 0.001 0.005 10 20";print "w Vel/"$1
}
END{print "q"}' | sac

cd Vel/
echo -e "\033[35m Remove the instrumental response \033[0m"

# --> interpolate or decimate/stretch --> delta--0.2
saclst delta f *.*.*.* | gawk '{
	if($2==0.1){
		print "r",$1;
		print "decimate 2";
		print "w over";
	}
	if($2==0.025){
		print "r",$1;
		print "decimate 2";
		print "decimate 2";
		print "decimate 2";
		print "w over"
	}
	if($2==0.01){
		print "r",$1;
		print "decimate 5";
		print "decimate 2";
		print "decimate 2";
		print "w over";
	}
	if($2==0.005){ 
		print "r",$1;
		print"decimate 2";
		print "decimate 5";
		print "decimate 2";
		print "decimate 2";
		print "w over";
	}
        if($2==0.05){
                print "r",$1;
		print "decimate 2";
		print "decimate 2";
                print "w over";
        }
	if($2==0.02){
		print "r",$1;
		print"decimate 2";
		print "decimate 5";
		print "w over";
	}
} END{print "quit"}' | sac
echo -e "\033[35m END interpolate or decimate to 0.05 \033[0m"

# --> Remove stations that do not have 3 components
ls -1 *.BH* | gawk '
BEGIN{
        FS=".";
        num=0;
        a="AA";b="AAA";c="AA";
        }
{
        if ($1==a&&$2==b&&$3==c)
                {num+=1;}
        else if(num<=2&&(a!="AA"||b!="AAA"))
                {       print "rm "a"."b"."c".*";
                        num=1;a=$1;b=$2;c=$3}
                else {  num=1;a=$1;b=$2;c=$3}
}
END{
        if (num<=2) print "rm "a"."b"."c".*";
}' | sh

# --> change those stations which CMPAZs have minor error
ls -1 *BH1 | gawk 'BEGIN{FS="."} {print $1"."$2"."$3".BH"}' > BH1list
ls -1 *BH2 | gawk 'BEGIN{FS="."} {print $1"."$2"."$3".BH"}' > BH2list
comm -12 BH1list BH2list > BH12list.dat
ls -1 *BHE | gawk 'BEGIN{FS="."} {print $1"."$2"."$3".BH"}' > BHElist
ls -1 *BHN | gawk 'BEGIN{FS="."} {print $1"."$2"."$3".BH"}' > BHNlist
comm -12 BHElist BHNlist > BHENlist.dat
cat BH12list.dat | gawk '{print "saclst cmpaz f "$1"1"}' | sh > 1az.dat
cat BH12list.dat | gawk '{print "saclst cmpaz f "$1"2"}' | sh > 2az.dat
cat BHENlist.dat | gawk '{print "saclst cmpaz f "$1"E"}' | sh > eaz.dat
cat BHENlist.dat | gawk '{print "saclst cmpaz f "$1"N"}' | sh > naz.dat
paste 1az.dat 2az.dat > 12az.dat
paste eaz.dat naz.dat > enaz.dat
cat 12az.dat | gawk '{
        diff1=($2+360-$4)%360;
        pota4=($2+270)%360;
        diff2=($4+360-$2)%360;
        potb4=($2+90)%360;
        if (diff1<diff2) {
                if (diff1!=90) print $3,pota4;
        }
        else {
                if (diff2!=90) print $3,potb4;
        }
}' | awk '{
        print "r "$1;
        print "ch cmpaz "$2;
        print "wh";
}
END{ print "q";}' | sac
cat enaz.dat | gawk '{
        diff1=($2+360-$4)%360;
        pota4=($2+270)%360;
        diff2=($4+360-$2)%360;
        potb4=($2+90)%360;
        if (diff1<diff2) {
                if (diff1!=90) print $3,pota4;
        }
        else {
                if (diff2!=90) print $3,potb4;
        }
}' | awk '{
        print "r "$1;
        print "ch cmpaz "$2;
        print "wh";
}
END{ print "q";}' | sac
rm *az.dat
echo -e "\033[35m END change CMPAZs \033[0m"

# --> Rot to gcp
mkdir rtr_mul_100/
ls *BHE | gawk '{split($1,aa,"BH");print"mv",$1,aa[1]"BH1";}' | sh
ls *BHN | gawk '{split($1,aa,"BH");print"mv",$1,aa[1]"BH2";}' | sh
saclst f *.BHZ | gawk '{ print "cut o 0 400"; a=$1;print "r "$1 ;print "rtr" ;print "mul 100"; gsub("BHZ","z",a);print "ch b 0";print "w rtr_mul_100/"a ;print "cut o 0 400"; a=$1; b=$1; gsub("BHZ", "BH1",a); gsub("BHZ","BH2", b); print "r",a,b; print "rtr" ; print "mul 100"; print "rot to gcp"; gsub("BH1","r",a); gsub("BH2","t",b); print "ch b 0"; print "w rtr_mul_100/"a, "rtr_mul_100/"b;} END {print "quit";}' | sac
echo -e "\033[35m Rotate to GCP && CUT \033[0m"

# --> Move repeating data away and start doing gCAP
mv rtr_mul_100/*.[rtz] ../data
cd ../data
mkdir others
mv *.[1-9]0.* others > /dev/null
cd ../

mv Vel/ RESP/ data/
mv data dataloc
rm -rf loc.seed II.*.???  *.HH?

# --> calculate the local Green's functions by FK
cp cmds/fk.sh .
sh fk.sh
rm -rf fk.sh
LeadDataTel.sh0000755002211700040140000001460014052773105012562 0ustar  zjiazwzhan#!/bin/sh

SAC_DISPLAY_COPYRIGHT=0

ls -1 *.SAC > 1SAC1
cat 1SAC1 | gawk '{split($1,m,".");print "mv "$1, m[1]"."m[2]"."m[3]"."m[4]}' | sh

###############    make directories    #################################
if [ -d RESP ]; then
rm RESP/*
echo "found RESP/"
else
mkdir RESP
fi
if [ -d data ]; then
echo "found Vel/";
else
mkdir data;
fi
if [ -d Vel ]; then
echo "found Vel/";
else
mkdir Vel;
fi

################     remove the instrumental response     ###############################
find . -name "SACPZ*" > 1PZs
cat 1PZs | gawk '{split($1,mm,"[_/.]"); if (mm[7]=="--") mm[7]="";print "mv "$1" RESP/"mm[5]"."mm[6]"."mm[7]"."mm[8]}'  | sh
sed -i '6,8s/Z//g' RESP/*
ls -1 *.BH? > BH1
saclst o b f *.BH? > otimelist.dat
cat otimelist.dat | gawk '{print "r "$1;print "ch o 0 b "$3-$2; print "wh"} END{print "q"}' | sac
cat BH1 |gawk '{
        print "r",$1;
        print "rmean";
        print "rtr";
        print "taper";
        print "transfer from polezero subtype RESP/"$1,"to vel freq 0.001 0.005 10 20";print "w Vel/"$1
}
END{print "q"}' | sac

cd Vel/

rm [12]PZs BH[12] 1SAC[12] [12]BH *.BH[12ENZ]
cd Vel/
echo -e "\033[35m Remove the instrumental response \033[0m"

# --> interpolate or decimate/stretch --> delta--0.2
saclst delta f *.*.*.* | gawk '{
	if($2==0.1){
		print "r",$1;
		print "decimate 2";
		print "w over";
	}
	if($2==0.025){
		print "r",$1;
		print "decimate 2";
		print "decimate 2";
		print "decimate 2";
		print "w over"
	}
	if($2==0.01){
		print "r",$1;
		print "decimate 5";
		print "decimate 2";
		print "decimate 2";
		print "w over";
	}
	if($2==0.005){ 
		print "r",$1;
		print"decimate 2";
		print "decimate 5";
		print"decimate 2";
		print"decimate 2";
		print "w over";
	}
        if($2==0.05){
                print "r",$1;
		print"decimate 2";
		print"decimate 2";
                print "w over";
        }
	if($2==0.02){
		print "r",$1;
		print"decimate 2";
		print "decimate 5";
		print "w over";
	}
} END{print "quit"}' | sac
echo -e "\033[35m Modify the sampling frequency to 0.2 \033[0m"

# --> Remove stations that do not have 3 components
ls -1 *.BH* | gawk '
BEGIN{
        FS=".";
        num=0;
        a="AA";b="AAA";c="AA";
        }
{
        if ($1==a&&$2==b&&$3==c)
                {num+=1;}
        else if(num<=2&&(a!="AA"||b!="AAA"))
                {       print "rm "a"."b"."c".*";
                        num=1;a=$1;b=$2;c=$3}
                else {  num=1;a=$1;b=$2;c=$3}
}
END{
        if (num<=2) print "rm "a"."b"."c".*";
}' | sh

# add P or S arrival times
#rm ./IM.*.?H?
sh ../cmds/addp.cmd *.BHZ
sh ../cmds/adds.cmd *.BH[EN12]
echo -e "\033[35m Add P or S arrival times \033[0m"

# --> change those stations which CMPAZs have minor error
########   change those stations which CMPAZs have minor error   ##########
ls -1 *BH1 | gawk 'BEGIN{FS="."} {print $1"."$2"."$3".BH"}' > BH1list
ls -1 *BH2 | gawk 'BEGIN{FS="."} {print $1"."$2"."$3".BH"}' > BH2list
comm -12 BH1list BH2list > BH12list.dat
ls -1 *BHE | gawk 'BEGIN{FS="."} {print $1"."$2"."$3".BH"}' > BHElist
ls -1 *BHN | gawk 'BEGIN{FS="."} {print $1"."$2"."$3".BH"}' > BHNlist
comm -12 BHElist BHNlist > BHENlist.dat
cat BH12list.dat | gawk '{print "saclst cmpaz f "$1"1"}' | sh > 1az.dat
cat BH12list.dat | gawk '{print "saclst cmpaz f "$1"2"}' | sh > 2az.dat
cat BHENlist.dat | gawk '{print "saclst cmpaz f "$1"E"}' | sh > eaz.dat
cat BHENlist.dat | gawk '{print "saclst cmpaz f "$1"N"}' | sh > naz.dat
paste 1az.dat 2az.dat > 12az.dat
paste eaz.dat naz.dat > enaz.dat
cat 12az.dat | gawk '{
        diff1=($2+360-$4)%360;
        pota4=($2+270)%360;
        diff2=($4+360-$2)%360;
        potb4=($2+90)%360;
        if (diff1<diff2) {
                if (diff1!=90) print $3,pota4;
        }
        else {
                if (diff2!=90) print $3,potb4;
        }
}' | awk '{
        print "r "$1;
        print "ch cmpaz "$2;
        print "wh";
}
END{ print "q";}' | sac
cat enaz.dat | gawk '{
        diff1=($2+360-$4)%360;
        pota4=($2+270)%360;
        diff2=($4+360-$2)%360;
        potb4=($2+90)%360;
        if (diff1<diff2) {
                if (diff1!=90) print $3,pota4;
        }
        else {
                if (diff2!=90) print $3,potb4;
        }
}' | awk '{
        print "r "$1;
        print "ch cmpaz "$2;
        print "wh";
}
END{ print "q";}' | sac
rm *az.dat

# --> Rot to gcp
mkdir rtr_mul_100/
ls *BHE | gawk '{split($1,aa,"BH");print"mv",$1,aa[1]"BH1";}' | sh
ls *BHN | gawk '{split($1,aa,"BH");print"mv",$1,aa[1]"BH2";}' | sh
saclst t1 b e f *.BHZ | gawk 'BEGIN{time_begin=t2-50;time_end=$2+300;} {if($3<time_begin&&$4>time_end) print "cut "$2-50,$2+300; a=$1;print "r "$1 ;print "rtr" ;print "mul 100"; gsub("BHZ","z",a);print "ch b 0 a -12345 t1 50 t2 250 t3 50 t4 250";print "w rtr_mul_100/"a ;print "cut t1 -50 t1 300"; a=$1; b=$1; gsub("BHZ", "BH1",a); gsub("BHZ","BH2", b);print "r",a,b; print "rtr" ; print "mul 100"; print "rot to gcp"; gsub("BH1","r",a); gsub("BH2","t",b); print "ch b 0 a -12345 t1 50 t2 250 t3 50 t4 250"; print "w rtr_mul_100/"a, "rtr_mul_100/"b;} END {print "quit";}' | tee | sac
echo -e "\033[35m Rotate to GCP && CUT \033[0m"

# --> Move repeating data away and start doing gCAP
cd rtr_mul_100
read -p " SNR? yes|no: " choose
case $choose in
	yes) echo 'You select yes'
# --> SNR: select data
perl ../../cmds/snr/SNR.pl
rm *.[rtz]
mv ./good_data/*.[rtz] ../../data/
rm -rf good_data
cd ../../
cd data/
mkdir others
mv *.[1-9]0.* others > /dev/null
cd ../
	;;

	no) echo 'You select no'
mv *.[rtz] ../../data/
cd ../../
rm -rf RESP
cd data/
mkdir others
mv *.[1-9]0.* others > /dev/null
cd ../
	;;
esac

# --> reserve the specified number of teleseismic stations 
cp cmds/snr/select.pl ./
perl select.pl data/

cd data/
mkdir all_stations/
cp *.[rtz] all_stations/
rm -rf select_stations/ excluding_stations/

mkdir select_stations 
saclst user0 f *.z | sort -r -nk2 | head -100 | gawk '{split($1,aa,".");
                                                                print "mv "aa[1]"."aa[2]"."aa[3]".r" ,"select_stations/";
                                                                print "mv "aa[1]"."aa[2]"."aa[3]".t" ,"select_stations/";
                                                                print "mv "aa[1]"."aa[2]"."aa[3]".z" ,"select_stations/"}' | sh
mkdir excluding_stations
mv *.[rtz] excluding_stations/
cp select_stations/*.[rtz] .
cd -
mv data datatel

# --> calculate the teleseismic Green's functions by TEL4
cp cmds/tel4.sh ./
sh tel4.sh

mv Vel/ RESP/ datatel/
rm -rf tel4.sh select.pl tel.seed