Rotations fixed

[kujaku11]

Rotation of the impedance and tipper were correct in the original versions. Plotting schemes were not and have been patched. All plotting schemes should now be coherent.

Description

Changed all the plotting schemes to now have a rotation_angle property that will rotate the impedance tensor and induction vectors clockwise positive.

Most of the plotting tools, especially polar plots assume positive counterclockwise and 0 = E. Adapted the plotting schemes to adhere to this by 90-strike.

Motivation and Context

Should have made the rotations consistent across the plotting schemes

How Has This Been Tested?

Have a test data set with a conductive square in a half space and made sure the ellipse and induction arrows were oriented properly. Visual tests, someone should verify.

Screenshots (if appropriate):

Types of changes

• [x ] Bug fix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to change)

Checklist:

• [x ] My code follows the code style of this project.
• My change requires a change to the documentation.
• [x ] I have updated the documentation accordingly.
• I have added tests to cover my changes.
• All new and existing tests passed.

[zhang01GA]

Thanks Jared for fixing the rotation issue. I am OK with merging.
Hi @alkirkby , if you have completed the correctness-testing?

[zhang01GA]

I am OK with the merging request if @alkirkby completed the correctness-testing.

[alkirkby]

I'm sorry I haven't yet, will get to this today and then merge if all is good.

[alkirkby]

I get this error when I try to run mtpy/examples/scripts/plot_phase_tensor_maps.py:
Traceback (most recent call last):

File "C:\mtpywin\mtpy\examples\scripts\plot_phase_tensor_map.py", line 15, in
import mtpy.imaging.phase_tensor_maps as pptmaps

File "C:\mtpywin\mtpy\mtpy\imaging\phase_tensor_maps.py", line 26, in
import mtpy.imaging.mtplottools as mtpl

File "C:\mtpywin\mtpy\mtpy\imaging\mtplottools.py", line 17, in
from mtpy.core import mt

File "C:\mtpywin\mtpy\mtpy\core\mt.py", line 20, in
import mtpy.analysis.pt as MTpt

File "C:\mtpywin\mtpy\mtpy\analysis\pt.py", line 28, in
class PhaseTensor(object):

File "C:\mtpywin\mtpy\mtpy\analysis\pt.py", line 328, in PhaseTensor
@z_err.setter

AttributeError: 'function' object has no attribute 'setter'

[alkirkby]

I seem to get this error below quite a bit, this should be fixed before we merge to develop. Let me know if you need the underlying code (I got a similar error when trying to run the phase tensor maps example in the examples/scripts directory

File "C:\mtpywin\mtpy\mtpy\analysis\pt.py", line 28, in
class PhaseTensor(object):

File "C:\mtpywin\mtpy\mtpy\analysis\pt.py", line 328, in PhaseTensor
@z_err.setter

AttributeError: 'function' object has no attribute 'setter'

[kujaku11]

Hi Alison and fei, I fixed the issue, z_err in pt was not decorated as a property. It's now @Property Def z_err() Hopefully that works now. Cheers

…

On Wed, May 20, 2020, 12:08 AM Alison Kirkby ***@***.***> wrote: ***@***.**** requested changes on this pull request. I seem to get this error below quite a bit, this should be fixed before we merge to develop. Let me know if you need the underlying code (I got a similar error when trying to run the phase tensor maps example in the examples/scripts directory File "C:\mtpywin\mtpy\mtpy\analysis\pt.py", line 28, in class PhaseTensor(object): File "C:\mtpywin\mtpy\mtpy\analysis\pt.py", line 328, in PhaseTensor @z_err.setter AttributeError: 'function' object has no attribute 'setter' — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#111 (review)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAQA5BS2JKU7OKC3U6JKCRLRSN6VNANCNFSM4M6MVMRQ> .

[alkirkby]

OK great, that works now. The strike still seems to be repeated in the 1st and 4th quadrants on the strike plots, is this correct?

Will look into the phase tensor plots now.

[alkirkby]

Would you be able to update the examples to show how to use the new rotation_angle property? I'm having trouble using it. I tried setting a rotation_angle property in examples/scripts/plot_phase_tensor_maps.py but it is giving me this error:

File "C:\mtpywin\mtpy\mtpy\imaging\phase_tensor_maps.py", line 518, in init
self.plot()

File "C:\mtpywin\mtpy\mtpy\imaging\phase_tensor_maps.py", line 834, in plot
**self.kwargs)

File "C:\W10DEV\Anaconda3\envs\mtpywin\lib\site-packages\matplotlib\patches.py", line 1394, in init
Patch.init(self, **kwargs)

File "C:\W10DEV\Anaconda3\envs\mtpywin\lib\site-packages\matplotlib\patches.py", line 96, in init
self.update(kwargs)

File "C:\W10DEV\Anaconda3\envs\mtpywin\lib\site-packages\matplotlib\artist.py", line 1006, in update
ret = [_update_property(self, k, v) for k, v in props.items()]

File "C:\W10DEV\Anaconda3\envs\mtpywin\lib\site-packages\matplotlib\artist.py", line 1006, in
ret = [_update_property(self, k, v) for k, v in props.items()]

File "C:\W10DEV\Anaconda3\envs\mtpywin\lib\site-packages\matplotlib\artist.py", line 1002, in _update_property
.format(type(self).name, k))

