click title to see Google Drive directory JParkCodes
The code agu.f (guess where I wrote it!) reads the cvel_YYYY.DDD.HH files and generates some plots of the phase-velocity as fct of period using a collection of algorithms: weighted sum, trimmed weighted sum, median. Version posted can process up to 100 cvel files. You can expand the arrays for more.
input:
cvel_2006.175.02 --> input file (one per event , just "ls cvel_20* >! c_agu"
cvel_2006.225.10
cvel_2006.225.19
cvel_2006.240.22
stop --------------------> end of input files
3 --------------------> make a plot
4 --------------------> make a plot and write to PS file
cvel_eastwest.ps ---------> PS plotfile
Tuesday, July 1, 2008
sacw_corr_many.f -- Using the wavelet-based surface wave code
click title to see Google Drive directory JParkCodes
sacw_corr_many.f does a 2-station estimation of phase velocity.
sacw --> sac-format data, analyzed with wavelets
corr --> cross-correlation is used to estimated phase shits and phase velocity
many --> code is set up to run a series of data-record pairs.
The posted code is designed for 10 sps data, and computes Slepian wavelets in 16 banks logarithmically spaced between 8 and 56 second center-period. Three complex-valued Slepian wavelets are used, in order to estimate phased cross-correlations with 2 and 4 degrees of freedom. The phase of the cross-correlation is scaled with the center-period of the wavelet to compute a phase-delay. This phase delay is computed at evenly spaced points in the seismogram. The station locations and the back-azimuth of the earthquake are used to estimate a phase velocity from the phase delay, but only at the time point where the maximum wavelet spectrum is found. The wavelet spectrum is computed only from the first Slepian wavelet, because experience demonstrates that this affords a more stable determination of the max-amplitude arrival time.
The input is at the terminal, so the command
sacw_corr_many < c_files_process
will be standard analysis method. This command file has the form (for stations PIIR and VOLR)
2003.325.04.12.54.PIIR.DHZ
2003.325.04.12.56.VOLR.DHZ
3 ---> plot to screen, could be "0" to skip plotting
4 ---> plot to screen and write output to PS file (also could skip plotting with "0")
2003.325.ps ---> name of the PS file for phase-velocity graph
2004.024.13.05.37.PIIR.DHZ
2004.024.13.05.40.VOLR.DHZ
3
4
2004.024.ps
stop ----> end of data input
The data files must have a common start time (despite the differing numbers in the filenames in this example) In practice, one can use SAC commands to achieve this easily. It is convenient to drop the minutes and seconds entries in the filenames, saving the year/julianday/hour filename format, because these parameters will match for different stations, and the minutes/seconds wont.. For instance, if you have a file "infiles" of the filename pairs to analyze, then
ls 20*VOLR*DHZ | sed 's/\./ /g' | awk '{print "read " $1 "." $2 "." $3 "*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off"}' >! c_synch
will generate the SAC macro to synchronize and cut to a common start point:
read 2003.325.04*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
read 2004.009.19*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
read 2004.024.13*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
read 2004.038.21*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
and so on . . . .
The code writes a small ASCII file of phase-vel information, cvel_YYYY.DDD.HH, which looks like this
cat cvel_2004.339.10
56.00 1.547 1.858 13.76 236.6 341.4 7.340
49.18 4.596 12.74 13.76 236.6 342.0 2.470
43.20 2.893 1.626 13.76 236.6 319.2 3.925
37.93 4.331 1.237 13.76 236.6 321.0 2.621
33.32 3.411 0.3311 13.76 236.6 322.8 3.328
29.27 3.258 0.9897E-01 13.76 236.6 334.8 3.484
25.70 2.934 0.2715 13.76 236.6 334.8 3.869
22.58 2.948 0.1586 13.76 236.6 347.4 3.850
19.83 2.939 0.2432 13.76 236.6 358.2 3.863
17.42 2.759 1.053 13.76 236.6 372.0 4.114
15.30 3.739 2.711 13.76 236.6 370.2 3.036
13.43 5.761 2.135 13.76 236.6 417.6 1.971
11.80 4.417 0.4976 13.76 236.6 456.0 2.570
10.37 1.556 0.1375 13.76 236.6 559.2 7.295
9.100 2.039 0.4498 13.76 236.6 725.4 5.568
8.000 -27.40 99.12 13.76 236.6 995.4 -0.4143
period cvel d(cvel) gcarc baz timeOcorrelation delay-time
sacw_corr_many.f does a 2-station estimation of phase velocity.
sacw --> sac-format data, analyzed with wavelets
corr --> cross-correlation is used to estimated phase shits and phase velocity
many --> code is set up to run a series of data-record pairs.
The posted code is designed for 10 sps data, and computes Slepian wavelets in 16 banks logarithmically spaced between 8 and 56 second center-period. Three complex-valued Slepian wavelets are used, in order to estimate phased cross-correlations with 2 and 4 degrees of freedom. The phase of the cross-correlation is scaled with the center-period of the wavelet to compute a phase-delay. This phase delay is computed at evenly spaced points in the seismogram. The station locations and the back-azimuth of the earthquake are used to estimate a phase velocity from the phase delay, but only at the time point where the maximum wavelet spectrum is found. The wavelet spectrum is computed only from the first Slepian wavelet, because experience demonstrates that this affords a more stable determination of the max-amplitude arrival time.
The input is at the terminal, so the command
sacw_corr_many < c_files_process
will be standard analysis method. This command file has the form (for stations PIIR and VOLR)
2003.325.04.12.54.PIIR.DHZ
2003.325.04.12.56.VOLR.DHZ
3 ---> plot to screen, could be "0" to skip plotting
4 ---> plot to screen and write output to PS file (also could skip plotting with "0")
2003.325.ps ---> name of the PS file for phase-velocity graph
2004.024.13.05.37.PIIR.DHZ
2004.024.13.05.40.VOLR.DHZ
3
4
2004.024.ps
stop ----> end of data input
The data files must have a common start time (despite the differing numbers in the filenames in this example) In practice, one can use SAC commands to achieve this easily. It is convenient to drop the minutes and seconds entries in the filenames, saving the year/julianday/hour filename format, because these parameters will match for different stations, and the minutes/seconds wont.. For instance, if you have a file "infiles" of the filename pairs to analyze, then
ls 20*VOLR*DHZ | sed 's/\./ /g' | awk '{print "read " $1 "." $2 "." $3 "*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off"}' >! c_synch
will generate the SAC macro to synchronize and cut to a common start point:
read 2003.325.04*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
read 2004.009.19*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
read 2004.024.13*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
read 2004.038.21*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
and so on . . . .
The code writes a small ASCII file of phase-vel information, cvel_YYYY.DDD.HH, which looks like this
cat cvel_2004.339.10
56.00 1.547 1.858 13.76 236.6 341.4 7.340
49.18 4.596 12.74 13.76 236.6 342.0 2.470
43.20 2.893 1.626 13.76 236.6 319.2 3.925
37.93 4.331 1.237 13.76 236.6 321.0 2.621
33.32 3.411 0.3311 13.76 236.6 322.8 3.328
29.27 3.258 0.9897E-01 13.76 236.6 334.8 3.484
25.70 2.934 0.2715 13.76 236.6 334.8 3.869
22.58 2.948 0.1586 13.76 236.6 347.4 3.850
19.83 2.939 0.2432 13.76 236.6 358.2 3.863
17.42 2.759 1.053 13.76 236.6 372.0 4.114
15.30 3.739 2.711 13.76 236.6 370.2 3.036
13.43 5.761 2.135 13.76 236.6 417.6 1.971
11.80 4.417 0.4976 13.76 236.6 456.0 2.570
10.37 1.556 0.1375 13.76 236.6 559.2 7.295
9.100 2.039 0.4498 13.76 236.6 725.4 5.568
8.000 -27.40 99.12 13.76 236.6 995.4 -0.4143
period cvel d(cvel) gcarc baz timeOcorrelation delay-time
Code for surface-wave phase-velocity estimation sacw_corr_many.f
click title to see Google Drive directory JParkCodes
c program sacw_corr_many
c JJP 04/13/07 updated 7/24/07 updated 12/11/07, again 6/25/08
c compute phase delays as function of frequency using the multiwavelet transform
c input two fil,es, typically BHZ components at different stations
c cross-correlate the wavelet transforms
c output is a colormap of time and frequency dependent time delay
c plotted with the wavelet transform of the two time series
c
c f77 -o /park/backup/bin/sacw_corr_many sacw_corr.f /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c xf77 -o /park/backup/bin/sacw_corr_many sacw_corr.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c f77 -o /Users/jjpark/bin/sacw_corr_many sacw_corr.f /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/sacw_corr_many sacw_corr_many.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
c to read SAC-format data, filename.bh? <-- br="" e="" for="" n="" routine="" subs="" z="">c
c generates wavelets with eigenvector decomposition + interpolation
c calculate padded FFT of data for convolution with wavelet FFTs
c general use - calculates set of wavelets for a specified number of points
c interpolates and fft other lengths on the fly
c max=12 real-wavelets max = 6 complex-valued wavelets
c
c code is based on the algorithms described in
c
c Lilly, J., and J. Park, Multiwavelet spectral and polarization analysis
c of seismic records, Geophysical J. International, v122, 1001-1021, 1995.
c
c Park, J. and M. E. Mann, Interannual temperature events and shifts in
c global temperature: A multiwavelet correlation approach, Earth Interactions,
c v. 4, online paper 4-001, 2000.
c
C SEE ALSO:
c Bear, L. K., and G. L. Pavlis, Estimation of slowness vectors and
c their uncertainties using multi-wavelet seismic array processing,
c Bull. Seismol. Soc. Am., v87, 755--769, 1997.
c
c this code generates 2-D arrays of spectra, polarization and statistical stuff
c it writes GMT scripts to display these arrays in dazzling color
c
c for an input filename sacfile.bh? the GMT scripts are written to
c GMD_sacfile.bh
c
c
c the ASCII files to plot are named
c cpt_nilcut.sh spw_nilcut.sh1 spw_nilcut.sh2
c dat_nilcut.sh
c cpd_nilcut.sh del_nilcut.sh
c
c dat_* is data to plot
c cpt_* and cohere.cpt are colorfiles for GMT
c del_* dep_* and dem_* are the time delay arrays
c spw_* wavelet spectra
c
c you can fill a lot of disk space if you are not tidy with old files!
c
c user must make scripts executable ("chmod +x GM*_sacfile.bh"),
c execute the scripts and display the postscript files
c program sacw_corr_many
c JJP 04/13/07 updated 7/24/07 updated 12/11/07, again 6/25/08
c compute phase delays as function of frequency using the multiwavelet transform
c input two fil,es, typically BHZ components at different stations
c cross-correlate the wavelet transforms
c output is a colormap of time and frequency dependent time delay
c plotted with the wavelet transform of the two time series
c
c f77 -o /park/backup/bin/sacw_corr_many sacw_corr.f /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c xf77 -o /park/backup/bin/sacw_corr_many sacw_corr.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c f77 -o /Users/jjpark/bin/sacw_corr_many sacw_corr.f /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/sacw_corr_many sacw_corr_many.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
c to read SAC-format data, filename.bh? <-- br="" e="" for="" n="" routine="" subs="" z="">c
c generates wavelets with eigenvector decomposition + interpolation
c calculate padded FFT of data for convolution with wavelet FFTs
c general use - calculates set of wavelets for a specified number of points
c interpolates and fft other lengths on the fly
c max=12 real-wavelets max = 6 complex-valued wavelets
c
c code is based on the algorithms described in
c
c Lilly, J., and J. Park, Multiwavelet spectral and polarization analysis
c of seismic records, Geophysical J. International, v122, 1001-1021, 1995.
c
c Park, J. and M. E. Mann, Interannual temperature events and shifts in
c global temperature: A multiwavelet correlation approach, Earth Interactions,
c v. 4, online paper 4-001, 2000.
c
C SEE ALSO:
c Bear, L. K., and G. L. Pavlis, Estimation of slowness vectors and
c their uncertainties using multi-wavelet seismic array processing,
c Bull. Seismol. Soc. Am., v87, 755--769, 1997.
c
c this code generates 2-D arrays of spectra, polarization and statistical stuff
c it writes GMT scripts to display these arrays in dazzling color
c
c for an input filename sacfile.bh? the GMT scripts are written to
c GMD_sacfile.bh
c
c
c the ASCII files to plot are named
c cpt_nilcut.sh spw_nilcut.sh1 spw_nilcut.sh2
c dat_nilcut.sh
c cpd_nilcut.sh del_nilcut.sh
c
c dat_* is data to plot
c cpt_* and cohere.cpt are colorfiles for GMT
c del_* dep_* and dem_* are the time delay arrays
c spw_* wavelet spectra
c
c you can fill a lot of disk space if you are not tidy with old files!
c
c user must make scripts executable ("chmod +x GM*_sacfile.bh"),
c execute the scripts and display the postscript files
Subscribe to:
Posts (Atom)