<jupyter_start><jupyter_text>Inversion of a Transect of MT Tipper data. Part IThis dataset consists of 29 Magnetotelluric stations arranged on a 2D transect. Our objective is to invert the tipper data on line 5000 using MtPy to (hopefully) reveal a conductive structure known to exist under the array. While the task is straightforward, there are a number of smaller tasks we'll have to complete in sequence:1) Loading and cleaning the data2) Parsing data into MtPy compatible formats3) Visualizing the data4) Set up Occam2D Inversion5) Explore Inversion Parameters6) Analyze ResultsBefore following this jupyter notebook tutorial, ensure you have the latest version of miniconda3 installed. Create an ```mtpy_inversion``` environment by installing the following packages in series```conda create --name mtpy_inversionsource activate mtpy_inversionconda install python=3.6.8conda install numpy scipy matplotlib pandas cartopypip install mtpy``` How to use this notebookCells in this notebook are divided into two types:- Helper Cells- Data Manipulation/Exploration CellsOnly the "Helper cells" are labeled for convenience. You may, of course, edit any cell you wish, or copy bits of code to use in your own analysis. There is one external module you need called ```cartopy_mapping.py```. Make sure this python file is in the same directory you run this notebook from.This notebook performs steps 1) and 2) as described above. It is meant to be used prior to executing Part_2.Data_Analysis. The first cell you execute will shoot out some warnings. Don't worry, this won't affect anything.After you're done with this notebook move on to Part II. Step 1. load the data. The data are provided in a .csv file. Since pandas is well-suited for loading csv's, I'll load it into a pandas dataframe.<jupyter_code>"""
HELPER CELL

for importing relevant modules

"""

import matplotlib.cm as cm
import matplotlib.pyplot as plt
import pandas as pd
from scipy import ndimage as nd
import cartopy_mapping
import numpy as np
from mtpy.core.z import Tipper, Z
from mtpy.core.mt import MT
import os
import glob
data = pd.read_csv('original_data/YGS-20200107-ELF_with_repeat_readings.csv')<jupyter_output><empty_output><jupyter_text>Let's Look at the data dimensions<jupyter_code>data<jupyter_output><empty_output><jupyter_text>I know the .csv data is structured Like so:* **Line**: Local grid's line value* **Station**: Local grid's station value* **Lat**: Latitude* **Longitude**: Longitude* **E_WGS84_Z7**: Easting* **N_WGS84_Z7z**: Northing* **ALT**: Station Altitude (m)* **DUR**: Duration of the reading (s)* **Noise**: Level of Industrial Noise* **nE_M**: East Tipper measurement of frequency M. if ```n=='r'```, data is the real component. If ```n=='i'```, data is the imaginary component.* **nN_M**: North Tipper measurement of frequency M. if ```n=='r'```, data is the real component. If ```n=='i'```, data is the imaginary component.* **SNI_M**: Signal-to-noise index of component M.**M** takes on the values of 11, 22, 45, 90, 180, 360, 720, and 1440 Hz.I will further assume the tipper is defined as:$$H_z(\omega) = \bf{[}T_n(\omega), T_e(\omega)\bf{]}\begin{bmatrix}           H_n(\omega) \\           H_e(\omega) \\         \end{bmatrix}$$         And individual components can be expressed as$$\frac{H_z(\omega)}{H_n(\omega)} = T_n(\omega)$$$$\frac{H_z(\omega)}{H_e(\omega)} = T_e(\omega)$$It appears as if each row represents the $T_x(\omega)$ data and metadata relevant to a single station. Lets check how many stations exist in line 5000.<jupyter_code>data_line_5000=data[data['Line']==5000]
len(data_line_5000.Station.unique())<jupyter_output><empty_output><jupyter_text>How many repeat measurements can we find in this line?<jupyter_code>ax = data_line_5000['Station'].value_counts().plot(kind='bar',
                                    figsize=(14,8),
                                    title="Number of measurements per station")
