problem with rotate method and function

[pavlis]

I discovered what could either be called a bug or a feature in the rotate method of Seismogram. It carries over to the rotate function which has a wrapper defined for map operators.

If could be called a "feature" if I make the documentation clearer. It turns out the angle of rotation is positive clockwise. That is the standard convention for azimuth angles, but is reversed for spherical coordinates. i.e. in very math course the angle phi spherical coordinates is positive in the anticlockwise direction. Stock geographic coordinates we use in seismology always have the x (x1) axis positive in the local east direction and y (x2) positive in the north direction. This small code fragment will show you that the convention used for rotate is positive clockwise and measured in radians:

from mspasspy.ccore.seismic import Seismogram
from mspasspy.algorithms.basic import rotate_to_standard,rotate
import numpy as np
import math

d=Seismogram(4)
d.dt=1.0
d.set_live()
phid=30.0
for j in range(d.npts):
    d.data[0,j]=math.cos(np.radians(phid))
    d.data[1,j]=math.sin(np.radians(phid))
phi=-30.0
print(d.data)
dtrans = rotate(d,np.radians(phi))

print(dtrans.data)

If you run that you'll see the input data has vectors that are [sqrt(2),0.5,0] and after rotate they become [1,0,0]. The math side is clear. For seismology this would come up if the "radial" direction for local wave propagation is at 60 degrees azimuth ( convention with R + in the direction of propagation. In that case P waves would have horizontal motion mainly in the R-Z plane. Standard vertical rotation convention is RTZ meaning x,y,z->R,T,Z. That means after the rotation x (x1) should point in the R direction. Confusing signs, I know, but the little test script shows the way you would transform the coordinates to RTZ would be to compute the azimuth difference between E (azimuth of 90) and the direction of propagation (+R). For above the azimuth of the propagation direction is 60 degrees so the azimuth relative to E is -30 degrees because azimuth is measured positive clockwise.

The other perspective is this is a bug from mixing up the forward an inverse transform of the transformation matrix for this problem. The standard formula for rotation by angle phi with the normal math convention is:

|  a   -b   0 |
|  b   a    0|
| 0    0     1|

where $a=cos( \phi )$ and $b=sin(\phi)$. The C++ code for rotate has the signs reversed. Why it is more that a little confusing is that if you use the azimuth sign convention the signs need to be reversed since sin is antisymmetric in phi. Furthermore, the transform of a rotation matrix is it's inverse, so yet another perspective is there is a bug using the transform of the matrix by mistake.

Why this isn't cut and dried is I'm not at all sure whether it is better to leave the code alone and make the documentation clear on this point or if change the C++ code and the make it clearer in the documentation that rotate uses math sign convention. Either way we very very much need to address this issue. @wangyinz let me know which I should do.

[wangyinz]

So, I don't have a preference here, but I wonder how this is handled in other packages, probably ObsPy.

Well, after looking at their rotate method, it seems we don't really have a direct comparison here...

[pavlis]

Good point to see what obspy does. They in fact, suggest what we should probably do to address this issue. That is, we probably need a wrapper that simplifies this confusing (to me anyway) problem that is a central issue in three component data. If you look what they do their "rotate" should actually probably be called "tranform" because they take the input and compute the transformation matrix to transform the data into some standard coordinate systems that have been defined over the years. Notably, RTZ, LQT, and some left handed version of same. I think I will write a pair of wrappers called RTZ_transform and LQT_transform that use the rotate and/or tranform methods of Seismogram to simplify the process. I think I'll make them so the required attributes can come from either the arg list or Metadata using the already defned keys "seaz" and "ema" (well, I think they are defined. They are standard CSS3.0 keys).

As to other implementations, from my experience we are not alone in having confusions on this point. The azimuth versus phi convention is the fundamental problem that confuses a lot of stuff like this in our field.

In summary, if you guys concur I'll go forward with building those wrappers and relegate the the lower level rotate and transform to the more primitive status they actually are.

[pavlis]

If we go ahead with the suggested wrappers I made above for RTZ and LQT coordinates, we should give some thought to what the names of the functions should be. Here are some options:

1. LQT_transform and RTZ_transform
2. Cryptic/symbolic: LQT and RTZ
3. Verbose: transform_to_LQT and transform_to_RTZ.

0ther ideas?

I personally suggest 3 as I can type fast, but some prefer the more cryptic. 3 is probably more consistent with the existing higher level functionality with the name rotate_to_standard that returns the data to standard ENZ coordinates. It also is consistent with the verbose name free_surface_transformation which is the 3rd common ray-based coordinate system.

[wangyinz]

I think verbose is better too

[pavlis]

I think this issue can be close.