AttributeError: 'Ellipse' object has no property 'rotation_angle'

[alkirkby]

I am also having trouble using the orthogonal and fold options in the strike roseplots, but that is a separate issue I think (I get the same issue in develop)

[kujaku11]

Sure, hats a good idea. I'll start on that tomorrow.

…

On Sun, May 24, 2020, 4:57 PM Alison Kirkby ***@***.***> wrote: I am also having trouble using the orthogonal and fold options in the strike roseplots, but that is a separate issue I think (I get the same issue in develop) — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#111 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAQA5BSP5FI6O2C3PZNLIH3RTGX53ANCNFSM4M6MVMRQ> .

[alkirkby]

Another thing I noticed is if you plot phase tensors using rotated edis, the phase tensors rotate by the rotation angle of the edi's. I guess this is what the rotation_angle property is for, to rotate it back? We should probably get the code to read in the rotation angle from the edi files and rotate back to zero for plotting.

[kujaku11]

Hi Alison, I think I fixed the bugs you found. I refactored the plot_pt_maps cause there was a lot to read through. I don't have a tif so couldn't check the raster. Bren, can you check that? I added to your examples/scripts/plot_strike_roseplots.py with examples of the various parameters, hopefully that helps. I set plot_orthogonal to False as the default. Good question on the rotation angles from the EDI files. Should the be default action be to rotate back to 0? Not sure, or if it is then there should be a warning that the data have been rotated back to 0. The rotation angle has been a bit of an enigma if it comes from someone else. What is that rotation angle relative to? Geographic, geomagnetic, something else. We just need to be sure that the user knows what happens when the rotation angle in the edi is not zero what happens. If we do things behind the scenes it can get confusing. Cheers

…

On Mon, May 25, 2020 at 3:58 PM Alison Kirkby ***@***.***> wrote: Another thing I noticed is if you plot phase tensors using rotated edis, the phase tensors rotate by the rotation angle of the edi's. I guess this is what the rotation_angle property is for, to rotate it back? We should probably get the code to read in the rotation angle from the edi files and rotate back to zero for plotting. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#111 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAQA5BXM2VJVJWMOGP7AKDTRTLZZTANCNFSM4M6MVMRQ> .

[alkirkby]

Sounds great, I will test it now.
Yeah I think the default should be to rotate back so that X is north and Y is east because this is the coordinate system the ellipses are plotted in. If we leave a rotation of say 30 degrees, the impedance values in the edi are relative to X = 30 degrees east of north so the phase tensor will be shown incorrectly if it's plotted with the X axis at 0 degrees (North). Does that sound right to you or am i missing something?

The information on what north is relative to in the data should be in the edi so can we get that out too? Usually it's geographic north I think. Our plots are generally shown with geographic north vertical, so I guess the data should be rotated so that the north axis in the edi matches north in the plot?

[kujaku11]

I think your correct in that if we plot an edi the orientation should be to geographic north. One issue I foresee is if we make that automatic, what if the user wants to rotate the data for some reason, will the plotting program automatically reset that rotation back to geographic? I guess there should be an option for plot geographic, if true automatically rotate the data to align with geographic. If false plots in data coordinates? Again, I think your assumption that in the edi files 0 is geographic north. I know that some older edi files have arbitrary points of references that are not documented. Mtpy formatted edi files should have that information in there so that's something we can pull. We should just notify the user that we are assuming 0 is geographic north and if something looks fishy check the point of reference?