ax.set_xlabel("Station")
ax.set_ylabel("Count")
ax.set_yticks([0,1,2])<jupyter_output><empty_output><jupyter_text>Stations 0, 3000,5000,6000,2000,1000,7000, and 4000 have 2 repeat occupations. All the other lines appear to only have one occupation. For now lets drop the first ocurrence of each repeat occupation.<jupyter_code>df_1            = data_line_5000.sort_values(by=['Station'],ascending = [True])
df_1['cumcount']= df_1.groupby('Station').cumcount()
data_line_5000  = data_line_5000[df_1['cumcount'] ==0]
data_line_5000  = data_line_5000.sort_values(by=['Station'],ascending = [True])
ax = data_line_5000['Station'].value_counts().plot(kind='bar',
                                    figsize=(14,8),
                                    title="Number of measurements per station")
ax.set_xlabel("Station")
ax.set_ylabel("Count")
ax.set_yticks([0,1,2])<jupyter_output><empty_output><jupyter_text>I'm not sure why the stations aren't sorted in ascending order, but our station-pruning operation has been successful.While we have the metadata nicely packed away in a dataframe, lets see how the stations are distributed<jupyter_code>cartopy_mapping.plot_map(data_line_5000)<jupyter_output><empty_output><jupyter_text>The data are arranged on a nice NE-SW transect. How about the data itself?<jupyter_code>"""
HELPER CELL

for plotting the data, separating by both direction and real/imaginary component.

"""
def sort_data_into_columns(df,direction='N',part='real'):
    data_columns = ['rE_0011', 'rE_0022', 'rE_0045', 'rE_0090',
       'rE_0180', 'rE_0360', 'rE_0720', 'rE_1440', 'rN_0011', 'rN_0022',
       'rN_0045', 'rN_0090', 'rN_0180', 'rN_0360', 'rN_0720', 'rN_1440','iE_0011', 'iE_0022', 'iE_0045', 'iE_0090', 'iE_0180', 'iE_0360',
       'iE_0720', 'iE_1440', 'iN_0011', 'iN_0022', 'iN_0045', 'iN_0090',
       'iN_0180', 'iN_0360', 'iN_0720', 'iN_1440']
    cols_to_use =[]
    for col in data_columns:
        if direction in col:
            if part=='real':
                if 'r' in col:
                    cols_to_use.append(col)
            else:
                if 'i' in col:
                    cols_to_use.append(col)
                    
    d_subset = df.loc[:,cols_to_use]
    d_subset = d_subset.reindex(sorted(d_subset.columns), axis=1)
    return d_subset

def plot_data(df,direction='N',part='real'):
    d_subset     = sort_data_into_columns(df,direction=direction,part=part)
    frequencies  = [int(x.split('_')[1]) for x in d_subset.columns]
    stations     = df['Station'].values
    data         = d_subset.values

    fig,ax = plt.subplots(figsize=(10,10))
    current_cmap = cm.get_cmap()
    current_cmap.set_bad(color='red')
    im = ax.imshow(data.T)
    cbar = ax.figure.colorbar(im, ax=ax,fraction=0.012, pad=0.04)
    cbar.ax.set_ylabel('value', rotation=-90, va="bottom")
    ax.set_xticks(np.arange(len(stations)))
    ax.set_yticks(np.arange(len(frequencies)))
    ax.set_xticklabels(stations)
    ax.set_yticklabels(frequencies)
    plt.setp(ax.get_xticklabels(), rotation=45, ha="right",
         rotation_mode="anchor")
    ax.set_title("{} Value of component {}".format(part,direction))
    fig.tight_layout()
    plt.show()
plot_data(data_line_5000,direction='N',part='real')
plot_data(data_line_5000,direction='N',part='imag')
plot_data(data_line_5000,direction='E',part='real')
plot_data(data_line_5000,direction='E',part='imag')<jupyter_output><empty_output><jupyter_text>Oops, looks like there are some ```NAN``` fields which show up as red. Strangely, there are specific frequencies which are ```NAN``` in the real component, but exist in the imaginary component (and vice versa). I'll have to make non-```NAN``` values ```NAN``` if the corresponding component is ```NAN```.<jupyter_code>"""
HELPER CELL

for assigning NAN values to the appropriate cell.

"""
def assign_NAN_values(df,direction='N',part='real'):
    frequencies = ['0011','0022','0045','0090','0180','0360','0720','1440']
    for freq in frequencies:
        for direction in ['N','E']:
            real_key = 'r'+direction+'_'+freq
            imag_key = 'i'+direction+'_'+freq
            df.loc[pd.isna(df[real_key]) | pd.isna(df[imag_key]), real_key] = np.nan
            df.loc[pd.isna(df[real_key]) | pd.isna(df[imag_key]), imag_key] = np.nan
                    
    return df
data_line_5000=assign_NAN_values(data_line_5000)
plot_data(data_line_5000,direction='N',part='real')
plot_data(data_line_5000,direction='N',part='imag')
plot_data(data_line_5000,direction='E',part='real')
plot_data(data_line_5000,direction='E',part='imag')<jupyter_output><empty_output><jupyter_text>Ahh much better. At least the data gaps make sense now. We now have a choice to make. I can either:* Ensure my inversion scheme can handle ```nan``` values.* Drop stations with ```nan``` values.* Restrict my frequency space to only contain non ```nan``` numbers.* Interpolate values within reason.There are drawbacks to each approach. Ensuring your inversion scheme can handle NaN's requires an intimate understanding of the program. My intended approach uses Occam2D, which does not handle nan values well. Dropping all stations with NaN values leaves us with a really sparse dataset. Restricting frequency space to non nan values means we now can only invert 3 of 8 frequencies. Interpolating values is like creating 'fake' data and will bias the inversion results.How about a compromise?First, I'll drop stations 1250, 1750, 4250, 5250, and 6500. These stations have the most nan values, and based on discussions with the team that acquired this data, nan values mean high uncertainty.<jupyter_code>bad_stations = [1250,1750,4250,5250,6500]
new_data  = data_line_5000[~data_line_5000.Station.isin(bad_stations)]
new_data  = new_data.reset_index()
plot_data(new_data,direction='N',part='real')
plot_data(new_data,direction='N',part='imag')
plot_data(new_data,direction='E',part='real')
plot_data(new_data,direction='E',part='imag')<jupyter_output><empty_output><jupyter_text>We still have NaN values on some of the channels. One approach to deal with this is to locally interpolate from the nearest value on each station. We could do this here OR we could do this when we convert to MtPy. I prefer the second option as it will allow us to use pandas ```fillna``` methods.  Step 2. Prepare Data for MtPyEvery row in this pandas dataframe represents spectral Tipper data collected from a station. Before we can do any inversion, it's prudent to convert these to xml files on disk. While this isn't necessary, it will allow us to separate data cleaning from data analysis.MtPy uses an object called ```MT``` contained in ```mtpy.core.mt```. An easy way to do the conversion is to use a function to assign every row's ```MT``` Object into an ```mt_obj``` column.Documentation for ```MtPy's``` ```MT```object is lacking, but from reading the git repository I can tell the following items should be incorporated into our object:* **Metadata**: Stuff like lat, lon, elevation, station. These are simple float/string attributes. Note that we can either provide MtPy with lat,lon or easting-northing information. I prefer using lat-lons in WGS84 format. It's a universal format used by most WMS, google maps, google earth, navigation software, and smartphone GPS apps. While it might be hard to convert lat-lon into meaningful distance measurements in the field, with the availability of packages like ```pyproj```, I find sending UTM coordinates are at best redundant and at worst, add unneeded complexity to data processing.* **Tipper Data**: This is a bit more complex. It requires an array of frequencies, complex numbers, and needs to be shaped appropriatey. There is an option for incorporating error into the Tipper array, but for now we'll just set those to zero. **THIS IS WHERE THE FILL_NA COMES IN**. I'll use a fillna operation that relies on a succinct stackoverflow solution using a scipy image processing routine.For ease of understanding, I'll make these separate functions which are then called by ```parse_row_into_mt_obj```.<jupyter_code>"""
HELPER CELL

for parsing pandas dataframe row data into an MT object

"""

def fill(data, invalid=None):
    """
    FROM STACKOVERFLOW: 
    https://stackoverflow.com/questions/3662361/fill-in-missing-values-with-nearest-neighbour-in-python-numpy-masked-arrays
    Replace the value of invalid 'data' cells (indicated by 'invalid') 
    by the value of the nearest valid data cell

    Input:
        data:    numpy array of any dimension
        invalid: a binary array of same shape as 'data'. True cells set where data
                 value should be replaced.
                 If None (default), use: invalid  = np.isnan(data)

    Output: 
        Return a filled array. 
    """

    if invalid is None: invalid = np.isnan(data)

    ind = nd.distance_transform_edt(invalid, return_distances=False, return_indices=True)
    return data[tuple(ind)]

def create_simple_metadata(row):
    """
    extracts metadata from the pandas dataframe row, specifically:
    site (lat, lon, elevation, station, line)
    acquisition (azimuth, recording time)
    """
    lat     = row['Lat']
    lon     = row['Long']
    elev    = row['ALT']
    station = row['Station']
    
    return lat, lon, elev, station

def generate_tipper_arrays(row):
    """
    function to populate numpy arrays from pandas dataframe
    """
    frequencies = np.asarray([11, 22, 45, 90, 180, 360, 720, 1440])
    t_n=[]
    t_e=[]
    for freq in frequencies:
        real_east_key = f"rE_{freq:04}"
        imag_east_key = f"iE_{freq:04}"
        real_north_key= f"rN_{freq:04}"
        imag_north_key= f"iN_{freq:04}"
        t_e.append( row[real_east_key]+row[imag_east_key]*1j)
        t_n.append( row[real_north_key]+row[imag_north_key]*1j)
        
    t_n = fill(np.asarray(t_n,dtype=np.complex64))
    t_e = fill(np.asarray(t_e,dtype=np.complex64))
    return frequencies, t_n, t_e

def create_tipper(row):
    """
    creates an MtPy Tipper array from source numpy arrays and relevant metadata
    
    ############################# BIG WARNING #################################
    Notice Ive done something funny with the error array. This is purely for MtPy compatibility.
    To incorporate error please assign other values.
    """
    frequencies, t_n, t_e = generate_tipper_arrays(row)
    
    array            = np.ones((len(frequencies),1,2),dtype=np.complex64)
    tipper_err_array = np.ones((len(frequencies),1,2)) #<- change this to include error
    array[:,0,0]=t_n
    array[:,0,1]=t_e
    
    tipper = Tipper(tipper_array=array, tipper_err_array=tipper_err_array,
                    freq=frequencies)
    return tipper

def create_dummy_impedanc_tensor():
    """
    creates an MtPy Impedance Tensor object from dummy data
    
    ############################# BIG WARNING #################################
    Both the values and the error are made up here because the data is only made up of tipper data.
    This is done to maintain compatibility with MtPY
    """
    frequencies = np.asarray([11, 22, 45, 90, 180, 360, 720, 1440])
    z_array     = np.ones((len(frequencies),2,2),dtype=np.complex64)
    impedance_tensor = Z(z_array=z_array, z_err_array=z_array, freq=frequencies)
    return impedance_tensor

def parse_row_into_mt_obj(row):
    lat, lon, elev, station = create_simple_metadata(row)
    tipper   = create_tipper(row)
    impedance_tensor  = create_dummy_impedanc_tensor()
    mt_obj   = MT(station=station,lat=lat,lon=lon,elev=elev,Tipper=tipper,Z=impedance_tensor)
    return mt_obj
new_data['mt_obj']=new_data.apply(parse_row_into_mt_obj, axis=1)
new_data<jupyter_output>/home/kevin/anaconda3/envs/mtpy_inv/lib/python3.6/site-packages/mtpy/core/z.py:125: ComplexWarning: Casting complex values to real discards the imaginary part
  self._resistivity[idx_f, ii, jj] * r_err<jupyter_text>One of the under appreciated aspects of pandas is that the live memory version of dataframes can store and index arbitrarily defined objects. I would not reccomend saving this again as a .csv to file, but it allows us to efficiently query, plot, and examine objects. Now that we seem to have created indexed ```MTPy``` objects, we need to verify that the attributes contained in the live objects match those from the dataframe. Trust, but verify. Another benefit of this approach is that I now have  a separate code section specifically dedicated to testin. If errors in the data cleaning are discovered at a later date, I can add additional checks here to prevent the same error from ocurring again.<jupyter_code>"""
HELPER CELL

for assuring the created MT object is equivalent to its dataframe source data

"""
def metadata_assertions(row):
    mt_obj = row['mt_obj']
    lat     = row['Lat']
    lon     = row['Long']
    elev    = row['ALT']
    station = row['Station']
    
    assert mt_obj.elev ==elev,"elevations don't match"
    assert mt_obj.station==station,"stations don't match"
    assert mt_obj.lon==lon,"longitudes don't match"
    assert mt_obj.lat==lat,"latitudes don't match"

def tipper_assertions(row):
    mt_obj = row['mt_obj']
    frequencies, t_n, t_e=generate_tipper_arrays(row)
    source_tipper = mt_obj.Tipper
    target_tipper = create_tipper(row)
    np.testing.assert_array_equal(source_tipper.tipper,target_tipper.tipper,'tipper arrays not equal')
    
def verify_row_corresponds_to_mt_obj(row):
    tipper_assertions(row)
    metadata_assertions(row)
for index, row in new_data.iterrows():
    verify_row_corresponds_to_mt_obj(row)
print('data is verified')<jupyter_output>data is verified<jupyter_text>MtPy supports three main file types often used in MT: .xml, .edi, and .Jfile. While the .edi format is probably the most popular at the moment, a wide variety of perspectives on where and how to store metadata within the format has motivated the community to explore other ways of storing data. My favorite candidate file format is the netcdf4-5 datatype widely used by the meteorological community to describe multidimensional datasets. The .xml format is also an attractive candidate due to its self-describing structure. For now, I'll stick to saving the data in an intermediate .edi format.I'll add some machinery to create and delete any intermediate files in the .edi directory. This is a best practice to ensure that any modification to data preprocessing is both reproducible and robust.The last thing we need to do here is save the ```MT_obj```'s to file and ensure their integrity by loading them up and verifying the data much like we did before.<jupyter_code>directory='edi_files'
# create temporary data directory
if not os.path.exists(directory):
    os.makedirs(directory)
    
# clear existing edi files in directory
files = glob.glob('edi_files/*.edi')
for f in files:
    os.remove(f)
# save every station as an .edi file into directory  
for index, row in new_data.iterrows():
    row['mt_obj'].write_mt_file(save_dir=directory, fn_basename=str(row['Station']), file_type='edi',
                      longitude_format='LON',latlon_format='dd')
"""
HELPER CELL

loads .edi files into a dataframe, then compares it to the previously loaded file

