MTpy Example 08 - Open-source Python packages for Magnetotellurics

Geoelectric Strike¶

Geoelectric strike indicates the preferred direction of electrical current flow. The strike direction becomes important when you want to figure out which modes of the impedance tensor are which, namely the TE and TM modes. TE is parallel to geoelectric strike and TM is perpendicular. When modeling in 1D and 2D orienting your modes properly is important. Modeling in 3D removes this dependency because all 4 components can be modeled.

MTpy provides 3 ways to estimate strike direction, through the impedance tensor invariants Liley, 1998, the phase tensor azimuth Caldwell et al., 2004, and the induction vector real azimuth.

The plots are rose diagrams where up is North (0 degrees strike) and East is the right (90 degrees strike). Therefore, strike is estimated as positive clockwise. The strike angle is relative to the coordinate system of the transfer functions, in this case geographic North.

You can plot the strike angles as a function of period (by decade), or you can use the working_dataframe to organize the strike angle by geographic region. We will show both.

NOTE: the induction vector or tipper angles should be perpendicular to the strike because they point towards a good conductor. I left it this way because how perpendicular the tippers can be a good indication of dimensionality.

from pathlib import Path
from mtpy import MTCollection
%matplotlib inline

Open MTCollection¶

In the previous example we created a MTH5 file from existing Yellowstone data. Let’s open that file here for plotting.

with MTCollection() as mc:
    mc.open_collection(Path().cwd().parent.parent.joinpath("data", "transfer_functions", "yellowstone_mt_collection_03.h5"))
    mt_data = mc.to_mt_data()

Plot Strike of All Stations¶

We have a large survey, so this may not be very useful as far as analysis, but useful to know we can do it.

All Periods Together¶

Here we are plotting the strike estimation for all station and all periods into one rose diagram. For an area this large, this plot is more or less useless. But good to see.

plot_strike_all = mt_data.plot_strike()

24:10:08T16:04:55 | INFO | line:861 |mtpy.imaging.plot_strike | _plot_all_periods | Note: North is assumed to be 0 and the strike angle is measured clockwise positive.

Plot Per Decade¶

A more useful way to plot the strike angle is per period. Again this is a large survey so these don’t have too much useful information. The numbers in the boxes under the plot are the median strike values for that decade.

plot_strike_all.plot_type = 1
plot_strike_all.fig_num = 2
# to fit on a laptop screen
plot_strike_all.font_size = 3
plot_strike_all.redraw_plot()

24:10:08T16:04:58 | INFO | line:765 |mtpy.imaging.plot_strike | _plot_per_period | Note: North is assumed to be 0 and the strike angle is measuredclockwise positive.

Plot DeGroot-Hedlin Line¶

The DeGroot-Hedlin line across the Snake River Plane is a good example to strike estimation. The profile line is aligned northwest, which was picked to be approximately perpendicular to the geoelectric strike and key geologic structures.

At most periods the geoelectric strike is roughly northeast with maybe 2 similar dominant strike angle, but notice the rotation between 100 - 1000 seconds, this could be near the Moho where there is a change from brittle to ductile material.

%%time
mc.working_dataframe = mc.master_dataframe[mc.master_dataframe.station.str.startswith("SR2")]
dgh_mt_data = mc.to_mt_data()

CPU times: total: 12.2 s
Wall time: 13.1 s

plot_strike_dgh = dgh_mt_data.plot_strike(plot_type=1, fig_num=4, font_size=3)

24:10:08T16:05:22 | INFO | line:765 |mtpy.imaging.plot_strike | _plot_per_period | Note: North is assumed to be 0 and the strike angle is measuredclockwise positive.

Rotate to Geoelectric Strike¶

It can be informative to understand how homogeneous or hetergeneous the geoelectric strike angle is for each period. If we rotate to the dominant strike direction N33E then will see how the strike angle varies with period. We will rotate by N33W so that N33E is up on the plot. The plot demonstrates that a change in geoelectric strike is observed between 100 - 1000 seconds.

