User Tools

Site Tools


wphase:documentation

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revision Previous revision
Next revision
Previous revision
Next revision Both sides next revision
wphase:documentation [2014/03/18 18:58]
wphase removed
wphase:documentation [2015/08/04 13:55]
wphase
Line 1: Line 1:
-==== Documentation ==== 
  
-=== Preparing directories,​ i_master and CMTSOLUTION ​=== +====== ​W-phase documentation ======
-Before each inversion, it is necessary to create a “run” directory containing two input ascii files:  +
-  * CMTSOLUTION:​ a file containing the PDE, the event centroid location and timing and optionally a “reference” moment tensor solution. +
-  * i_master: a file containing other parameters such as the band-pass filter parameters, ​ minimum and maximum distances, etc.                          ​+
  
-The format of these files are described in +In addition to this page, you should probably have a look to [[wphase:​tutorial|this tutorial page]].
  
-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. +===== Installation =====
  
 +==== Getting the code ====
  
-=== Extracting data from (mini)SEED file(s) === +Currently, ​the W-phase package is hosted as github repositoryTo check out the W-phase repository
-Once the i_master and CMTSOLUTION files are created, we can extract waveforms and instrument response parameters and perform ​rough screening by epicentral distanceThis can be done using the extract.csh script (usually in bin) :  +<​code>​ 
-- extract.csh extract Z, N, E, 1, 2  channels in SEED file(s), +git clone https://​github.com/​eost/​wphase.git wphase_package 
-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.+</​code>​
  
-=== Calculating Synthetics, deconvolution and filtering === +To update your W-phase repository (pull changes) 
-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 the prepare_wp_*.csh routines.+<​code>​ 
 +cd /to/the/​wphase/​directory/​ 
 +git pull origin master 
 +</​code>​
  
-  *  "​prepare_wp.csh" ​Z, N, E, 1, 2 channels, +For more details on using git for W-phase, [[wphase:repository|you can read this page.]]
-  *  "​prepare_wp.csh Z" : 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+==== Dependencies ====
  
 +The w-phase package have only been tested on Unix and Linux computers. You will need the following:
 +  - csh shell
 +  - gcc and gfortran
 +  - python2.7 (or later)
 +  - You have to install numpy, matplotlib, basemap and netCDF4 to run some python scripts which make figures.
 +  - rdseed
 +  - GMT4
  
-=== Inversion ​=== +To install these dependencies on MacOs, [[wphase:​macos|you can refer to this page]]. 
-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 wpinversion -h.+==== Building the code ====
  
-By defaultthis program assumes ​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. +To install the codewe must first setup few environment variablesIf you use csh or tcsh: 
-For a complete and detailed list of options, please run wpinversion -h .+<​code>​ 
 +setenv GMT_BIN ​     /path/to/gmt/bin 
 +setenv RDSEED ​      /​path/​to/​rdseed/​executable 
 +setenv GF_PATH ​     /​path/​to/​greens/​functions/​database 
 +setenv WPHASE_HOME ​ /​path/​to/​wphase/​package 
 +</​code>​ 
 +If you use bash: 
 +<​code>​ 
 +export GMT_BIN=/​path/​to/​gmt/​bin 
 +export RDSEED=/​path/​to/​rdseed/​executable 
 +export GF_PATH=/​path/​to/​greens/​functions/​database 
 +export WPHASE_HOME=/​path/​to/​wphase/​package 
 +</​code>​
  
-Several output files are created after running the inversion:​ +These variables ​are necessary both at installation time and at run timeIn addition to these variablesit is handy to include ​the wphase "​bin"​ directory into the PATH environment variable: 
-- WCMTSOLUTION (ASCII) is the solution obtained after inversion. +<​code>​ 
-- p_wpinversion (postscript) which can be visualized with gv (or gs, evince ...) and printed with lp or lpr   +setenv PATH         ​${PATH}:​$WPHASE_HOME/​bin 
-- o_wpinversion (ASCII) is the output list of stations. The 1st column is the list of sac filesthe 2nd column ​is the azimuth (degrees), the 3rd column is the epicentral distance (degrees), the 4th column is the index of the first sample of each channel in the concatenated data vector, the 5th column is the index of the first sample of the next channel, the 6th column is the rms misfit, the 7th column is number of samples the 8th column is the concatenated sampthe normalized rms misfit, the 8th column is the  +</​code>​ 
-- fort.15* (ASCII) is a comparison between observed and predicted concatenated W-phase traces. The first column is the data vector, the second column is the predicted data vector after inversion and the (optional) third column is the data predicted from the (optional) reference solution in the CMTSOLUTION file.  +or in bash: 
-- fort.15_LH? files (ASCII) where '?'​ is Z, N, E, 1 or 2 are the same concatenated traces but separated by channels. ​ +<​code>​ 
-- *.synth.sac (SAC files) which are the predicted data after inversion.+export PATH=${PATH}:​$WPHASE_HOME/​bin 
 +</​code>​
  