…

On Tue, May 26, 2020, 2:25 PM Alison Kirkby ***@***.***> wrote: Sounds great, I will test the new branch. Yeah I think the default should be to rotate back so that X is north and Y is east because this is the coordinate system the ellipses are plotted in. If we leave a rotation of say 30 degrees, the impedance values in the edi are relative to X = 30 degrees east of north so the phase tensor will be shown incorrectly if it's plotted with the X axis at 0 degrees (North). Does that sound right to you or am i missing something? The information on what north is relative to in the data should be in the edi so can we get that out too? Usually it's geographic north I think. Our plots are generally shown with geographic north vertical, so I guess the data should be rotated so that the north axis in the edi matches north in the plot? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#111 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAQA5BXTUMYRC4IS5IAYW4TRTQXTXANCNFSM4M6MVMRQ> .

[alkirkby]

Yeah, I guess there should be a plot_geographic option. But regardless of how the data are rotated, I thought the phase tensor azimuth should always be the same relative to north? If the data are rotated by 30 degrees the phase tensors should plot relative to the 30 degree axis, shouldn't they, so they should look the same as unrotated data would look like plotted normally? It doesn't make sense to me that the rotation angle should affect how they plot?

[alkirkby]

I guess if the reference is not recorded in the edi file we have to assume something, and geographic north is probably our best bet in terms of a default.

[alkirkby]

OK almost there with testing. I just tried the rotation_angle property in /examples/scripts/plot_phase_tensor_section.py and when I set the rotation_angle property I get an error:

from mtpy.imaging.phase_tensor_pseudosection import PlotPhaseTensorPseudoSection
import os.path as op
import os

path to edis

edi_path = r'C:\mtpywin\mtpy\examples\data\edi_files_2'

save path

savepath = r'C:\tmp'

edi list

elst=[op.join(edi_path,edi) for edi in os.listdir(edi_path) if ((edi.endswith('.edi')))]# and edi.startswith('GB')

create a plot object

plotObj = PlotPhaseTensorPseudoSection(fn_list = elst,
linedir='ns', # 'ns' if the line is closer to north-south, 'ew' if line is closer to east-west
stretch=(17,8), # determines (x,y) aspect ratio of plot
station_id=(0,10), # indices for showing station names
plot_tipper = 'yri', # plot tipper ('y') + 'ri' means real+imag
font_size=5,
lw=0.5,
rotation_angle=0,
ellipse_dict = {'ellipse_colorby':'skew_seg',# option to colour by phimin, phimax, skew, skew_seg
'ellipse_range':[-12, 12, 3]} # set color limits - default 0,90 for phimin or max,
# # [-12,12] for skew. If plotting skew_seg need to provide
# # 3 numbers, the 3rd indicates interval, e.g. [-12,12,3]
)

update ellipse size (tweak for your dataset)

plotObj.ellipse_size = 2.5

plotObj.plot()

Traceback (most recent call last):

File "C:\Users\u64125\OneDrive - Geoscience Australia\AusLAMP_NSW\mtpy_test\plot_phase_tensor_section_2.py", line 33, in
'ellipse_range':[-12, 12, 3]} # set color limits - default 0,90 for phimin or max,

File "C:\mtpywin\mtpy\mtpy\imaging\phase_tensor_pseudosection.py", line 403, in init
for key, value in kwargs:

ValueError: too many values to unpack (expected 2)

[kujaku11]

Hi Alison, I fixed the bug in the program, easy fix Line 404 in imaging/phase_tensor_pseudosection.py Changed - for key, value in kwargs: to - for key, value in kwargs.items():

…