"""

def check_one_edi_per_station():
    files = glob.glob('edi_files/*.edi')
    assert len(files)==24,"differing number of .edi files and stations"
    
def convert_tipper_to_csv_columns(mt_obj):
    new_row      = {}
    frequencies  = mt_obj.Tipper.freq
    tipper_array = mt_obj.Tipper.tipper
    t_n=[]
    t_e=[]
    for index,freq in enumerate(frequencies):
        int_freq = int(freq)
        real_east_key = f"rE_{int_freq:04}"
        imag_east_key = f"iE_{int_freq:04}"
        real_north_key= f"rN_{int_freq:04}"
        imag_north_key= f"iN_{int_freq:04}"
        
        new_row[real_east_key]  = tipper_array[index,0,1].real
        new_row[imag_east_key]  = tipper_array[index,0,1].imag
        new_row[real_north_key] = tipper_array[index,0,0].real
        new_row[imag_north_key] = tipper_array[index,0,0].imag
    return new_row
    
def extract_edi_to_dataframe_row(row):
    mt_obj  = MT(row['filepath'])
    new_row = convert_tipper_to_csv_columns(mt_obj)
    new_row['Lat']    = mt_obj.lat
    new_row['Long']   = mt_obj.lon
    new_row['ALT']    = mt_obj.elev
    new_row['Station']= int(float(mt_obj.station))
    new_row['mt_obj'] = mt_obj
    return new_row
    
    
def convert_edi_files_into_dataframe():
    files = glob.glob('edi_files/*.edi')
    df    = pd.DataFrame(files,columns=['filepath'])
    new_rows = []
    for index,row in df.iterrows():
        new_rows.append(extract_edi_to_dataframe_row(row))
    df    = pd.DataFrame(new_rows)
    df  = df.sort_values(by=['Station'],ascending = [True])
    return df

def check_both_dataframes_are_equivalent(source_df,target_df):
    columns = ['Station', 'Lat', 'Long', 'ALT','rE_0011', 'rE_0022', 'rE_0045', 'rE_0090',
       'rE_0180', 'rE_0360', 'rE_0720', 'rE_1440', 'rN_0011', 'rN_0022',
       'rN_0045', 'rN_0090', 'rN_0180', 'rN_0360', 'rN_0720', 'rN_1440',
       'iE_0011', 'iE_0022', 'iE_0045', 'iE_0090', 'iE_0180', 'iE_0360',
       'iE_0720', 'iE_1440', 'iN_0011', 'iN_0022', 'iN_0045', 'iN_0090',
       'iN_0180', 'iN_0360', 'iN_0720', 'iN_1440']
    for col in columns:
        target_series = target_df[col]
        source_series = source_df[col]
        try:
            np.testing.assert_array_equal(target_series.values,
                                          source_series.values,'{} columns not equal'.format(col))
        except:
            print('******************************')
            print(col)
            print('******************************')
            for s,t in zip(target_series.values,source_series.values):
                print('source:{} target:{} diff:{}'.format(s,t,s-t))
    
check_one_edi_per_station()
df = convert_edi_files_into_dataframe()
check_both_dataframes_are_equivalent(df,new_data)<jupyter_output>******************************
Lat
******************************
source:61.33295316 target:61.332953 diff:1.5999999902760464e-07
source:61.33458935 target:61.334589 diff:3.500000005374204e-07
source:61.33623517 target:61.336235 diff:1.6999999985500835e-07
source:61.33802624 target:61.338026 diff:2.3999999854140697e-07
source:61.33945685 target:61.339457 diff:-1.500000053056283e-07
source:61.34292716 target:61.342927 diff:1.5999999902760464e-07
source:61.34618612 target:61.346186 diff:1.199999957179898e-07
source:61.34794011 target:61.34794 diff:1.1000000199601345e-07
source:61.34967013 target:61.34967 diff:1.2999999654539351e-07
source:61.35133223 target:61.351332 diff:2.2999999771400326e-07
source:61.35316032 target:61.35316 diff:3.199999980552093e-07
source:61.35469588 target:61.354696 diff:-1.199999957179898e-07
source:61.35643206 target:61.356432 diff:6.000000496442226e-08
source:61.35801488 target:61.358015 diff:-1.2000000282341716e-07
source:61.35982270000001 target:61.359823 di[...]<jupyter_text>The utility of doing these kinds of checks is immediately demonstrated. Every column value test failed except for station. However, the tipper data appears to be identical within floating point precision ($\approx$1e-16). Lat-Lon data appears to have dropped off a decimal place ($\approx$1e-7), perhaps as part of MtPy's internal handling of decimal degrees. Still, the data are close enough to be trusted.**WAIT**. Remember when I said we needed to fill those nan values? lets check really quick to see if that happened.<jupyter_code>new_data  = new_data.reset_index()
plot_data(df,direction='N',part='real')
plot_data(df,direction='N',part='imag')
plot_data(df,direction='E',part='real')
plot_data(df,direction='E',part='imag')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11715711144934889, 'fraction_numerical': 0.0504397091155082, 'mean_word_length': 4.374460349920473, 'pattern_counts': {'":': 0, '<': 51, '<?xml version=': 0, '>': 50, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 3}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '209791', 'license': 'MIT', 'n_tokens_mistral': 8325, 'n_tokens_neox': 7440, 'n_words': 2289, 'path': 'Part_1.Data_Processing.ipynb', 'repo_name': 'alexandarj/ELF2020', 'size': '23652'}