plot_strike_dgh.rotation_angle = 30
plot_strike_dgh.fig_num = 6
plot_strike_dgh.redraw_plot()

24:10:08T16:06:57 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:06:57 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:06:57 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:06:58 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:06:58 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:06:58 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:06:58 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:06:58 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:06:58 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:06:59 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:06:59 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:06:59 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:06:59 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:06:59 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:00 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:00 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:00 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:00 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:00 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:00 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:01 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:01 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:01 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:01 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:01 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:02 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:02 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:02 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:02 | INFO | line:219 |mtpy.core.mt | rotate | Rotated transfer function by: 30.000 degrees clockwise in reference frame NED.
24:10:08T16:07:08 | INFO | line:765 |mtpy.imaging.plot_strike | _plot_per_period | Note: North is assumed to be 0 and the strike angle is measuredclockwise positive.

Plot by Geologic Region¶

We can estimate various geologic regions from the map. Looks like:

Stations "SR205 to “SR232” are in mountains to the north
Stations “SR232” to “SR255” are in the Snake River Plane
Stations “SR255” to “SR290” are in mountains to the south

plot_stations = dgh_mt_data.plot_stations(fig_num=8)

24:10:08T16:05:40 | WARNING | line:163 |mtpy.imaging.plot_stations | plot | Could not add base map because HTTPSConnectionPool(host='basemap.nationalmap.gov', port=443): Max retries exceeded with url: /arcgis/rest/services/USGSTopo/MapServer/tile/8/92/46 (Caused by SSLError(SSLCertVerificationError(1, '[SSL: CERTIFICATE_VERIFY_FAILED] certificate verify failed: self-signed certificate in certificate chain (_ssl.c:1129)')))

# Need to make this from the master and make it a list so that it doesn't get over written 
# when we set the working dataframe
dgh_stations = list(mc.master_dataframe[mc.master_dataframe.station.str.startswith("SR2")].station.unique())

mc.master_dataframe.survey.unique()

array(['Transportable_Array', 'YSBB', 'Yellowstone-Snake_River_Plain'],
      dtype=object)

Northern Stations¶

northern_stations = [f"Yellowstone-Snake_River_Plain.{station}" for station in dgh_stations if int(station[2:]) < 233]
northern_mt_data = dgh_mt_data.get_subset(northern_stations)

plot_northern_strike = northern_mt_data.plot_strike(fig_num=9, plot_type=1)

24:10:08T16:10:02 | INFO | line:765 |mtpy.imaging.plot_strike | _plot_per_period | Note: North is assumed to be 0 and the strike angle is measuredclockwise positive.

Snake River Plane Stations¶

snp_stations = [f"Yellowstone-Snake_River_Plain.{station}" for station in dgh_stations if int(station[2:]) > 231 and int(station[2:]) < 259]
snp_mt_data = dgh_mt_data.get_subset(northern_stations)

plot_snp_strike = snp_mt_data.plot_strike(fig_num=11, plot_type=1)

24:10:08T16:10:35 | INFO | line:765 |mtpy.imaging.plot_strike | _plot_per_period | Note: North is assumed to be 0 and the strike angle is measuredclockwise positive.

Southern Stations¶

southern_stations = [f"Yellowstone-Snake_River_Plain.{station}" for station in dgh_stations if int(station[2:]) > 257]
southern_mt_data = dgh_mt_data.get_subset(southern_stations)

plot_southern_strike = southern_mt_data.plot_strike(fig_num=12, plot_type=1)

24:10:08T16:10:51 | INFO | line:765 |mtpy.imaging.plot_strike | _plot_per_period | Note: North is assumed to be 0 and the strike angle is measuredclockwise positive.

References¶

Lilley, F. E. M. (Ted). (1998). Magnetotelluric tensor decomposition: Part I, Theory for a basic procedure. GEOPHYSICS, 63(6), 1885–1897. 10.1190/1.1444481
Caldwell, T. G., Bibby, H. M., & Brown, C. (2004). The magnetotelluric phase tensor. Geophysical Journal International, 158(2), 457–469. 10.1111/j.1365-246x.2004.02281.x