<jupyter_start><jupyter_text>Example of loading MT Stations<jupyter_code># For updating bezpy without needing to restart the entire kernel
%load_ext autoreload
%autoreload 2

import glob
# NOTE: Must load Basemap before geopandas...
from mpl_toolkits.basemap import Basemap

import numpy as np

%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib as mpl

import bezpy<jupyter_output><empty_output><jupyter_text>Set up the area of interest<jupyter_code>frd_lat = 38.205
frd_lon = -77.373
lat_bounds = [frd_lat-3, frd_lat+3]
lon_bounds = [frd_lon-3, frd_lon+3]<jupyter_output><empty_output><jupyter_text>Set up the paths to all the data<jupyter_code>MT_data_folder = "/Users/glucas/work/data/MT_sites/"
# Get all sub folders downloaded from the USArray website
xml_files = sorted(glob.glob(MT_data_folder + "MT_TF_*/*.xml"))
MT_sites = {}
MT_xy_locs = []
ratings = []

for fname in xml_files:
    
    # Reads in an IRIS XML file and returns a Site object
    site = bezpy.mt.read_xml(fname)
    # Only keep within our bounds
    if site.longitude < lon_bounds[0] or site.longitude > lon_bounds[1] or \
       site.latitude < lat_bounds[0] or site.latitude > lat_bounds[1]:
        continue
        
    MT_sites[site.name] = site
    MT_xy_locs.append((site.longitude, site.latitude))
    ratings.append(site.rating)
    
MT_xy_locs = np.array(MT_xy_locs)
print("There are {0} MT sites".format(len(MT_sites)))

fig, ax = plt.subplots(figsize=(12,10))

m = Basemap(ax=ax, llcrnrlon=lon_bounds[0],
            llcrnrlat=lat_bounds[0],
            urcrnrlon=lon_bounds[1],
            urcrnrlat=lat_bounds[1],
            resolution='i', projection='merc',
            lon_0=(lon_bounds[0]+lon_bounds[1])/2+lon_bounds[0],
            lat_0=(lat_bounds[0]+lat_bounds[1])/2+lat_bounds[0])


#m = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
#            projection='lcc',lat_1=33,lat_2=45,lon_0=-95)

# draw coastlines, state and country boundaries, edge of map.
m.drawcoastlines()
m.drawstates()
m.drawcountries()

cmap = plt.cm.Reds

m.scatter(MT_xy_locs[:,0], MT_xy_locs[:,1], latlon=True, s=40, c=cmap(np.array(ratings)/5.))
m.scatter(frd_lon, frd_lat, latlon=True, s=50, marker='x')


bounds = np.arange(7)
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
# fake up the array of the scalar mappable.
sm._A = []
cbar = fig.colorbar(sm, spacing='proportional', ticks=bounds-0.5)
cbar.set_ticklabels(bounds)
cbar.set_label("Station Quality")<jupyter_output>There are 73 MT sites<jupyter_text>Inspecting a Site object<jupyter_code>site = MT_sites["VAP57"]
site
fig = site.plot_apparent_resistivity()<jupyter_output><empty_output><jupyter_text>Now look at what an interpolation through the data looks like<jupyter_code>freqs = 1./np.logspace(np.log10(5), np.log10(1e5), 1e4)
fig = site.plot_apparent_resistivity(interp_freqs=freqs)<jupyter_output><empty_output><jupyter_text>Convolution with an artifical magnetic field<jupyter_code>ts = np.arange(0, 10001)
# 200 second period
magX = np.sin(2*np.pi*ts/200)
magY = np.zeros(len(magX))
plt.plot(ts, magX)
Ex, Ey = site.convolve_fft(magX, magY, dt=1)
plt.plot(ts, magX, label=r'B$_x$', c='C0')
plt.plot(ts, Ex, label=r'E$_x$', c='C1')
plt.plot(ts, Ey, label=r'E$_y$', c='C2')
plt.legend()
# Now flip the magX/magY
magX, magY = magY, magX
Ex, Ey = site.convolve_fft(magX, magY, dt=1)
plt.plot(ts, magY, label=r'B$_y$', c='C0')
plt.plot(ts, Ex, label=r'E$_x$', c='C1')
plt.plot(ts, Ey, label=r'E$_y$', c='C2')
plt.legend()<jupyter_output><empty_output><jupyter_text>Impulse Response<jupyter_code>dtir = bezpy.mt.DTIR(nmin=-600, nmax=3600, dt=1,
                     decay_factor=1.e-12, model_space=True, verbose=True)
dtir
ns, zn = dtir.calculate(site=site)
zf = dtir.get_zf(periods=1./freqs)
# Data
plt.errorbar(site.periods, site.Z[1,:].real, site.Z_var[1,:], c='C0', fmt='o', marker='o')
plt.errorbar(site.periods, site.Z[2,:].real, site.Z_var[2,:], c='C1', fmt='o', marker='x')
# DTIR
plt.semilogx(1./freqs, zf[1,:].real, c='C0')
plt.semilogx(1./freqs, zf[2,:].real, c='C1')
Ex, Ey = dtir.convolve(magX, magY)
plt.plot(ts, magX, label=r'B$_x$', c='C0')
plt.plot(ts, Ex, label=r'E$_x$', c='C1')
plt.plot(ts, Ey, label=r'E$_y$', c='C2')
plt.legend()
magX, magY = magY, magX
Ex, Ey = dtir.convolve(magX, magY)
plt.plot(ts, magY, label=r'B$_y$', c='C0')
plt.plot(ts, Ex, label=r'E$_x$', c='C1')
plt.plot(ts, Ey, label=r'E$_y$', c='C2')
plt.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16637050199911152, 'fraction_numerical': 0.029098178587294535, 'mean_word_length': 5.297902097902098, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 30, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 4}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '154780', 'license': 'MIT', 'n_tokens_mistral': 2017, 'n_tokens_neox': 1928, 'n_words': 352, 'path': 'notebooks/MT_Example.ipynb', 'repo_name': 'erigler-usgs/bezpy', 'size': '4502'}