writing catalog to NLL obs doesn't handle 12 (etc) channels properly

[filefolder]

Avoid duplicates

• I searched existing issues

Bug Summary

Noticing that when one exports a catalog object to NLL's obs format, it converts any channel outside of ZNEH to a ?. Not really an issue as NLL still handles it fine.

e.g. obspy/io/nlloc line 400 or so

    for pick in catalog[0].picks:
        wid = pick.waveform_id
        station = wid.station_code or "?"
        component = wid.channel_code and wid.channel_code[-1].upper() or "?"
        if component not in "ZNEH":
            component = "?"
        onset = ONSETS_REVERSE.get(pick.onset, "?")
        phase_type = pick.phase_hint or "?"
        polarity = POLARITIES_REVERSE.get(pick.polarity, "?")

However when reading NLL's output .hpy file back in, this (presumably horizontal) channel is now converted to N which then becomes a problem if the station's actual channels are 12.

It seems like a trivial fix to change if component not in "ZNEH" to at least ZNE123H (or even ZNE123HABCUVWTR?) unless I am missing something. I know NLL will handle 123 but not completely sure about the more esoteric channels.

Code to Reproduce

No response

Error Traceback

No response

ObsPy Version?

1.3.1.post0+390.g3ad5eee621.obspy.master

Operating System?

No response

Python Version?

No response

Installation Method?

developer installation / from source / git checkout

[filefolder]

2nd issue that is a bit more serious but also an easy fix in the same file...

obspy is reading in the phase information in hyp files incorrectly. For example, the top two lines correspond to the given arrivals.

PHASE ID    Ins Cmp   On Pha  FM   Date    HrMn   Sec     Err  ErrMag    Coda      Amp       Per       PriorWt  >   TTpred    Res       Weight    StaLoc(X  Y           Z)      SDist    SAzim  RAz  RDip RQual   Tcorr     TTerr
NWAO         ?    Z    ? P      ? 20221031 0502   41.3602 GAU -1.00e+00 -1.00e+00 -1.00e+00 -1.00e+00    1.0000 >   12.7040   -0.3041    1.2104  -24.3880 -103.3008   -0.3800   73.3529 107.42  40.0 200.0  0     0.0000    0.6352
NWAO         ?    ?    ? S      ? 20221031 0502   51.4602 GAU -1.00e+00 -1.00e+00 -1.00e+00 -1.00e+00    1.0000 >   22.7013   -0.2013    1.1613  -24.3880 -103.3008   -0.3800   73.3529 107.42  40.0 200.0  0     0.0000    1.1351

[Arrival
	   resource_id: ResourceIdentifier(id="smi:local/a9cfe6f2-8a0a-4030-aa61-3ec2601dd4e6")
	       pick_id: ResourceIdentifier(id="smi:local/cdffadb1-9954-4eaa-bbf9-57f1c8e23011")
	         phase: 'P'
	       azimuth: 107.42
	      distance: -0.003417422102491176
	 takeoff_angle: 40.0
	 time_residual: 12.704
	   time_weight: -0.3041,
 Arrival
	   resource_id: ResourceIdentifier(id="smi:local/c423d45d-94cb-4797-80b8-ae5cb4789a83")
	       pick_id: ResourceIdentifier(id="smi:local/f93e7b3e-0435-4c36-8f98-1fd2004093c5")
	         phase: 'S'
	       azimuth: 107.42
	      distance: -0.003417422102491176
	 takeoff_angle: 40.0
	 time_residual: 22.7013
	   time_weight: -0.2013,

last.hyp.txt

The distances, residuals, and weights have been jumbled. Possibly obspy is following an older format version of NLL HYP file (e.g. v6? I'm using v7)

In any event the following three changes in core.py produces the expected results

        arrival.distance = kilometer2degrees(float(line[22]))
        arrival.time_residual = float(line[17])
        arrival.time_weight = float(line[18])

[Arrival
	   resource_id: ResourceIdentifier(id="smi:local/58af4411-5af4-4a7c-952c-4ebd6b279f0e")
	       pick_id: ResourceIdentifier(id="smi:local/366a99cb-2f4e-4598-83ea-0879e5e14987")
	         phase: 'P'
	       azimuth: 107.42
	      distance: 0.6596784782679606
	 takeoff_angle: 40.0
	 time_residual: -0.3041
	   time_weight: 1.2104,
 Arrival
	   resource_id: ResourceIdentifier(id="smi:local/887cd0db-2ab7-404c-a250-744d86b0e75b")
	       pick_id: ResourceIdentifier(id="smi:local/da2aa00c-a3dc-4059-a706-d9e1c005c944")
	         phase: 'S'
	       azimuth: 107.42
	      distance: 0.6596784782679606
	 takeoff_angle: 40.0
	 time_residual: -0.2013
	   time_weight: 1.1613,

continuing my deep dive into this (previously I was converting these all by hand!) so possibly more to come

[megies]

Regarding the first issue

I totally agree, this probably should be changed to just use the last character of channel code unchanged. Looking at NLLoc docs, I probably was confused by component identification for the trace for which the time pick corresponds (i.e. Z, N, E, H), with "i.e." implying that it should be exactly one the mentioned characters, but probably what is meant is "e.g.". I would also toss out the capitalization, just use last character as is.

Actually, I'm looking at some of our test files from @claudiodsf (io/nlloc/tests/data/nlloc.hyp) and I see that it has all three letters of the channel code in this "component" field, so probably we can just do the same.

Probably just make it:

    component = wid.channel_code or "?"

--

Regarding the second issue

well yeah, thats just very unfortunate to change the NLLOC hyp file format which is free form whitespace separated values in a breaking fashion. There apparently was one value added in the middle of the phase lines (PriorWt).

Original

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

New (NLL_FORMAT_VER_2 ?)

PHASE ID Ins Cmp On Pha FM Date HrMn Sec Err ErrMag Coda Amp Per PriorWt > TTpred Res Weight StaLoc(X Y Z) SDist SAzim RAz RDip RQual Tcorr TTerr

Is there a way to tell from the file that it is "version 2"? Can you post one of your hyp files? If yes, then we can determine the indices of items on the phase lines from it, if not we might have to check at what position in the line this magic ">" separator is positioned to judge from that what indices to look at.

[filefolder]

yep i was afraid of that...

I did post the ver2 hyp file above, it's easy to miss https://github.com/obspy/obspy/files/10028900/last.hyp.txt

if V2 has more elements i suppose that is one tell

[megies]

I guess we could check the NLLoc version as stated in the SIGNATURE line, I think it should always be there. But then we would need to check at what exact nonlinloc version this changed.. to be honest, I'd probably check if "PriorWt" in line PHASE.. or well, actually, it feels safest to just look at position of ">" magic character in the PHASE line

[filefolder]

agree, probably easiest and fastest to if line[112] == '>' then for the v2 shifts