-Since i_wpinversion includes a lot of noisy traces, the result is usually noisy+All theses variable assignment can be included in your .tcshrc or your .bashrc file.
  
-1. A first screening can be performed using the option ​-med : this will reject data with an anomalously small or large interval between maximum ​and minimum amplitudes.+Once you have downloaded ​the last version of the W-phase package ​and that the above environment variables are properly definedThen you can proceed as follows to compile the package: 
 +<​code>​ 
 +cd ${WPHASE_HOME}/​src 
 +make -B 
 +</​code>​
  
-2. After a first run of the inversion, it is possible to perform another inversion throwing out the stations with a relative mismatch exceeding 3.0 (i.e. really bad data). The file o_wpinversion which was created by the first inversion will be used an input. After inversion with the better dataset, a new o_wpinversion file will be created. For example using a threshold of 3.0 : +-----
->​$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 +
  
-3. Finally, it is possible ​to remove the channels showing a large rms ratios (observed/​predicted and predicted/​observed) using the option –nr (e.g. –nr 2.0): +===== How to run W-phase  =====
->​$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 : +You should also have a look to [[wphase:​tutorial|the W-phase tutorial page]].
-> mv o_wpinversion o_wpinv.noth +
-> cp o_wpinv.th0.9 o_wpinversion+
  
 +==== Preparing directories,​ i_master and CMTSOLUTION ====
 +Before each inversion, it is necessary to create a “run” directory containing two input ascii files: ​
 +  * CMTSOLUTION:​ a file containing the PDE, the event centroid location and timing and optionally a “reference” moment tensor solution.
 +  * i_master: a file containing other parameters such as the band-pass filter parameters, ​ minimum and maximum distances, etc.                          ​
  
-=== RUNA scripts === +The format of these files are described in [[wphase:documentation#​Data formats|the file formats ​section]]
-Steps 2 to 4 can be performed by running RUNA3*csh scripts (usually in bin). The following scripts launch the routines ​described in section 2 - 4 and calculate a moment tensor solution after median data screening and after rejecting the worst traces using the threshold 5.0 3.0 0.9 (i.e. -th, see 2. in section 4): +
-- RUNA3.csh : Z, N, E, 1 and 2  channels, +
-- RUNA3_only_Z.csh : only Z channel. +
-These scripts don’t perform ​the “ratio screening” (i.e. 3. in section ​4)To run the same routines with an additional screening based on the ratio (observed/​predicted and predicted/​observed),​ the following scripts can be used: +
-- RUNA3r.csh : Z, N, E, 1, 2 channels,+
  
-Remark: +Path to SEED file(s) must be correctly specified after the keyword 'SEED:' in the i_master fileIf 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
- RUNAr.csh is not yet fully tested.+
  
 +-----
  
-Several optional ​parameters can be used when running RUNA3*csh scripts. These parameters will apply to the wpinversion program. They are described when running wpinversion -h.  +==== Extracting data from SEED ==== 
-e.g.:  RUNA3.csh -wz 1.0 -wn 0.3 -we 0.5  ​will give a weights of 1.0, 0.3 and 0.5 for Z, N, E components respectively.+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: ​ 
 +<​code>​${WPHASE_HOME}/​bin/​extract.csh</​code>​  
 +which will extract ​Z, N, E, 1, 2  channels in SEED file(s).
  
-To perform step 3 and 4 only, using previously ​extracted ​data (from step 2 or after using your own extraction routines), it is possible to use:  +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 scriptscreened 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.
-- RUNA3_lite.csh :  ZN, E, 1, 2 channels.+
  
 +-----
  
-=== Grid searches ​=== +==== Calculating Synthetics, deconvolution and filtering ==== 
-In the grid-search scheme, there is a first global rough exploration which is followed by finer samplings around minimal pointsIf the optimum is found near a bound, ​the explored space is extended.+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: 
 +<​code>​${WPHASE_HOME}/​bin/​prepare_wp.csh</​code>​  
 +for Z, N, E, 1, 2 channels or 
 +<​code>​${WPHASE_HOME}/​bin/​prepare_wp.csh Z</​code>​  
 +for Z channel only
  
-Grid searches can be performed using the script : +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
-- wp_grid_search.py +
-A detailed help can be found by using  '​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 : +----- 
-By default, the ts grid-search only shift the Green'​s functions to find an optimum centroid timing. At the end of the grid-search, a new inversion is performed by fixing the source half duration (hd) to the optimal value of ts.  +==== Inversion ==== 
- Using the option -s, the grid-search is preformed while considering ts=hd which is more time consuming. ​ +The inversion ​can then be performed using the Kernels functions ​in ./GF and data files listed in i_wpinversion
- Grid-search results ​can be found in grid_search_ts_out (see note 3)A postscript file (ts_p_wpinversion) as well as a CMTSOLUTION file (ts_WCMTSOLUTION) corresponding to the optimum ts will also be generated automatically. In greneral all the output files containing the prefix ts_ are associated to the time-shift grid search. +
-- Then, using the optimum value of ts, it perform a 2D centroid position grid-search (lat/lon). The centroid position grid-search generates a text file named grid_search_xy_out with a structure similar to grid_search_ts_out (see note 3). As for time-shift grid-search,​ you will also find the usual output ​files with the prefix xy_. +
-- It is also possible to perform a 3D grid-search using the option –z.  +
-7. Plot routines +
-All the plotting scripts are coded using python and the module pylab which have to be installed before using them. The module basemap is also needed for plotting maps but it is optional (even if we recommend to install it). +
-There is basically 3 script which can be used to plot the W phase inversion results:+
  
-- make_grids.py can be used to display grid-searches results, ​run (use make_grids.py ​-h to see the options).+For more details on options ​run  
 +<​code>​${WPHASE_HOME}/​bin/​wpinversion ​-h</​code>​
  
-The 2 other python scripts plot observed and synthetic seismograms after W phase inversion +By default, this program assumes a zero trace moment tensor (i.edeviatoric 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.
-In order to plot complete seismograms individually and place station on a map : +
-- traces_global.py : draw complete seismograms and show station location on a map (if basemap is available). Please use traces_global.py -h to see the options+
-make_cwp.py : display ​the concatenated W phase traces (use make_cwp.py -h to see the options).+
  
--------+Several output files are created after running the inversion:​ 
 +  * WCMTSOLUTION (ASCII) is the solution obtained after inversion. 
 +  * p_wpinversion (postscript) which can be visualized with gv (or gs, evince ...) and printed with lp or lpr. 
 +  * o_wpinversion (ASCII) is the output list of stations. The 1st column is the list of sac files, the 2nd column is the azimuth (degrees), the 3rd column is the epicentral distance (degrees), the 4th column is the index of the first sample of each channel in the concatenated data vector, the 5th column is the index of the first sample of the next channel, the 8th column is the rms misfit 
 +  * fort.15* (ASCII) is a comparison between observed and predicted concatenated W-phase traces. The first column is the data vector, the second column is the predicted data vector after inversion and the (optional) third column is the data predicted from the (optional) reference solution in the CMTSOLUTION file.  
 +  * fort.15_LH? files (ASCII) where '?'​ is Z, N, E, 1 or 2 are the same concatenated traces but separated by channels.  
 +  * *.synth.sac (SAC files) which are the predicted data after inversion.
  
-==== Example of run ==== +Since i_wpinversion includes ​lot of noisy traces, the result is usually uncertainTo improve the solution, we must remove noisy data. This can be done as follows: 
-This section provide ​typical example ​of W phase inversion run.  +  - A first screening can be performed using the option -med : this will reject data with an anomalously small or large interval between maximum ​and minimum amplitudes. 
-The first step is to create i_master ​and CMTSOLUTION filesLet’s consider ​the Mw8.8 2010 Maule earthquakeWe can use as input the following CMTSOLUTION ​file +  - After a first run of the inversion, it is possible to perform another inversion throwing out the stations with a relative mismatch exceeding 3.0 (i.e. really bad data). The file o_wpinversion which was created by the first inversion will be used an input. After inversion with the better dataset, a new o_wpinversion ​file will be created. For example using a threshold of 3.0 :<​code>​ 
- PDEW2010 ​ 2 27  6 34 15.60 -35.8500  ​-72.7100  44.8 8.9 8.9 NEAR COAST OF CENTRAL CH                 +$WPHASE_HOME/​bin/​wpinversion -ifil o_wpinversion -ofil o_wpinv.th3.-th 3.0  
-event name    201002270634A +</​code>​ (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.)<​code>​ 
-time shift: ​    81.7621 +$WPHASE_HOME/​bin/​wpinversion -ifil o_wpinv.th3.-ofil o_wpinv.th1.3 -th 1. 
-half duration: ​ 81.7621 +$WPHASE_HOME/​bin/​wpinversion ​-ifil o_wpinv.th1.-ofil o_wpinv.th0.9 -th 0. 
-latitude: ​     ​-35.8500 +</​code>​ 
-longitude: ​    -72.7100 +  Finally, it is possible to remove the channels showing a large rms ratios (observed/​predicted and predicted/​observed) using the option –nr (e.g. –nr 2.0)<​code>​ 
-depth: ​         44.8000 +$WPHASE_HOME/​bin/​wpinversion ​-ifil o_wpinv.th0.9 -ofil o_wpinv.th1.3 -nr 2.0 
-Mrr:       1.040000e+29 +</​code>​
-Mtt:      ​-3.920000e+27 +
-Mpp:      ​-1.000000e+29 +
-Mrt:       ​3.040000e+28 +
-Mrp:      ​-1.520000e+29 +
-Mtp     -1.190000e+28+
  
-This file correspond to the Global CMT solution except that we fixed the centroid latitutelongitude and depth to the PDE values and that the time-shift and half duration have been fixed a priori using some initial magnitude estimate.+At the end of the screening processit 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 : 
 +<​code>​ 
 +mv o_wpinversion o_wpinv.noth 
 +cp o_wpinv.th0.9 o_wpinversion 
 +</​code>​
  
-The i_master file is the following ​:+----- 
 +==== RUNA scripts ==== 
 +Data extraction, screening and inversion described above can be performed by running one of the RUNA3*csh scripts (usually in bin). These scripts perform data extraction/​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 section 4): 
 +<​code>​ 
 +${WPHASE_HOME}/​bin/​RUNA3.csh 
 +</​code>​  
 +will perform median and rms misfit screening for Z, N, E, 1 and 2  channels, 
 +<​code>​ 
 +${WPHASE_HOME}/​bin/​RUNA3_only_Z.csh  
 +</​code>​ 
 +will perform median and rms misfit screening for Z channels only.
  
-EVNAME: ​    ​NEAR_COAST_OF_CENTRAL_CH +These two scripts don’t perform the “ratio screening” (i.e3above). To run the same routines with an additional screening based on the ratio (observed/predicted and predicted/observed), the following scripts can be used
-SEED:       ../SEEDS/201002270634A_eqdat.seed +<​code>​ 
-DMIN       5.00 +${WPHASE_HOME}/​bin/​RUNA3r.csh  
-DMAX:       55.00 +</​code>​ 
-CMTFILE: ​   CMTSOLUTION +will perform median, rms misfit and rms ratio for channels Z, N, E, 1channels,
-filt_order: 4 +
-filt_cf1: ​  0.00100 +
-filt_cf2: ​  ​0.00500 +
-filt_pass:  ​1 +
-IDEC_2:  ​ ​280 ​ 0.1 +
-IDEC_3: ​ 0.001  0.1  100  0.03 +
-GFDIR: ​    GF +
-WP_WIN: ​   15.0+
  
-The dataset is contained in the seed file ../​SEEDS/​201002270634A_eqdat.seed (in this case the seed volume was obtained using the NetDC request tool of the IRIS DMC). It contains a large dataset of LHZ channels from 0 to 180degrees of epicentral distances for a distribution of stations corresponding to FDSN, GSN_BROADBAND and STS1 virtual networks. For the purpose of this example, we decided to restrict the our dataset to the stations having epicentral distances within 5.0 to 55.0 degrees. Of course it is possible to consider more stations by increasing DMAX since the green function database allow us to perform the W phase inversion for stations within 90 degrees of epicentral distances. The filter pass band is 1mHz-5mHz and the time window extent ​is from the P wave until 15*∆.+Remark: RUNAr.csh is not yet fully tested.
  
-Once these files are ready, we can proceed with a moment tensor inversion at the centroid location specified in CMTSOLUTION (i.e. the PDE location in this example): 
-> ${WPHASE_HOME}/​bin/​RUNA3.csh ​ 
-This command will extract the data from seed volume, calculate the kernel functions corresponding to each moment tensor element, filter the waveforms, screens the data and perform inversion. The solution is displayed on the terminal and is given in the WCMTSOLUTION file: 
  
- ​PDEW2010 ​ 2 27  6 34 15.60 -35.8500  ​-72.7100  44.8 8.9 8.9 NEAR COAST OF CENTRAL CH                 +Several optional parameters can be used when running RUNA3*csh scriptsThese parameters will apply to the wpinversion programThey are described when running wpinversion ​-h 
-event name:     ​201002270634A +e.g.:  ​RUNA3.csh -wz 1.-wn 0.3 -we 0.5  will give a weights of 1.0, 0.3 and 0.5 for Z, N, E components respectively.
-time shift: ​    ​81.7621 +
-half duration:  ​81.7621 +
-latitude: ​     ​-35.8500 +
-longitude: ​    -72.7100 +
-depth: ​         44.8000 +
-Mrr:       ​9.106371e+28 +
-Mtt:       5.055543e+27 +
-Mpp:      -9.611925e+28 +
-Mrt:      -8.117296e+28 +
-Mrp:      -1.862157e+29 +
-Mtp:      -1.167994e+28+
  
-Other output files are provided, such as p_wpinversion ​and fort.15 as detailed in previous sections. The ASCII beach ball can be drawn using the cmtascii routine:+To perform data screening ​and inversion only, using previously extracted SAC data, it is possible to use 
 +<​code>​ 
 +${WPHASE_HOME}/​bin/​RUNA3_lite.csh 
 +</​code>​ 
 +that is designed for Z, N ,Z ,1, 2 channels
  
->​${WPHASE_HOME}/​bin/​cmtascii WCMTSOLUTION +----- 
-Moment mag.  :  8.83 +==== Grid searches ==== 
-PDE location : Lat35.85S; Lon ​72.71W;​ Dep44.8 km +In the grid-search scheme, there is a first global rough exploration which is followed by finer samplings around minimal pointsIf the optimum is found near a bound, the explored space is extended.
-Centroid loc.: Lat35.85S; Lon ​72.71W;​ Dep44.8 km +
-Origin time  : 2010/02/27 06:34:15.60 +
-Time delay   : 81.8  sec +
-Half duration: 81.8  sec +
-(…) +
-Best Double Couple: M0=2.24E+29 dyn.cm +
-NP1: Strike=351 ; Dip=14 ; Slip= 60  +
-NP2: Strike=202 ; Dip=78 ; Slip= 97 +
  
-             ######### ​             +Grid searches can be performed using the script <​code>​${WPHASE_HOME}/​bin/​wp_grid_search.py</​code>​. A detailed help can be found by using <​code>​${WPHASE_HOME}/​bin/​wp_grid_search.py ​-h</​code>​. The grid-search computation is parallelized using openMP (if available, it takes advantage of muliple cores).
-         ​----------------- ​         +
-      ----------------####​--- ​      +
-    ----------------########​--- ​    +
-   ​----------------##########​--- ​   +
-  ----------------############​--- ​  +
- ​---- ​  ​---------##############​---  +
------ P --------################​--- +
------   ​--------################​--- +
----------------#################​--- +
---------------######## ​  #######​--- +
--------------#########​ T #######​--- +
- ​------------######### ​  #######​--  +
-  ----------###################​-- ​  +
-   ​---------##################​-- ​   +
-    --------################​--- ​    +
-      -----################​-- ​      +
-         ​--#############​-- ​         +
-             ######​---+
  
-This solution ​is computed ​with a fixed location ​(PDE). It is also possible to determine an optimum centroid position and timing by running ​the  ​wp_grid_search.py script:+This script run a time-shift (ts) grid search followed by a centroid position grid-search : 
 +   * By default, the ts grid-search only shift the Green'​s functions to find an optimum centroid timing. At the end of the grid-search,​ a new inversion ​is performed by fixing the source half duration (hd) to the optimal value of ts. \\ Using the option -s, the grid-search is preformed while considering ts=hd which is more time consuming. \\ Grid-search results can be found in grid_search_ts_out (see note 3). A postscript file (ts_p_wpinversion) as well as a CMTSOLUTION file (ts_WCMTSOLUTION) corresponding to the optimum ts will also be generated automatically. In gereneral all the output files containing the prefix ts_are associated to the time-shift grid search. 
 +   * Then, using the optimum value of ts, it perform a 2D centroid position grid-search (lat/lon). The centroid position grid-search generates a text file named grid_search_xy_out ​with a structure similar to grid_search_ts_out ​(see note 3). As for time-shift grid-search,​ you will also find the usual output files with the prefix xy_. 
 +   ​* ​It is also possible to perform a 3D grid-search using the option –z
  
->​${WPHASE_HOME}/​bin/​wp_grid_search.py –z+----- 
 +==== Plot routines ====
  
-Without ​the '​-z'​ flag, this command performs a 2D grid search with a fixed depth. If the option –z is given ,   the depth is also explored. This command will perform the time-shift grid-search as well as the centroid position grid-search. The grid-search is first performed with a large sampling interval ​and refined around ​the parameters corresponding ​to minimal rms misfits. The solution after time-shift grid-search ​is given in ts_WCMTSOLUTION and the final optimum solution after centroid position grid-search ​is given in xy_WCMTSOLUTION:​+All the plotting scripts are coded using python ​and the module pylab which have to be installed before using them. The module basemap ​is also needed for plotting maps but it is optional (even if we recommend to install it).
  
- ​PDEW2010 ​ 2 27  6 34 15.60 -35.8500 ​ -72.7100 ​ 44.8 8.9 8.9 NEAR COAST OF CENTRAL CH                 +There are 3 script which can be used to plot the W phase inversion results:
-event name    201002270634A +
-time shift: ​    ​60.0000 +
-half duration: ​ 60.0000 +
-latitude: ​     -35.4500 +
-longitude: ​    ​-72.8325 +
-depth: ​         25.5000 +
-Mrr:       ​9.393492e+28 +
-Mtt:       ​2.707494e+26 +
-Mpp:      -9.420567e+28 +
-Mrt:       ​1.669174e+28 +
-Mrp:      -1.567680e+29 +
-Mtp:      -1.197536e+28+
  
-Again, this solution can be drawn using cmtascii ​:+   * The first script is <​code>​${WPHASE_HOME}/​bin/​make_grids.py</​code>​ which be used to display grid-searches results. Use <​code>​${WPHASE_HOME}/​bin/​make_grids.py -h</​code>​ to have more details on available options and arguments. The 2 other scripts plot observed and synthetic seismograms after W phase inversion.  
 +   * In order to plot complete seismograms individually and place station on a map<​code>​${WPHASE_HOME}/​bin/​traces_global.py</​code>​ which draw complete seismograms and show station location on a map (if basemap is available). Please use <​code>​${WPHASE_HOME}/​bin/​traces_global.py -h</​code>​ for more details on available options. 
 +   * To plot concatenated waveforms:<​code>​${WPHASE_HOME}/​bin/​make_cwp.py</​code>​ For more details on available options and arguments see <​code>​${WPHASE_HOME}/​bin/​make_cwp.py -h</​code>​ to see the options.
  
->​${WPHASE_HOME}/​bin/​cmtascii WCMTSOLUTION+-------
  
-Moment mag.  :  8.78 +==== Using ETOPO01 Global Relief Model ====
-PDE location : Lat35.85S; Lon ​72.71W;​ Dep44.8 km +
-Centroid loc.: Lat35.45S; Lon ​72.83W;​ Dep25.5 km +
-Origin time  : 2010/02/27 06:​34:​15.60 +
-Time delay   : 60.0  sec +
-Half duration: 60.0  sec +
-(…) +
-Best Double Couple: M0=1.84E+29 dyn.cm +
-NP1: Strike18 ; Dip=16 ; Slip=111  +
-NP2: Strike=176 ; Dip=75 ; Slip= 84  +
-             ​----###​-- ​             +
-         ​-------########​-- ​         +
-      ---------############​-- ​      +
-    -----------##############​-- ​    +
-   ​-----------################​-- ​   +
-  ------------#################​-- ​  +
- ​-------------#################​---  +
---------------######## ​  #######​--- +
-----   ​-------########​ T #######​--- +
----- P -------######## ​  #######​--- +
-----   ​-------##################​--- +
----------------################​---- +
- ​--------------################​---  +
-  -------------###############​--- ​  +
-   ​-------------#############​--- ​   +
-    ------------###########​---- ​    +
-      -----------########​---- ​      +
-         ​----------##​----- ​         +
-             #####​----+
  
 +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>​
  
-Wphase solution after grid search+If Basemap 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>​
  
-Each spatial and temporal location explored during the grid-search is detailed in grid_search_ts_out (time-shift grid search) grid_search_xyz_out (lat, lon, dep grid-search) and grid_search_xy_out (finer lat, lon grid-search which is performed at the optimum depth found during the lat, lon, dep global exploration). 
  
-The grid-search exploration result can be displayed using the python script make_grids.py :+-------
  
->​${WPHASE_HOME}/​bin/​make_grids.py –b 
  
-The optional parameter –b activates the use of the basemap toolkit in order to draw the topography and coastlines. The resulting figures shown here in Fig. 1 and Fig. 2 are written in the pdf files grid_search_ts.pdf and grid_search_xy.pdf.+===== Data formats =====
  
-There are two script which allows ​to compare predicted and observed waveforms.  +You can also have a look to [[wphase:tutorial|the W-phase tutorial page]].
-The first one shows the comparison between concatenated W phase traces using fort.15 files: +
- +
-> ${WPHASE_HOME}/​bin/​make_cwp.py +
- +
-The resulting figures are given in files CWP_R.pdf ​ CWP_W.pdf where the predicted waveforms in CWP_R.pdf are computed from the reference solution (e.g. CMTSOLUTION) ​ and the predicted waveforms in CWP_W.pdf are calculated from the W phase solution (either in WCMTSOLUTION,​ ts_WCMTSOLUTION or xy_WCMTSOLUTION)The LHZ traces contained in CWP_W.pdf are drawn in Fig. 3.  +
-The second script plot the waveform fit: +
- +
-> ${WPHASE_HOME}/​bin/​traces_global.py +
- +
-The resulting figure is given in file wp_pages.pdf which is shown here in Fig. 4. If the basemap python module is available, a map is drawn showing the station location with respect to the centroid location. Otherwise, a simple polar representation is used to display the station location. +
- +
-------- +
- +
-==== Note on file formats ==== +
  
 === CMTSOLUTION FILE === === CMTSOLUTION FILE ===
  
 example: example:
 +<​code>​
  PDE 2003  9 25 19 50  6.40  41.8100 ​ 143.9100 ​ 27.0 6.9 8.1 HOKKAIDO, JP  PDE 2003  9 25 19 50  6.40  41.8100 ​ 143.9100 ​ 27.0 6.9 8.1 HOKKAIDO, JP
 event name:     ​092503C ​       ​ event name:     ​092503C ​       ​
Line 288: Line 229:
 Mrp:       ​2.590000e+28 Mrp:       ​2.590000e+28
 Mtp:      -6.620000e+27 Mtp:      -6.620000e+27
 +</​code>​
  
 The first line of this file is the PDE.  The first line of this file is the PDE. 
-The first 4 characters of this line (including the first white space in the example above) generaly indicates the agency which provide the origin time, the hypocenter location and preliminary magnitude estimates.  +  * The first 4 characters of this line (including the first white space in the example above) generaly indicates the agency which provide the origin time, the hypocenter location and preliminary magnitude estimates.  
-Characters from 6 to 27 correspond to the PDE origin time +  ​* ​Characters from 6 to 27 correspond to the PDE origin time 
-Characters from 29 to 52 are the PDE latitude, longitude and depth. This hypocenter location will be used to define the time window in order to select the part of the waveform which will be inverted +  ​* ​Characters from 29 to 52 are the PDE latitude, longitude and depth. This hypocenter location will be used to define the time window in order to select the part of the waveform which will be inverted 
-Characters from 54 to 60 are PDE magnitude estimates The rest of the line provides some details on epicenter location (region, country, ...)+  ​* ​Characters from 54 to 60 are PDE magnitude estimates The rest of the line provides some details on epicenter location (region, country, ...)
  
 The second line, 'event name' correspond to the event id The second line, 'event name' correspond to the event id
 +
 Lines 3 and 4 are the parameters of the STF. 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). 
  
 +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).
  
  
Line 308: Line 250:
  
 This file is composed of several fields (some of them are optional): This file is composed of several fields (some of them are optional):
-EVNAME: Arbitrary event name +   ​* ​EVNAME: Arbitrary event name 
-SEED:  path to SEED file. This field must be used several times if there are multiple SEED files. +   * SEED:  path to SEED file. This field must be used several times if there are multiple SEED files. 
-DMAX:  Maximum distance to be used.  +   * DMAX:  Maximum distance to be used.  
-DMIN:  Minimum distance to be used (optional). ​ If omitted, it is defaulted to 0. +   * DMIN:  Minimum distance to be used (optional). ​ If omitted, it is defaulted to 0. 
-CMTFILE: Path to the CMTSOLUTION file +   * CMTFILE: Path to the CMTSOLUTION file 
- +   * filt_order: order of the butterworth filter. 
-filt_order: order of the butterworth filter. +   * filt_cf1: ​ low-frequency cut –off. 
-filt_cf1: ​ low-frequency cut –off. +   * filt_cf2: ​ high-frequency cut –off. 
-filt_cf2: ​ high-frequency cut –off. +   * filt_pass: unused parameter (let it to 1). 
-filt_pass: unused parameter (let it to 1). +   * IDEC_2: the second number is the number of seconds just before P wave over which the baseline is defined (let the other parameters to 2 and 0.1) . 
-IDEC_2: the second number is the number of seconds just before P wave over which the baseline is    +   * IDEC_3: Instrument response fit parameters for the deconvolution. The two first parameters define the frequency band in which the fit is measured. The third parameter is the number of samples used in this frequency range and the last parameter is a maximum misfit above which the channel is rejected. 
-      ​defined (let the other parameters to 2 and 0.1) . +   * WP_WIN: time window definition.  ​ 
-IDEC_3: Instrument response fit parameters for the deconvolution. The two first parameters define the frequency band in which the fit is measured. The third parameter is the number of samples used in this frequency range and the last parameter is a maximum misfit above which the channel is rejected. +     ​* ​If one value $B(eg. WP_WIN: 15) is used, then the window is defined by $[P_{tt},P_{tt}+B\times\Delta]where $P_{tt}$ ​is the P-wave travel-time.  ​ 
-WP_WIN: time window definition. ​ If one number ​B (eg. WP_WIN: 15) is used, then the window is defined by [Ptt,Ptt+B*∆] where Ptt is the P-wave travel-time. ​ If two numbers ​A B are specified (e.g., WP_WIN: 0. 15.), then the window is [Ptt+A*∆,Ptt+B*∆].  If three numbers ​A B C are used (e.g., WP_WIN: 0. 15. 12.) the time window is defined as:   +     ​* ​If two values $A$ $Bare specified (e.g., WP_WIN: 0. 15.), then the window is $[P_{tt}+A\times\Delta,P_{tt}+B\times\Delta]$.  ​ 
- [Ptt+A*∆,Ptt+B*∆] if >C  +     ​* ​If three values $A$, $B$ and $Care used (e.g., WP_WIN: 0. 15. 12.) the time window is defined as: 
- [Ptt+A*C,Ptt+B*C] if <C . +        * $[P_{tt}+A\times\Delta,P_{tt}+B\times\Delta]if $\Delta>C$ 
-If Four numbers ​A B C D are given (e.g., WP_WIN: 0. 15. 12. 50.) then the window is: +        * $[P_{tt}+A\times ​C,P_{tt}+B\times ​C]if $\Delta<C 
- [Ptt+A*∆,Ptt+B*∆] if C<<D  +     * If Four values, $A$, $B$, $C$ and $Dare given (e.g., WP_WIN: 0. 15. 12. 50.) then the window is: 
- [Ptt+A*C,Ptt+B*C] if <C , +        * $[P_{tt}+A\times \Delta,P_{tt}+B\times\Delta]if $C<\Delta<D$ 
- [Ptt+A*D,Ptt+B*D] if >D . +        * $[P_{tt}+A\times ​C,P_{tt}+B\times ​C]if $\Delta<C$ 
-GFDIR: Name of the Green'​s function directory (optional)+        * $[P_{tt}+A\times ​D,P_{tt}+B\times ​D]if $\Delta>D$ 
 +   * GFDIR: Name of the Green'​s function directory (optional)
  
 example: example:
 +<​code>​
 EVNAME: ​  ​Tokachi-Oki_2003 EVNAME: ​  ​Tokachi-Oki_2003
 SEED:    ../​../​WP6/​SEEDS/​2003_tokachi_oki_LH.SEED ​ SEED:    ../​../​WP6/​SEEDS/​2003_tokachi_oki_LH.SEED ​
