Fix rotations again

[kujaku11]

This PR is to address issue #55 where rotations are currently set to rotate counterclockwise. The correct rotation should be clockwise. This seems to be a reoccurring issue and needs to be addressed properly so it doesn't' happen again, or at least there are some checks.

• Change rotation matrix to be clockwise positive [cos sin][-sin cos]
• Update tests to test for correct rotation
  • Rotation of Z
  • Rotation of Tipper
• Add documentation

[codecov]

Codecov Report

Attention: Patch coverage is 87.80488% with 30 lines in your changes missing coverage. Please review.

Project coverage is 33.66%. Comparing base (8a2cb64) to head (ad24576).
Report is 47 commits behind head on updates.

Files with missing lines	Patch %	Lines
mtpy/core/mt_data.py	63.63%	8 Missing ⚠️
mtpy/core/mt.py	79.31%	6 Missing ⚠️
mtpy/utils/calculator.py	83.87%	5 Missing ⚠️
mtpy/processing/aurora/process_aurora.py	50.00%	4 Missing ⚠️
mtpy/imaging/plot_phase_tensor_maps.py	0.00%	2 Missing ⚠️
mtpy/core/transfer_function/base.py	88.88%	1 Missing ⚠️
...py/core/transfer_function/z_analysis/distortion.py	50.00%	1 Missing ⚠️
mtpy/imaging/plot_phase_tensor_pseudosection.py	0.00%	1 Missing ⚠️
mtpy/processing/kernel_dataset.py	91.66%	1 Missing ⚠️
tests/utils/test_calculator.py	96.87%	1 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           updates      #57      +/-   ##
===========================================
+ Coverage    33.24%   33.66%   +0.42%     
===========================================
  Files          146      146              
  Lines        23983    24144     +161     
===========================================
+ Hits          7972     8129     +157     
- Misses       16011    16015       +4     

