kujaku11 commented 4 years ago

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 commented 4 years ago

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

alkirkby commented 4 years ago

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

alkirkby commented 4 years ago

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'

kujaku11 commented 4 years ago

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 notifications@github.com wrote:

@alkirkby 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 https://github.com/MTgeophysics/mtpy/pull/111#pullrequestreview-415047836, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAQA5BS2JKU7OKC3U6JKCRLRSN6VNANCNFSM4M6MVMRQ .

alkirkby commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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

On Sun, May 24, 2020, 4:57 PM Alison Kirkby notifications@github.com 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 https://github.com/MTgeophysics/mtpy/pull/111#issuecomment-633319088, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAQA5BSP5FI6O2C3PZNLIH3RTGX53ANCNFSM4M6MVMRQ .

alkirkby commented 4 years ago

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 commented 4 years ago

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 notifications@github.com 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 https://github.com/MTgeophysics/mtpy/pull/111#issuecomment-633736939, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAQA5BXM2VJVJWMOGP7AKDTRTLZZTANCNFSM4M6MVMRQ .

alkirkby commented 4 years ago

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 commented 4 years ago

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 notifications@github.com 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 https://github.com/MTgeophysics/mtpy/pull/111#issuecomment-634287718, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAQA5BXTUMYRC4IS5IAYW4TRTQXTXANCNFSM4M6MVMRQ .

alkirkby commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 notifications@github.com 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 https://github.com/MTgeophysics/mtpy/pull/111#issuecomment-634435786, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAQA5BXX5G6OMPNPDZHUW4TRTSPVFANCNFSM4M6MVMRQ .

kujaku11 commented 4 years ago

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 notifications@github.com 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 https://github.com/MTgeophysics/mtpy/pull/111#issuecomment-634419369, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAQA5BT75CCD6CZYTZCDQA3RTSIR5ANCNFSM4M6MVMRQ .

alkirkby commented 4 years ago

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 commented 4 years ago

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_funrotate ===== 4 failed, 188 passed, 6 skipped, 3179 warnings in 1498.80s (0:24:58) =====

kujaku11 commented 4 years ago

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

On Sun, May 31, 2020, 5:36 PM Fei Zhang notifications@github.com 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 https://github.com/MTgeophysics/mtpy/pull/111#issuecomment-636556694, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAQA5BVZS336GNHHNTYE4JDRULZXNANCNFSM4M6MVMRQ .

kujaku11 commented 4 years ago

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 peacock.jared@gmail.com wrote:

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

On Sun, May 31, 2020, 5:36 PM Fei Zhang notifications@github.com 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 https://github.com/MTgeophysics/mtpy/pull/111#issuecomment-636556694, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAQA5BVZS336GNHHNTYE4JDRULZXNANCNFSM4M6MVMRQ .

MTgeophysics / mtpy

Rotations fixed #111

Description

Motivation and Context

How Has This Been Tested?

Screenshots (if appropriate):

Types of changes

Checklist:

path to edis

save path

edi list

create a plot object

[-12,12] for skew. If plotting skew_seg need to provide

update ellipse size (tweak for your dataset)

[-12,12] for skew. If plotting skew_seg need to provide

3 numbers, the 3rd indicates interval, e.g. [-12,12,3]

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

=========================== short test summary info

=========================== short test summary info