Line 348: Line 292:
 WP_WIN: 15.  ​ WP_WIN: 15.  ​
 GFDIR: ./GF  GFDIR: ./GF 
- +</​code>​
  
  
Line 358: Line 301:
  
 Example: Example:
 +<​code>​
    ​60.0000 ​  ​0.20140377    ​60.0000 ​  ​0.20140377
    ​81.7621 ​  ​0.32070120    ​81.7621 ​  ​0.32070120
Line 383: Line 326:
 047 002    62.0000 ​   81.7621 ​  ​-35.8500 ​  ​-72.7100 ​   44.8000 ​  ​0.20324196 ​  ​0.31038266 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 048 002    64.0000 ​   81.7621 ​  ​-35.8500 ​  ​-72.7100 ​   44.8000 ​  ​0.20771868 ​  ​0.31790207
 +</​code>​
  
 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 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.
Line 397: Line 341:
  
 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 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
-> make_grids.py –t +<code> 
 +${WPHASE_HOME}/​bin/​make_grids.py –t  
 +</​code>​
 to create the file grid_search_ts.pdf containing Fig. 1 which corresponds to the above example. 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. 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.
 +<​code>​
   -35.4500 ​  ​-72.8325 ​   25.5000 ​  ​0.19118874   -35.4500 ​  ​-72.8325 ​   25.5000 ​  ​0.19118874
   -35.2500 ​  ​-72.7100 ​   25.5000 ​  ​0.19314386   -35.2500 ​  ​-72.7100 ​   25.5000 ​  ​0.19314386
Line 422: Line 365:
 105 002    60.0000 ​   60.0000 ​  ​-35.3500 ​  ​-72.7100 ​   25.5000 ​  ​0.19296721 ​  ​0.29330241 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 106 002    60.0000 ​   60.0000 ​  ​-35.5500 ​  ​-72.7100 ​   25.5000 ​  ​0.19360157 ​  ​0.29435001
 +</​code>​
  
 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 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. ​
Line 437: Line 381:
  
 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 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
-> make_grids.py –p -b +<code> 
 +${WPHASE_HOME}/​bin/​make_grids.py –p -b  
 +</​code>​
 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 basemap module in order to plot coastlines and topography. 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 basemap module in order to plot coastlines and topography.
- 
wphase/documentation.txt · Last modified: 2022/01/10 07:39 by wphase