FocMec setup

[AbeJones]

Hi Tobais,

I have downloaded FocMec from IRIS and placed in within my Obspyck dir. I was just wondering what needs to be done in terms of bash scripts so that focmec can be called simply by clicking on do focmec in the Obspyck interface ?

[megies]

Hi Abraham,

you need three pieces in a focmec directory in your plugin_dir:

• the focmec binary (or a symlink to it)
• a (e.g.) bash script that handles calling focmec (mainly because focmec is called with a "maximum polarity errors" option and the solutions with less errors simply get drowned in a multitude of worse solutions; so we need to do several focmec runs starting at 0 maximum polarity errors and then going higher until at least one solution is found; e.g. here)
• an focmec control/input file (with commands/settings to be piped into focmec as stdin, e.g. here)
• the file with picks/polarities (will be written by obspyck from the picks currently in the GUI)

I've added my setup in e9c1fe9. Let me know if anything is unclear.

(Edit: added some more clarifications)

[klaaur]

Does obspyck create the input files from the picks in the GUI or does the user need to supply focmec input files for each event?

[megies]

Input files with picks are of course generated from what is currently set on the waveforms in the GUI.
Some notes:

• after setting picks that should be used for focmec a location run has to be done (focmec needs azimuth and takeoff angles)
• a pick that is changed (polarity changed, time changed) becomes an entirely new pick, so after any changes run location routine again
• to set S-wave polarities on rotated components (polarities left/right, forward/backward on T/R components)..
  • set some picks and make an initial location
  • remove any unwanted picks again (e.g. S picks on unrotated components)
  • click "rotate ZNE"
  • set any wanted S picks with polarities on rotated components
  • run location again
  • run focmec

I've added some clarifications in my above post

[klaaur]

Thanks for the update above, that helped clarify things greatly.

[klaaur]

When I run 'show FocMec' I get the following error:

//anaconda/lib/python2.7/site-packages/obspy/core/util/attribdict.py:97: UserWarning: Setting attribute "_beachball" which is not a default attribute ("creation_info", "method_id", "moment_tensor", "resource_id", "evaluation_status", "station_distribution_ratio", "comments", "triggering_origin_id", "waveform_id", "nodal_planes", "misfit", "station_polarity_count", "azimuthal_gap", "principal_axes", "evaluation_mode").
warnings.warn(msg)
//anaconda/lib/python2.7/site-packages/obspy/core/util/attribdict.py:97: UserWarning: Setting attribute "_beachball2" which is not a default attribute ("creation_info", "method_id", "moment_tensor", "resource_id", "evaluation_status", "station_distribution_ratio", "comments", "triggering_origin_id", "waveform_id", "nodal_planes", "misfit", "station_polarity_count", "azimuthal_gap", "principal_axes", "evaluation_mode").
warnings.warn(msg)
Traceback (most recent call last):
File "obspyck.py", line 508, in on_qToolButton_showFocMec_toggled
self.drawFocMec()
File "obspyck.py", line 1992, in drawFocMec
self.canv.draw()
File "//anaconda/lib/python2.7/site-packages/matplotlib/backends/backend_qt5agg.py", line 159, in draw
FigureCanvasAgg.draw(self)
File "//anaconda/lib/python2.7/site-packages/matplotlib/backends/backend_agg.py", line 474, in draw
self.figure.draw(self.renderer)
File "//anaconda/lib/python2.7/site-packages/matplotlib/artist.py", line 62, in draw_wrapper
draw(artist, renderer, _args, *_kwargs)
File "//anaconda/lib/python2.7/site-packages/matplotlib/figure.py", line 1159, in draw
func(_args)
File "//anaconda/lib/python2.7/site-packages/matplotlib/artist.py", line 62, in draw_wrapper
draw(artist, renderer, *args, *_kwargs)
File "//anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.py", line 2319, in draw
a.draw(renderer)
File "//anaconda/lib/python2.7/site-packages/matplotlib/artist.py", line 62, in draw_wrapper
draw(artist, renderer, _args, *_kwargs)
File "//anaconda/lib/python2.7/site-packages/matplotlib/text.py", line 757, in draw
raise ValueError("posx and posy should be finite values")
ValueError: posx and posy should be finite values

I realize this is a matplotlib/text.py error during drawing but do you have any idea why this would be occurring and how to fix it? Lastly, I noticed that it always picks the first solution for the focal mechanism to show, but you can choose to push the 'next FocMec' button to see the next solution. When I push the button I get the following error:

selecting Focal Mechanism No. 2 of 20:
Strike: 180.00 Dip: 85.00 Rake: -90.00 Misfit: 0.00
Traceback (most recent call last):
File "obspyck.py", line 528, in on_qToolButton_nextFocMec_clicked
self.drawFocMec()
File "obspyck.py", line 1923, in drawFocMec
ax.add_collection(fm_._beachball2)
File "//anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.py", line 1668, in add_collection
self._set_artist_props(collection)
File "//anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.py", line 860, in _set_artist_props
a.axes = self
File "//anaconda/lib/python2.7/site-packages/matplotlib/artist.py", line 240, in axes
raise ValueError("Can not reset the axes. You are "
ValueError: Can not reset the axes. You are probably trying to re-use an artist in more than one Axes which is not supported

Am I misusing the button or is there another way to show another solution? If I'm not misusing it do you know how to fix this error? Thanks.

[megies]

Hmm.. I honestly haven't used focmec in a while.. Do you have an example with openly accessible data that I could run on my machine?

[klaaur]

Sure. I'm pulling from the FDSN server and using stations for the Tohuku event (as is your default time in the .obspyckrc file. Thus use your example time in the example.cfg file ):

Change your EXAMPLE1 to the following stations:

EXAMPLE1=IU.DAV.00.BH_,IU.MIDW.00.BH_,IU.TIXI.00.BH_,IU.BILL.00.BH_,IU.COLA.10.BH_,IU.RSSD.10.BH_,IU.SFJD.00.BH_,IU.JOHN.00.BH_,IU.FUNA.00.BH_,IU.CTAO.00.BH_,II.TLY.00.BH_,II.AAK.00.BH_,II.OBN.00.BH_,II.MSVF.00.BH_,II.WRAB.10.BH_,II.ARU.00.BH_,G.INU.00.BH_,IU.MAJO.00.BH_,II.ERM.00.BH_,JP.JHJ2..BH_,G.KIP.00.BH_,IU.YSS.00.BH_,IC.MDJ.00.BH_,IU.INCN.00.BH_,IU.TATO.00.BH_,JP.JOW..BH_,IU.GUMO.00.BH_,HK.HKPS.00.BH_,IC.BJT.00.BH_,IU.PET.00.BH_,IU.ULN.00.BH_,IU.KWAJ.00.BH_,MI.ANNE..BH_,IU.YAK.00.BH_,IU.MA2.11.BH_,IU.POHA.10.BH_,IU.PMG.00.BH_,AT.SMY..BH_,IC.HIA.00.BH_,IC.ENH.00.BH_,AK.ATKA..BH_,AK.SII..BH_,CI.PHL..BH_,AU.MANU..BH_,JP.CBIJ..BH_,JP.YOJ..BH_,KC.TARG..BH*

You could shorten the station list and just use the first 5 stations if you want. Either way with fewer picks and/or more picks I get the same error. I'm trying to troubleshoot the output more and see if I can narrow down what might be causing it. I'm using a global nlloc configuration file (ak135 velocity model) from the NNLOC examples but can pass my configuration file to you if you'd like too. Let me know.

[klaaur]

I believe I've narrowed it down to the fact that the incidence angles aren't being created, if I'm correct in the assumption that the phase file (input focmec file) should have 4 columns with the following information:

col1 col2 col3 col4
sta azimuth incidence angle polarity

I'm missing the incidence angle, and so the array of incidence angles is full of nans, hence probably why the scatter plot and self.canv.draw() lines have issues with the incis[mask]. Do the incidence angles come from NLLOC or are they internal to the obspyck? E.g., Here's an output of the file:

Phases for focmec: 5

AAK 297.48 nanU
ARU 318.08 nanU
ATKA 411.34 nanD
BILL 16.65 nanU
BJT 281.82 nanU

The question though is how can I fix this?

[megies]

incidence angles aren't being created

Hmm.. you only get incidence/takeoff angles and backazimuth if you also calculate the angle cubes for NonLinLoc, that is you need ....angle.buf files in addition to .....time.buf files. Maybe that's the problem, that you don't have those?

But there should be a better error message and an appropriate Exception being raised by obspyck..

[klaaur]

You're correct. I emailed Anthony Lomax to ask how to add the takeoff angles in global mode of nlloc, as the LOCANGLES_YES parameter doesn't calculate them without the angle files. If anyone else has this problem they can do the following: To get angles you need to add the following lines in /nlloc_global_sample/taup/TauP_Table_NLL.sh
echo First arriving P angles
java edu.sc.seis.TauP.TauP_Table_NLL -nll_angles ${NLL_GRID} -mod ${MODEL} -ph P,Pg,p,Pn,Pdiff,PKP,PKiKP,PKIKP -o ${OUTPATH}/${MODEL}
then re-run TauP_Table_NLL.sh (or the whole example script).

Then you need to set LOCANGLES_YES parameter in the input configuration file. Everything else should work within obspyck thereafter.

[klaaur]

In terms of my other issue raised above, when selecting the 'next FocMec' button after running the 'show FocMec' button I get the following error:

selecting Focal Mechanism No. 2 of 20:
Strike: 180.00 Dip: 85.00 Rake: -90.00 Misfit: 0.00
Traceback (most recent call last):
File "obspyck.py", line 528, in on_qToolButton_nextFocMec_clicked
self.drawFocMec()
File "obspyck.py", line 1923, in drawFocMec
ax.add_collection(fm_._beachball2)
File "//anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.py", line 1668, in add_collection
self._set_artist_props(collection)
File "//anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.py", line 860, in _set_artist_props
a.axes = self
File "//anaconda/lib/python2.7/site-packages/matplotlib/artist.py", line 240, in axes
raise ValueError("Can not reset the axes. You are "
ValueError: Can not reset the axes. You are probably trying to re-use an artist in more than one Axes which is not supported

Am I misusing the button or is there another way to show another solution? E.g., if I don't want the first solution do I select 'next Focmec' prior to pushing the 'show FocMec' button? If I'm not misusing it do you know how to fix this error? Thanks.

[megies]

How big is the nonlinloc cube for one station? Is it some manageable file size? (Kinda lazy to go through those nlloc steps..)

[klaaur]

NLLOC "./nlloc.20110311.055059.grid0" "LOCATED" "Location completed."
SIGNATURE "Anthony Lomax (www.alomax.net) NLLoc:v6.00.0 24Oct2016 14h35m22"
COMMENT ""
GRID 361 181 2 -180 -90 0 1 1 1 PROB_DENSITY
SEARCH OCTREE nInitial 27648 nEvaluated 50016 smallestNodeSide 104.166667/104.166667/0.083333 oct_tree_integral 1.223155e+06 scatter_volume 1.222790e+06
HYPOCENTER x 142.969 y 37.9688 z 0.0416667 OT 17.7816 ix -1 iy -1 iz -1
GEOGRAPHIC OT 2011 03 11 05 46 17.7816 Lat 37.9688 Long 142.969 Depth 0.0416667
QUALITY Pmax 0.362529 MFmin 1.01998 MFmax 4.96789 RMS 1.3632 Nphs 5 Gap 268.084 Dist 20.8674 Mamp -9.9 0 Mdur -9.9 0
VPVSRATIO VpVsRatio -1 Npair 0 Diff -2e+30
STATISTICS ExpectX 72.2923 Y 54.1571 Z 1.00591 CovXX 4.0843e+07 XY -5.08622e+06 XZ 543.128 YY 7.12582e+06 YZ 14.1591 ZZ 0.333243 EllAz1 49.7212 Dip1 -89.9989 Len1 1.07108 Az2 188.394 Dip2 -0.000838462 Len2 4743.91 Len3 1.211715e+04
STAT_GEOG ExpectLat 54.1571 Long 72.2923 Depth 1.00591
TRANSFORM GLOBAL
QML_OriginQuality assocPhCt 5 usedPhCt 5 assocStaCt -1 usedStaCt 5 depthPhCt -1 stdErr 1.3632 azGap 268.084 secAzGap 282.304 gtLevel - minDist 2318.6 maxDist 6258.24 medDist 3648.46
QML_OriginUncertainty horUnc -1 minHorUnc 3829.25 maxHorUnc 9780.86 azMaxHorUnc 98.3943
FOCALMECH Hyp 37.9688 142.969 0.0416667 Mech 0 0 0 mf 0 nObs 0
PHASE ID Ins Cmp On Pha FM Date HrMn Sec Err ErrMag Coda Amp Per > TTpred Res Weight StaLoc(X Y Z) SDist SAzim RAz RDip RQual Tcorr
BJT ? ? ? P ? 20110311 0550 59.67 GAU 9.8 -1 -1 -1 > 283.5067 -1.6152 0.9205 116.1679 40.0183 -0.1970 20.8674 284.04 0.0 0.0 0 0.0000
BILL ? ? ? P ? 20110311 0552 52.4 GAU 2.4 -1 -1 -1 > 394.9209 -0.3041 1.2594 166.4531 68.0653 -0.3200 32.8012 15.96 0.0 0.0 0 0.0000
ATKA ? ? ? P ? 20110311 0552 56.31 GAU 8.9 -1 -1 -1 > 395.2278 3.2990 0.9484 -174.1975 52.2016 -0.0550 32.8362 410.26 0.0 0.0 0 0.0000
AAK ? ? ? P ? 20110311 0555 21.92 GAU 17 -1 -1 -1 > 543.7225 0.4154 0.6518 74.4942 42.6375 -1.6330 51.0235 298.26 0.0 0.0 0 0.0000
ARU ? ? ? P ? 20110311 0556 1.278 GAU 4.7 -1 -1 -1 > 582.5642 0.9326 1.2199 58.5625 56.4302 -0.2500 56.3242 318.58 0.0 0.0 0 0.0000
END_PHASE
END_NLLOC

My oct tree search is set to go through LOCSEARCH OCT 96 48 6 0.05 50000 10000 4 0

Was this the information you're looking for or simply how big the buffer/header files are?

[megies]

or simply how big the buffer/header files are?

Yep, I was interested in file size of the travel time cubes (*.buf).