Flag	Coverage Δ
tests	33.66% <87.80%> (+0.42%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[kujaku11]

@alkirkby I've update the rotation matrix and added tests to make sure the rotation is correct. Can you have a look when you get a chance. Rotation of Z and Tipper should be correct now, could you verify that they are correct now?

[alkirkby]

@kujaku11 Cool, thanks for fixing- that seems to work now.

One additional weird thing that might be worth addressing now is that the mt_object.rotation_angle just gets updated to the last angle that it was rotated to. And rotation angles are not read from file. So if I read in an MT object rotated to 115, I get mt_obj.rotation_angle = 0.

Then if I read in an unrotated MT file and rotate by 115, I get mt_obj.rotation_angle=115. Then if I rotate by -115 I get mt_obj.rotation_angle = -115, where it should really be 0.

[kujaku11]

@alkirkby thanks, I need to have another look at the math to make sure that the rotations are proper, but good to hear it seems to work. Just to confirm you are comparing rotation of a transfer function using mtpy and the geotools, and now they look the same?

Good catch, so what the rotation angle should do is update from the previous angle. Like you say if you start out with 10 then rotate by 20, the rotation angle should be 30 instead of 20. I can fix that too.

As for reading in the rotation angle, I'll have a look at that as well.

[kujaku11]

@alkirkby or @kkappler Could you confirm the math here.

For rotation of the impedance tensor we want to apply a similarity transformation so that the properties remain the same.

We set the rotation matrix R to be the clockwise rotation matrix:

$$ \mathbf{R} = \begin{bmatrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{bmatrix}$$

Then if we formalize the similarity transformation as:

$$ \mathbf{A}' = \mathbf{R}^{T} \mathbf{A} \mathbf{R}$$

If we have an impedance tensor with a strike of zero and apply a rotation of 30 degrees the new strike is 30 degrees.

If we formalize the similarity transformation as:

$$ \mathbf{A}' = \mathbf{R} \mathbf{A} \mathbf{R}^{T}$$

Then if we have a strike of 0 and apply a rotation of 30 degrees the new strike is -30.

Which seems incorrect. I would think that if you rotate by a positive number then you should add that number to the existing strike angle. But maybe I'm wrong. @alkirkby can you confirm with this version of the branch that the rotation is the same as that from what you get from geotools?

From what I've seen in other codes is that they use the 2nd similarity transformation. that is what the current branch has.

[alkirkby]

@alkirkby thanks, I need to have another look at the math to make sure that the rotations are proper, but good to hear it seems to work. Just to confirm you are comparing rotation of a transfer function using mtpy and the geotools, and now they look the same?

Good catch, so what the rotation angle should do is update from the previous angle. Like you say if you start out with 10 then rotate by 20, the rotation angle should be 30 instead of 20. I can fix that too.

As for reading in the rotation angle, I'll have a look at that as well.

Hi - yes that's right, I rotated an edi file in geotools, and compared it to the result in mtpy, and they were the same

[alkirkby]

@alkirkby or @kkappler Could you confirm the math here.

For rotation of the impedance tensor we want to apply a similarity transformation so that the properties remain the same.

We set the rotation matrix R to be the clockwise rotation matrix:

R = [ cos ⁡ θ sin ⁡ θ − sin ⁡ θ cos ⁡ θ ]

Then if we formalize the similarity transformation as:

A ′ = R T A R

If we have an impedance tensor with a strike of zero and apply a rotation of 30 degrees the new strike is 30 degrees.

If we formalize the similarity transformation as:

A ′ = R A R T

Then if we have a strike of 0 and apply a rotation of 30 degrees the new strike is -30.

Which seems incorrect. I would think that if you rotate by a positive number then you should add that number to the existing strike angle. But maybe I'm wrong. @alkirkby can you confirm with this version of the branch that the rotation is the same as that from what you get from geotools?

From what I've seen in other codes is that they use the 2nd similarity transformation. that is what the current branch has.

Hi- yep that's right. Based on my testing yesterday I got the same from geotools as MTPy.
According to Simpson and Bahr, the transform to rotate clockwise should be:
A ′ = R A R T

Yes I agree, if you rotate by a number then that number should be added to the strike.

[kujaku11]

@alkirkby Great, thanks for the confirmation. I'm still gonna do some testing and add some tests so that we can confidently say its correct.

In Geotools can you estimate geoelectric strike? Just curious what the strike angle is if you rotate say by 30. Is the new strike angle plus 30?

[kujaku11]

@alkirkby, @kkappler So with the current formulation, that appears in nearly all MT basics

$$Z_{rotated} = \mathbf{R} \mathbf{Z} \mathbf{R}^T$$

where:

$$ \mathbf{R} = \begin{bmatrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{bmatrix}$$

Now if we apply this with a positive angle we would expect that the new strike angle would be original + new. However, it seems to be the opposite.

For example if we have a general 2x2 matrix with a strike of 45:

$$ \mathbf{Z} = \begin{bmatrix} \sqrt{2} & 1\\ 1 & \sqrt{2} \end{bmatrix}$$

Then rotate by 30 degrees:

$$ \mathbf{Z} = \begin{bmatrix} 2.28 & .5 \\ .5 & .55 \end{bmatrix}$$

the strike angle is now 15.

Strike is calculated as

$$\phi_{strike} = \dfrac{1}{2} \tan{\dfrac{a[0, 1] + a[1, 0]}{a[0, 0] - a[1, 1]}}$$

And even when I do this for a known model response I get the same, rotation counterclockwise

Original data with a NE striking feature at depth

Rotated by +30: expect the strike to rotate clockwise not counterclockwise

So what I'm I missing here, is it how the strike angles are defined? But those are the same as what's in the literature. Is it the reference frame that we are using to measure the angle from? Any thoughts would be helpful.

[kkappler]

@kujaku11 @alkirkby

I committed a notebook with some rotations notes that we can use as a starting point for conventions and maybe a tutorial in future.

There are a few conventions here, that I label in the notebook (coordinate frame, etc.).

There are two functions in the notebook that seem important rotate_field_to_cardinal, and rotate_cardinal_to_field, which are just inverses of each other.

An excerpt, that may key in on the confusion:

Writing down the definition of the strike angle is would help to totally dial this in.

[alkirkby]

Thanks @kkappler, that clears things up nicely. We need to be careful when talking about strike and rotation angle. My earlier comment "Yes I agree, if you rotate by a number then that number should be added to the strike." is incorrect. What I should have said is ".....that number should be added to the rotation angle"

So, if an edi file is read in with a rotation angle of 0, then we rotate it by 30 degrees, then mtobj.rotation_angle goes from 0 to 30.
If the edi file is read in with a rotation angle of 30, then this information should be read into the mt object such that it also has a rotation angle of 30.

The way I see it, the strike of that mt object shouldn't change, no matter what the rotation angle. If our hypothetical mt object had a strike of 20, then we rotate it 30 degrees, the strike is now -10 (relative to our new rotated measurement frame) but it's still 20 degrees relative to north. So maybe we need to include rotation angle in our calculation of strike.

And yes I agree, whatever we do we should document it in the code!

[alkirkby]

Hi @kkappler, @kujaku11
OK, understood. But what if I just want to rotate by an angle? The way I normally use MTPy for rotations is to simply rotate the data. I determined my data has a strike angle of 101 degrees (east of north) and wanted to rotate the data by 101 degrees east. Can I just go mt_obj.rotate(101) and it will apply the correct rotation matrix? Which I believe is:

$$ \mathbf{R} = \begin{bmatrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{bmatrix}$$

??

[kujaku11]

@alkirkby , @kkappler That's a good question. I would try using this branch and use the following code. And maybe try the same with Geotools and see which one looks correct. This is my confusion as well. So it would be nice to get some verification. I hope that the NED is correct. [Be sure to pull this branch before you try it, I just updated it.]

from mtpy import MT

m1 = MT()
m1.read(r"edi_file")
m1.coordinate_reference_frame = "ned"

m_rot_ned = m1.rotate(110, inplace=False)

m1.coordinate_reference_frame = "enu"

m_rot_enu = m1.rotate(110, inplace=False)

m_rot_ned.plot_phase_tensor(fig_num=1)
m_rot_enu.plot_phase_tensor(fig_num=2)

The correct one should have a strike of 0 in the new reference frame.

Strike angle confusion

Ah, I think I may have confused myself. I was thinking that if you rotate a coordinate reference frame clockwise that you would add the rotated angle to the strike. But I don't think that's the case. As Alison had said the strike angle should stay the same its just the reference frame is rotating. So if you rotate by $\alpha$ and the original strike is $\theta_{strike}$ then the new strike angle will be $\theta_{new} = \theta_{strike} - \alpha$ which means the strike angle will be smaller in the new reference frame. Is that correct? This would mean an R as Alison defines would be correct (clockwise rotation). I changed the code to mirror that.

So if you want to rotate to strike you rotate by the strike angle? If you have 110 as strike you rotate the reference frame to 110. Does that work?

[alkirkby]

Hi @kujaku11 -
Just pulled the latest updates and tested.
Confirming that if I run the following:

mto_unrotated = MT()
mto_unrotated.from_edi('unrotated edi file')
mto_unrotated.coordinate_reference_frame = "ned"
mto_unrotated.rotate(101)

mto_rotated = MT()
mto_rotated.from_edi('edi file rotated to 101 degrees in geotools')

The impedances of the two objects are the same (within 1e-6). Note there is a difference as mto_unrotated.Z.z seems to be rounded to 4-6 d.p.?

The new rotation angle of mto_unrotated after rotation is 101, so that seems correct.

Two comments:

• The rotation angle of the mto_rotated is not read from file (so mto_rotated.rotation_angle returns 0 when really it is 101)
• The code throws an error if I don't set mto_unrotated.coordinate_reference_frame. Would it be better to have a default set so the user doesn't have to set this? If I hadn't read @kkappler's notebook I would have no idea what ned or enu are.

[kujaku11]

@alkirkby Great, at least there is confirmation that we are doing the correct rotation.

1. I'll have a look at the rotation angle read from a file. That is in mt-metadata.

2. The default should be ned, but your probably right there should be some more obvious option. I'll work on the logic for making sure there is a default value always set. Could you post the header of the file you read in and the error? Just curious why it didn't set the default value.

[alkirkby]

Hi @kujaku11 I don't know if I can post the whole header as it's confidential data - was there a field in there that I can check for?

[kujaku11]

@alkirkby no problem. Is there any type of sign convention parameter in the header or notes section?

[kujaku11]

@alkirkby Can you pull the main branch of mt-metadata and confirm that the rotation angles are being pulled from the EDI file into the MT object correctly?

I also changed the logic so that the coordinate_reference_frame is default NED but add other options that may make sense to the user. See the doc string for MT if you want more options. But under the hood there are only 2 NED or ENU. I added a warning if the user sets to ENU cause if you assume that coordinate reference frame then the phase tensor, and plots will likely be incorrect and I'm not going to accommodate ENU at this time.

If that works I'll make this pull request live for you and Karl to review.

[alkirkby]

Hmm, I think I have mt-metadata installed through conda and not git. Is there a branch of mtpy I can pull that will do this? Or should I install mt-metadata and somehow link it to mtpy?