In addition to this page, you should probably have a look to this tutorial page.
Currently, the W-phase package is hosted as a github repository. To check out the W-phase repository:
git clone https://github.com/eost/wphase.git wphase_package
To update your W-phase repository (pull changes)
cd /to/the/wphase/directory/ git pull origin master
For more details on using git for W-phase, you can read this page.
The w-phase package have only been tested on Unix and Linux computers. You will need the following:
To install these dependencies on MacOs, you can refer to this page.
To install the code, we must first setup a few environment variables. If you use csh or tcsh:
setenv RDSEED /path/to/rdseed/executable setenv GF_PATH /path/to/greens/functions/database setenv WPHASE_HOME /path/to/wphase
If you use bash:
export RDSEED=/path/to/rdseed/executable export GF_PATH=/path/to/greens/functions/database export WPHASE_HOME=/path/to/wphase/package
These variables are necessary both at installation time and at run time. In addition to these variables, it is handy to include the wphase “bin” directory into the PATH environment variable:
setenv PATH ${PATH}:$WPHASE_HOME/bin
or in bash:
export PATH=${PATH}:$WPHASE_HOME/bin
All theses variable assignment can be included in your .tcshrc or your .bashrc file.
Once you have downloaded the last version of the W-phase package and that the above environment variables are properly defined. Then you can proceed as follows to compile the package:
cd ${WPHASE_HOME}/src make -B
You should also have a look to the W-phase tutorial page.
Before each inversion, it is necessary to create a “run” directory containing two input ascii files:
The format of these files are described in the file formats section.
Path to SEED file(s) must be correctly specified after the keyword 'SEED:' in the i_master file. If multiple SEED files are used for the same inversion, each of them must be referenced properly in i_master using one 'SEED:' line per file.
Once the i_master and CMTSOLUTION files are created, we can extract waveforms and instrument response parameters and perform a rough screening by epicentral distance. This can be done using:
${WPHASE_HOME}/bin/extract.csh
which will extract Z, N, E, 1, 2 channels in SEED file(s).
Data is extracted from SEED volumes and written in SAC format in the subdirectory 'DATA'. PZ files containing instrument responses are also located in this directory. After running an extract*.csh script, screened and windowed data files are listed in 'scr_dat_fil_list' while 'coeffs_rec_lut' contains the coefficients for the time domain deconvolution of the instrument response.
As SEED volumes are no longer supported, waveforms and instrument response can also be extracted from a gzipped tar archive (hereafter referred to as “tgz”) including both SAC and SAC_PZ (poles and zeros) files. Instead of using the “extract.csh” script to extract data from SEED, we can then extract data from a tgz archive using:
${WPHASE_HOME}/bin/tgz2DATA_org.py path_to_tgz_archive.tgz
where “path_to_tgz_archive.tgz” is the path to the tgz archive including SAC and SAC_PZ files. This script will extract SAC and SAC_PZ files a the directory named “DATA_org”.
Two tgz archive formats are supported: WILBER (default) and NIED.
Using the WILBER format, the tgz archive should include SAC files named with a “.SAC” or “.sac” extension and SAC_PZ file names starting with “SAC_PZ”. If data is ordered using the IRIS Wilber 3 system, this format can be obtained directly by selecting “Output format: SAC binary” (little or big-endian) and “Bundle As: gzipped tar archive” in the final Request Data page (see below).
Once the tgz archive is prepared, extraction from a WILBER formatted archive is done by running:
${WPHASE_HOME}/bin/tgz2DATA_org.py path_to_tgz_archive.tgz wilber
As this is the default format, the “wilber” option can be removed from the command line.
Using the NIED tgz archive format, SAC_PZ files should end with “.zp” and SAC file names can have any format (but should not end with “.zp”). Extraction from a NIED archive can be done using
${WPHASE_HOME}/bin/tgz2DATA_org.py path_to_tgz_archive.tgz nied
—
The next step is to calculate the kernel functions associated with the 6 elements of the seismic moment tensor for the stations listed in 'scr_dat_fil_list' and to convolve them with the moment rate function (MRF) specified in the CMTSOLUTION file (time shift and half duration). We must then deconvolve the instrument response from the data and band pass filter each waveform in the frequency pass band specified in i_master. This is performed using:
${WPHASE_HOME}/bin/prepare_wp.csh
for Z, N, E, 1, 2 channels or
${WPHASE_HOME}/bin/prepare_wp.csh Z
for Z channel only
The scripts listed above will generates i_wpinversion which is a list of sac files to be used in the inversion. This original input file usually includes lots of noisy channels which should be removed further. Output kernel functions will be found generally in ./GF
The inversion can then be performed using the Kernels functions in ./GF and data files listed in i_wpinversion.
For more details on options run
${WPHASE_HOME}/bin/wpinversion -h
By default, this program assumes a zero trace moment tensor (i.e. deviatoric moment tensor). This constraint can be disabled by adding the option “-nont”. Another interesting option is “-dc” which enable the inversion for a pure double couple. The options “-stike”, “-dip”, “-rake” and “-mom” provide the possibility to run double couple inversions at fixed strike, dip, rake and/or scalar moment.
Several output files are created after running the inversion:
Since i_wpinversion includes a lot of noisy traces, the result is usually uncertain. To improve the solution, we must remove noisy data. This can be done as follows:
$WPHASE_HOME/bin/wpinversion -ifil o_wpinversion -ofil o_wpinv.th3.0 -th 3.0
(in this example, the new o_wpinversion file is named o_wpinv.th3.0). This step can be repeated two or more times with different thresholds using the option -th (e.g 1.3, 0.9, etc.):
$WPHASE_HOME/bin/wpinversion -ifil o_wpinv.th3.0 -ofil o_wpinv.th1.3 -th 1.3 $WPHASE_HOME/bin/wpinversion -ifil o_wpinv.th1.3 -ofil o_wpinv.th0.9 -th 0.9
$WPHASE_HOME/bin/wpinversion -ifil o_wpinv.th0.9 -ofil o_wpinv.th1.3 -nr 2.0
At the end of the screening process, it is necessary to clean up the run directory so that o_wpinversion correspond to the optimum dataset which can be used later by grid-search and plot routines :
mv o_wpinversion o_wpinv.noth cp o_wpinv.th0.9 o_wpinversion
Data screening and inversion can all be done at once using the RUNA3_lite.csh script:
${WPHASE_HOME}/bin/RUNA3_lite.csh
This script can be used after data extraction in “DATA_org” using either “tgz2DATA_org.py” or “extract.csh” to extract data from a tgz archive or a SEED volume, respectively (see above). Once the “DATA_org” directory is populated, the “RUNA3_lite.csh” script will perform data screening and calculate a moment tensor solution after median data screening and after rejecting traces associated with large misfit using the threshold 5.0 3.0 0.9 (i.e. -th, see 2. in the “Inversion” section above).
If data is extracted directly from SEED volume(s), data extraction/screening and inversion can be done directly by running:
${WPHASE_HOME}/bin/RUNA3.csh
To do the inversion on “Z” components only (i.e., not using the horizontals), use:
${WPHASE_HOME}/bin/RUNA3_only_Z.csh
These two scripts don’t perform the “ratio screening” (i.e. 3. in the “Inversion” section above). To run the same routines with an additional screening based on the rms amplitude ratio (observed/predicted and predicted/observed), the following script can be used:
${WPHASE_HOME}/bin/RUNA3r.csh
will perform median, rms misfit and rms ratio for channels Z, N, E, 1, 2 channels, Warning: RUNAr.csh is not yet fully tested.
In the grid-search scheme, there is a first global rough exploration which is followed by finer samplings around minimal points. If the optimum is found near a bound, the explored space is extended.
Grid searches can be performed using the script
${WPHASE_HOME}/bin/wp_grid_search.py
. A detailed help can be found by using
${WPHASE_HOME}/bin/wp_grid_search.py -h
. The grid-search computation is parallelized using openMP (if available, it takes advantage of muliple cores).
This script run a time-shift (ts) grid search followed by a centroid position grid-search :
All the plotting scripts are coded using python and the module pylab which have to be installed before using them. The module cartopy is also needed for plotting maps but it is optional (even if we recommend to install it).
There are 3 script which can be used to plot the W phase inversion results:
${WPHASE_HOME}/bin/make_grids.py
which be used to display grid-searches results. Use
${WPHASE_HOME}/bin/make_grids.py -h
to have more details on available options and arguments. The 2 other scripts plot observed and synthetic seismograms after W phase inversion.
${WPHASE_HOME}/bin/traces.py
which draw complete seismograms and show station location on a map (if cartopy is available). Please use
${WPHASE_HOME}/bin/traces.py -h
for more details on available options.
${WPHASE_HOME}/bin/make_cwp.py
For more details on available options and arguments see
${WPHASE_HOME}/bin/make_cwp.py -h
to see the options.
ETOPO01 is a 1 arc-minute global relief model of the Earth's surface provided by Amante et al. (2009): <http://www.ngdc.noaa.gov/mgg/global/global.html>
If Cartopy is installed, you can optionally draw ETOPO01 topography and bathymetry. To do so, the path to the ETOPO01 netCDF file must be assigned to the environment variable $ETOPOFILE which can be done in your .tcshrc fil (or .bachrc, .cshrc, etc.). The make_grids.py script have been tested only for grid-registered netCDF file of the ETOPO1 bedrock file available at: <http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/bedrock/grid_registered/netcdf/ETOPO1_Bed_g_gmt4.grd.gz>
You can also have a look to the W-phase tutorial page.
example:
PDE 2003 9 25 19 50 6.40 41.8100 143.9100 27.0 6.9 8.1 HOKKAIDO, JP event name: 092503C time shift: 19.8100 half duration: 33.5000 latitude: 42.2100 longitude: 143.8400 depth: 28.2400 Mrr: 7.770000e+27 Mtt: -4.110000e+27 Mpp: -3.660000e+27 Mrt: 1.320000e+28 Mrp: 2.590000e+28 Mtp: -6.620000e+27
The first line of this file is the PDE.
The second line, 'event name' correspond to the event id
Lines 3 and 4 are the parameters of the STF.
Lines 5, 6 and 7 are the 'latitude', 'longitude' and 'depth' of the centroïd wich will be used for the rotation of horizontal components and for the computation of synthetic kernel functions.
The last 6 lines are optional, they correspond to a moment tensor solution (e.g. a reference GCMT for comparision with the W phase inversion).
This file is composed of several fields (some of them are optional):
example:
EVNAME: Tokachi-Oki_2003 SEED: ../../WP6/SEEDS/2003_tokachi_oki_LH.SEED DMAX: 90. DMIN: 0. CMTFILE: CMTSOLUTION # data deconvolution filt_order: 4 filt_cf1 : 0.001 filt_cf2 : 0.005 filt_pass : 1 IDEC_2: 2 280 0.1 IDEC_3: 0.001 0.1 100 0.03 #inversion WP_WIN: 15. GFDIR: ./GF
grid_search_ts_out corresponds to the time-shift grid-search.
Example:
60.0000 0.20140377 81.7621 0.32070120 000 000 1.0000 81.7621 -35.8500 -72.7100 44.8000 0.60186830 1.83279165 001 000 5.0000 81.7621 -35.8500 -72.7100 44.8000 0.58186860 1.60449952 002 000 9.0000 81.7621 -35.8500 -72.7100 44.8000 0.55865993 1.40553383 003 000 13.0000 81.7621 -35.8500 -72.7100 44.8000 0.53241366 1.23245105 004 000 17.0000 81.7621 -35.8500 -72.7100 44.8000 0.50339138 1.08142045 (…) 014 000 57.0000 81.7621 -35.8500 -72.7100 44.8000 0.20380139 0.31131967 015 000 61.0000 81.7621 -35.8500 -72.7100 44.8000 0.20198340 0.30827735 016 000 65.0000 81.7621 -35.8500 -72.7100 44.8000 0.21088280 0.32324637 017 000 69.0000 81.7621 -35.8500 -72.7100 44.8000 0.22888998 0.35415847 018 000 73.0000 81.7621 -35.8500 -72.7100 44.8000 0.25358477 0.39808714 (…) 039 000 157.0000 81.7621 -35.8500 -72.7100 44.8000 0.65903444 3.48513694 040 000 161.0000 81.7621 -35.8500 -72.7100 44.8000 0.65708573 3.35654281 041 000 165.0000 81.7621 -35.8500 -72.7100 44.8000 0.65410621 3.18301823 042 001 59.0000 81.7621 -35.8500 -72.7100 44.8000 0.20151331 0.30749196 043 001 63.0000 81.7621 -35.8500 -72.7100 44.8000 0.20516181 0.31360143 044 001 67.0000 81.7621 -35.8500 -72.7100 44.8000 0.21889318 0.33689015 045 002 58.0000 81.7621 -35.8500 -72.7100 44.8000 0.20231439 0.30883067 046 002 60.0000 81.7621 -35.8500 -72.7100 44.8000 0.20140377 0.30730901 047 002 62.0000 81.7621 -35.8500 -72.7100 44.8000 0.20324196 0.31038266 048 002 64.0000 81.7621 -35.8500 -72.7100 44.8000 0.20771868 0.31790207
The two first lines correspond respectively to the optimum and to the initial centroid time-shifts, the first column being the time-shift itself and the second column being the associated rms misfit.
The following lines provide some details on the grid-search : 1st col. : index for this time-shift value, 2nd col. : iteration number, 3rd col. : time-shift 4rd col. : centroid latitude 5th col. : centroid longitude 6th col. : centroid depth 7th col. : rms 8th col. : normalized rms
In the example above you can see that a global grid-search is performed from 1 to 165 sec (until computation 41) at the 0th iteration with a time step of 4sec. The 1st iteration corresponds to computations 42-44 which extend the grid-search to unexplored time-shift between 57sec and 69sec with a sampling interval of 2sec. The 2nd iteration perform an even finer sampling to guarantee a 1sec grid-search between 57sec and 65sec. Once this file is available to us (after running the time-shift grid-search), we can run
${WPHASE_HOME}/bin/make_grids.py –t
to create the file grid_search_ts.pdf containing Fig. 1 which corresponds to the above example.
The output centroid position grid-search file grid_search_xy_out as a quite similar same structure replacing the time-shit in the 2 first lines by optimum and initial latitudes, longitudes and depth.
-35.4500 -72.8325 25.5000 0.19118874 -35.2500 -72.7100 25.5000 0.19314386 000 000 60.0000 60.0000 -36.4500 -74.1794 25.5000 0.24214954 0.37750789 001 000 60.0000 60.0000 -36.4500 -73.6896 25.5000 0.22150727 0.34137918 002 000 60.0000 60.0000 -36.4500 -73.1998 25.5000 0.20980594 0.32142471 003 000 60.0000 60.0000 -36.4500 -72.7100 25.5000 0.20938393 0.32071163 (…) 067 001 60.0000 60.0000 -35.6500 -72.4651 25.5000 0.19798867 0.30161961 068 001 60.0000 60.0000 -35.8500 -72.4651 25.5000 0.19830051 0.30213800 069 001 60.0000 60.0000 -35.8500 -72.7100 25.5000 0.19640545 0.29899114 070 001 60.0000 60.0000 -34.6500 -72.9549 25.5000 0.19787217 0.30142600 (…) 104 002 60.0000 60.0000 -35.3500 -72.5875 25.5000 0.19326491 0.29379393 105 002 60.0000 60.0000 -35.3500 -72.7100 25.5000 0.19296721 0.29330241 106 002 60.0000 60.0000 -35.5500 -72.7100 25.5000 0.19360157 0.29435001
The first line corresponds to the optimum centroid location and the second line indicates the a priori location specified in the CMTSOLUTION file. In these two lines, the 1st, 2nd and 3rd column correspond respectively to the centroid latitute, longitude and depth. The 4th column presents the associated rms misfits.
The following lines provide some details on the grid-search : 1st col. : index of the explorated centroid position, 2nd col. : iteration number, 3rd col. : time-shift 4th col. : half duration 5th col. : centroid latitude 6th col. : centroid longitude 7th col. : centroid depth 8th col. : rms 9th col. : normalized rms
Once the grid-search is performed, the grid_search_xy_out is available to us (after running the time-shift grid-search), we can run
${WPHASE_HOME}/bin/make_grids.py –p -b
to create the file grid_search_xy.pdf containing Fig. 1 which corresponds to the above example. The option –b is optional: it allows to use the cartopy module in order to plot coastlines and topography.