MTData: Grid Example

It is becoming more common for MT data to be collected in a grid rather than a profile because of 3D inversions. In this example we can take a look at how to work with MT data collected on a grid. We will use the same MTCollection from the first example.

1. Load data

First we will open the MTCollection and change the working dataframe to get only those station in the ‘grid’ survey.

[1]:

from pathlib import Path
from mtpy import MTCollection

[2]:

with MTCollection() as mtc:
    mtc.open_collection(Path().cwd().joinpath("test_mt_collection.h5"))
    mtc.working_dataframe = mtc.master_dataframe.loc[mtc.master_dataframe.survey == "grid"]
    mtd = mtc.mt_data

26:05:08T16:56:50 | INFO | line:1035 |mth5.mth5 | close_mth5 | Flushing and closing c:\Users\peaco\OneDrive\Documents\GitHub\mtpy-v2\docs\source\notebooks\test_mt_collection.h5

2. Plot Station Locations

To make sure we got the data we expected, we can plot the station locations. You can change the basemap, see providers for more details

[3]:

station_plot = mtd.plot_stations()
../_images/notebooks_3_MTData_grid_5_0.png

2a. Change basemap

Here is an example of how to change the basemap. We it to use the ESRI terrain map.

[4]:

import contextily as cx

[5]:

station_plot.cx_source = cx.providers.Esri.NatGeoWorldMap
station_plot.redraw_plot()
../_images/notebooks_3_MTData_grid_8_0.png

3. Plot Phase Tensor Map

Now that we seem to have the correct data, lets plot phase tensor maps. These are broadband data with induction vectors, so lets plot those too. The phase tensor map is busy so lets identify key aspects.

  • Ellipse shape: the ellipses often elongate in the preferred direction of current flow

  • Ellipse face color: by default are colored by the parameter phimin but can be changed by setting the attribute pt_map.ellipse_colorby. The phimin gives the lower bounds on how the subsurface resistivity is changing, where reds are becoming more conductive and blues are becoming more resistive.

  • Ellipse edge color: by default is colored by the skew angle which is indicative of dimensionality with high skew being 3D. Also the color indicates in which direction currents are being skewed.

  • Wedges: the long axis is phimax and the short axis are phimin.

  • Arrows: black are real induction vectors and blue are imaginary induction vectors.

[6]:

pt_map = mtd.plot_phase_tensor_map(
    plot_tipper="yri",
    cx_source=cx.providers.Esri.NatGeoWorldMap,
    ellipse_size=.02,
    arrow_size=.05
)
../_images/notebooks_3_MTData_grid_10_0.png 
[7]:

pt_map.plot_period = 100
pt_map.redraw_plot()
../_images/notebooks_3_MTData_grid_11_0.png

4. Plot Resistivity and Phase Maps

An informative first order check on the data is to plot apparent resistivity and phase maps of the data. This can help identify odd stations for further investigations and get a general idea of subsurface resistivity changes. These plots are a little messy in a Jupyter Notebook.

[8]:

rp_map = mtd.plot_resistivity_phase_maps(plot_det=True, subplot_wspace=.45, marker_size=4)
../_images/notebooks_3_MTData_grid_13_0.png

4a. Plot Diagonal Components

It can be informative to plot the diagnoal components as well, in this case they outline the basins nicely.

[9]:

rp_map.plot_xx = True
rp_map.plot_yy = True
rp_map.plot_det = False
rp_map.redraw_plot()
../_images/notebooks_3_MTData_grid_15_0.png

5. Plot Depth of Investigation

It can also be informative to understand how deep the measurements are sensitive to. Here a Niblett-Bostick approximation is used to estimate the depth of penetration for each station at a single period.

[10]:

depth_of_penetration = mtd.plot_penetration_depth_map(subplot_wspace=.35)
../_images/notebooks_3_MTData_grid_17_0.png

6. Plot Strike

Here we will plot all periods of estimated strike. Notice that the plot includes the strike as estimated from the invariants (left) of Weaver et al. (2002), the phase tensor (middel) of Caldwell et al. (2004), and the induction vector strike.

Important: The induction strike points towards good conductors so should therefore be perpendicular to the impedance strike. We left it this way as a sanity check on strike angles.

6a. Plot all periods in one plot

This plot shows all strike estimations for all stations for all periods.

[11]:

strike_plot = mtd.plot_strike()

26:05:08T16:58:52 | INFO | line:860 |mtpy.imaging.plot_strike | _plot_all_periods | Note: North is assumed to be 0 and the strike angle is measured clockwise positive.
../_images/notebooks_3_MTData_grid_19_1.png

6b. Strike per period

[12]:

strike_plot.plot_type = 1
strike_plot.redraw_plot()

26:05:08T16:58:54 | INFO | line:765 |mtpy.imaging.plot_strike | _plot_per_period | Note: North is assumed to be 0 and the strike angle is measuredclockwise positive.
../_images/notebooks_3_MTData_grid_21_1.png 
[13]:

strike_plot.plot_orientation = "vertical"
strike_plot.redraw_plot()

26:05:08T16:58:57 | INFO | line:765 |mtpy.imaging.plot_strike | _plot_per_period | Note: North is assumed to be 0 and the strike angle is measuredclockwise positive.
../_images/notebooks_3_MTData_grid_22_1.png

7. Extract a Profile

One adventageous way of looking at grid data is to extract profiles across interesting structures. This can be done by knowing the start and end points of the profile you would like to look at and providing a distance away from the profile line to include.

By analyzing the strike direction we can see that a profile approximately N75E would be perpendicular to geoelectric strike. Lets cut across the basins.

[14]:

mtd.utm_crs = 32611

[15]:

%%time
profile = mtd.get_profile(-118.45, 38.78, -117.85, 38.95, 5000)

CPU times: total: 1.75 s
Wall time: 1.77 s

[16]:

profile_stations_plot = profile.plot_stations()
../_images/notebooks_3_MTData_grid_26_0.png 
[ ]: