azimuth calculated wrongly in flinn (signal/polarization)

[ftilmann]

In function flinn there is the following line

        if azimuth > 180.0:
        azimuth -= 180.0

As azimuth can legitimately vary between 0-360 (-180 - 180) I believe this should read

        if azimuth > 180.0:
        azimuth -= 360.0

Here is the link to the line in the code https://github.com/obspy/obspy/blob/1c656830a33e72b7d335f01d28c38b97d8c4d58a/obspy/signal/polarization.py#L167

There is another problem with this routine: it measures incidence as the angle from vertical up, so that (for a P wave) the azimuth is the (backazimuth +- 180). For seismological analysis (e.g. array seismology) measured azimuth is generally expected to point back at the source, and incidence is the angle with respect to vertical down. This effectively means that the 180° needs to be added/subtracted from the azimuth returned (after correction above). In my view this would be achieved most transparently by adding a minus sign in either https://github.com/obspy/obspy/blob/1c656830a33e72b7d335f01d28c38b97d8c4d58a/obspy/signal/polarization.py#L153 or https://github.com/obspy/obspy/blob/1c656830a33e72b7d335f01d28c38b97d8c4d58a/obspy/signal/polarization.py#L155

e.g. write

    incidence = math.degrees(math.atan2(eve, -eigvec[2][0]))

instead of

    incidence = math.degrees(math.atan2(eve, eigvec[2][0]))

[megies]

CC @jwassermann

[jwassermann]

Well, it is a pretty old code but I think it is doing the right thing (what I intended). First of all I use in all codes (or try to) the azimuth and incidence as ray-based measures. Azimuth is in this case measured as angle via north of the vector pointing from source to the station and vertical up for the incidence. Back-azimuth is then defined as the angle via north of the vector pointing from the station back to the source. In case of unknown polarity (compression or dilatation) and unknown source position these different definition do not matter as there is a clear 180 deg ambiguity. And exactly for this situation the code was written originally (source tracking at a volcano). You might be right that in obspy a more general definition serves better to a larger audience, but as I mentioned the code is quite old and was adapted to python more than a decade ago. Probably this should be changed or better described in a new release.

[ftilmann]

@jwassermann I can follow your explanations with regard to the definition (the second point of the issue starting with "There is another problem ...")--although I think it should be changed---, but it does not address the azimuth calculation, and I doubt the code as currently written does what you want. For a P wave, there is no 180° ambiguity even in the case of unknown first motion polarity - the vertical down part of the particle motion will have an azimuth pointing towards the source, and the vertical up part will have an azimuth pointing away from the source. So if you limit angle of incidence to 0-90°, as in https://github.com/obspy/obspy/blob/1c656830a33e72b7d335f01d28c38b97d8c4d58a/obspy/signal/polarization.py#L160 and following lines, azimuth must be able to vary between 0° and 360° (or -180 to 180°) equivalently to azimuth and plunge for line feature in geology). [An alternative representation would have azimuth between 0-180° and incidence between 0 and 180° but this is less useful for seismological applications, and also not used in flinn]. So subtracting 180° from azimuth without flipping incidence is just wrong I think.

[ftilmann]

If you would like to keep the definition with incidence vertical up, I propose the following explanatory text:

"Note that incidence is measured with respect to vertical up and can take values between 0-90°. This implies that when applied to an ideal P wave with linear particle motion, the azimuth returned will be the backazimuth±180°"

[Note that it will generally NOT be the source azimuth due to Earth curvature. Of course in volcano seismology the curvature can be neglected. ] The sentence will only be true after correcting the azimuth calculation.

[zenytaa]

i have a problem when calculating rockfalls azimuth ( that confirmed head into southwest based on visual observation)*, the result is in the range 165-167 degrees using the "stock" calculation from obspy, which is wrong because the station is in the northeast from the event. but when I change calculation just like @ftilmann does, the result is still the same (of course). is it necessary to add a minus in both azimuth and incidences calculation? or just one of them?

[megies]

So would it be acceptable to add the following clarifying sentence @ftilmann proposed and then close this, @jwassermann?

"Note that incidence is measured with respect to vertical up and can take values between 0-90°. This implies that when applied to an ideal P wave with linear particle motion, the azimuth returned will be the backazimuth±180°"

[ftilmann]

@megies There are two sub-issues I raised:

1. The azimuth calculation as currently implemented is wrong - a code change is needed here, and this cannot be addressed with a documentation change alone.
2. I proposed a different definition of incidence angle (with 0° representing vertical down instead of vertical up), which I feel is more geared to general seismological applications, as the azimuth will then represent the backazimuth, which I believe is also how array processing is usually done. This is ultimately a matter of opinion, though, so if obspy developers do not agree, then I proposed the quoted documentation change.

[megies]

To be honest, I don't know how to resolve this at the moment. To my impression there are two different opinions of the use case of the code and the definitions going along with it. I don't know this part of the code, so I'm not keen to touch anything and might have to delay this to some later release.