On Tue, May 26, 2020 at 10:23 PM Alison Kirkby ***@***.***> wrote: OK almost there with testing. I just tried the rotation_angle property in /examples/scripts/plot_phase_tensor_section.py and when I set the rotation_angle property I get a mysterious error: from mtpy.imaging.phase_tensor_pseudosection import PlotPhaseTensorPseudoSection import os.path as op import os path to edis edi_path = r'C:\mtpywin\mtpy\examples\data\edi_files_2' save path savepath = r'C:\tmp' edi list elst=[op.join(edi_path,edi) for edi in os.listdir(edi_path) if ((edi.endswith('.edi')))]# and edi.startswith('GB') create a plot object plotObj = PlotPhaseTensorPseudoSection(fn_list = elst, linedir='ns', # 'ns' if the line is closer to north-south, 'ew' if line is closer to east-west stretch=(17,8), # determines (x,y) aspect ratio of plot station_id=(0,10), # indices for showing station names plot_tipper = 'yri', # plot tipper ('y') + 'ri' means real+imag font_size=5, lw=0.5, rotation_angle=0, ellipse_dict = {'ellipse_colorby':'skew_seg',# option to colour by phimin, phimax, skew, skew_seg 'ellipse_range':[-12, 12, 3]} # set color limits - default 0,90 for phimin or max, # # [-12,12] for skew. If plotting skew_seg need to provide # # 3 numbers, the 3rd indicates interval, e.g. [-12,12,3] ) update ellipse size (tweak for your dataset) plotObj.ellipse_size = 2.5 plotObj.plot() Traceback (most recent call last): File "C:\Users\u64125\OneDrive - Geoscience Australia\AusLAMP_NSW\mtpy_test\plot_phase_tensor_section_2.py", line 33, in 'ellipse_range':[-12, 12, 3]} # set color limits - default 0,90 for phimin or max, File "C:\mtpywin\mtpy\mtpy\imaging\phase_tensor_pseudosection.py", line 403, in *init* for key, value in kwargs: ValueError: too many values to unpack (expected 2) — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#111 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAQA5BXX5G6OMPNPDZHUW4TRTSPVFANCNFSM4M6MVMRQ> .

[kujaku11]

maybe we should handle it as follows: - read in transfer function - If the rotation angle is not zero then notify the user that the data will be rotated to by the -rotation angle to 0 assuming 0 is geographic north - search for a key word that indicates the rotation angle is not relative to geographic north, maybe geomagnetic or principle axis and adjust the azimuth accordingly. If not it is up to the user to figure it out. I don't think there will be many cases where it isn't geographic north. Rotations are confusing and get me every time. As I understand it, the phase tensor azimuth is relative to the coordinate system in which the x and y align to. Commonly those are how the survey was set up, typically geomagnetic. So if you estimate the strike angle from phase tensor of the 'raw' data, the azimuth will be the angle relative to geomagnetic north. If you want to put that into geographic coordinates you need to add on the declination. The ellipse shape should not change though because phimin, phimax, and skew are all invariant. Does that make sense?

…

On Tue, May 26, 2020 at 9:22 PM Alison Kirkby ***@***.***> wrote: Yeah, I guess there should be a plot_geographic option. But regardless of how the data are rotated, I thought the phase tensor azimuth should always be the same relative to north? If the data are rotated by 30 degrees the phase tensors should plot relative to the 30 degree axis, shouldn't they, so they should look the same as unrotated data would look like plotted normally? It doesn't make sense to me that the rotation angle should affect how they plot? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#111 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAQA5BT75CCD6CZYTZCDQA3RTSIR5ANCNFSM4M6MVMRQ> .

[alkirkby]

Cool, yep that sounds like a good approach. If the correct metadata is recorded in the edi's, the datum for the rotation should be there. But yes, that's not always the case.

And yeah, that's my understanding with rotations too. And since all the map plotting functions in MTPy (correct me if I'm wrong) have geographic north as vertical on the plot, we should always "correct" or rotate the data so that the X axis of the impedance tensor is also vertical. If the data have north as geomagnetic north, they need to be rotated back to geographic north before plotting...

So if the user feeds in rotated edi files, we need to "correct" them so that north is vertical. Your approach described above would do that, so I think that's the best way (with plenty of warnings saying that's what we're doing).

[zhang01GA]

Hi @alkirkby and @kujaku11 Recently we have switched to Python-3 auto unit-testing in Github + Travis CI. So could you please merge the updated "develop" into your branch "rotation_test" and then run pytest locally in your machine using Python-3 environment. I have done such a test in my Ubuntu18.04 with Anaconda virtual env Python 3.7.6.

It would be great if you can investigate and fix the unit tests before merge.

Run pytest -v --ignore=tests/SmartMT tests

Observed 4 failed unittests:

=========================== short test summary info ============================
FAILED tests/analysis/test_geometry.py::Test_Geometry::test_get_eccentricity_from_edi_file
FAILED tests/core/test_z.py::TestZ::test_rotate - AssertionError: False is no...
FAILED tests/modeling/test_occam2d.py::TestOccam2D::test_fun - ValueError: ca...
FAILED tests/modeling/ModEM/test_modem_inputfiles_builder.py::TestModemInputFilesBuilder::test_fun_rotate
===== 4 failed, 188 passed, 6 skipped, 3179 warnings in 1498.80s (0:24:58) =====

[kujaku11]

Thanks Fei , I will try that tonight or tomorrow morning.

…

On Sun, May 31, 2020, 5:36 PM Fei Zhang ***@***.***> wrote: Hi @alkirkby <https://github.com/alkirkby> and @kujaku11 <https://github.com/kujaku11> Recently we have switched to Python-3 auto unit-testing in Github + Travis CI. So could you please merge the updated "develop" into your branch "rotation_test" and then run pytest locally in your machine using Python-3 environment. I have done such a test in my Ubuntu18.04 with Anaconda virtual env Python 3.7.6. It would be great if you can investigate and fix the unit tests before merge. Run *pytest -v --ignore=tests/SmartMT tests* Observed 4 failed unittests: =========================== short test summary info ============================ *FAILED tests/analysis/test_geometry.py::Test_Geometry::test_get_eccentricity_from_edi_file FAILED tests/core/test_z.py::TestZ::test_rotate - AssertionError: False is no... FAILED tests/modeling/test_occam2d.py::TestOccam2D::test_fun - ValueError: ca... FAILED tests/modeling/ModEM/test_modem_inputfiles_builder.py::TestModemInputFilesBuilder::test_fun_rotate* ===== 4 failed, 188 passed, 6 skipped, 3179 warnings in 1498.80s (0:24:58) ===== — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#111 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAQA5BVZS336GNHHNTYE4JDRULZXNANCNFSM4M6MVMRQ> .

[kujaku11]

Hi Fei, I pulled from develop and completed the testing. There were 2 errors that I fixed with rotations: - core/test_z had the wrong rotation matrix, it was positive counterclockwise, now it is positive clockwise and passes. - analysis/test_geometry had the wrong values for the eccentricity for the default. Changed those to what they should be and passes. The other errors that I got have to do with the ModuleNotFoundError: No module named 'icons_qt5_rc' and I don't know how to fix that. I pushed to rotation_test. I am working on how to deal with plotting transfer functions that are already rotated, something Alison had mentioned. Cheers On Sun, May 31, 2020 at 7:27 PM Jared Peacock <[email protected]> wrote:

…

Thanks Fei , I will try that tonight or tomorrow morning. On Sun, May 31, 2020, 5:36 PM Fei Zhang ***@***.***> wrote: > Hi @alkirkby <https://github.com/alkirkby> and @kujaku11 > <https://github.com/kujaku11> Recently we have switched to Python-3 auto > unit-testing in Github + Travis CI. So could you please merge the updated > "develop" into your branch "rotation_test" and then run pytest locally in > your machine using Python-3 environment. I have done such a test in my > Ubuntu18.04 with Anaconda virtual env Python 3.7.6. > It would be great if you can investigate and fix the unit tests before > merge. > > Run *pytest -v --ignore=tests/SmartMT tests* > > Observed 4 failed unittests: > > =========================== short test summary info > ============================ > > > > *FAILED > tests/analysis/test_geometry.py::Test_Geometry::test_get_eccentricity_from_edi_file > FAILED tests/core/test_z.py::TestZ::test_rotate - AssertionError: False is > no... FAILED tests/modeling/test_occam2d.py::TestOccam2D::test_fun - > ValueError: ca... FAILED > tests/modeling/ModEM/test_modem_inputfiles_builder.py::TestModemInputFilesBuilder::test_fun_rotate* > ===== 4 failed, 188 passed, 6 skipped, 3179 warnings in 1498.80s > (0:24:58) ===== > > — > You are receiving this because you were mentioned. > Reply to this email directly, view it on GitHub > <#111 (comment)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AAQA5BVZS336GNHHNTYE4JDRULZXNANCNFSM4M6MVMRQ> > . >