Detectability App¶
This app provides a capability to explore detectability of a target layer using an electromagnetic (EM) experiment. Detectability can be defined as
where is the number of data; is a EM simulation operator calculating EM datum at i-th time channel; and are vertical resistivity profiles. For instance, and could correspond , a vertical resistivity profile with a layer, and , a homogenous background resistivity, respectively. In such a case, the detectability is an anomalous layer response normalized by the background response in percentage.
App parameters¶
EM system:AEMortTEMUse Profile 2: if checked, the app allow user to edit the second vertical resistivity profile; if not, the second resistvity profile is set to be homogeneous background.z: depth to the layerh: thickness of the layerrho_b: background resistivityrho: layer resistivity
Plotting parameters¶
Resistivity range: range of resistivity shown in the left panelDepth range: range of depth shown in the left panel
import numpy as np
import matplotlib as mpl
from matplotlib import pyplot as plt
mpl.rcParams.update({'font.size':12})
from SimPEG import maps
import SimPEG.electromagnetics.time_domain as tdem
from SimPEG import utils
from ipywidgets import interactive, widgets
from IPython.display import display, HTML, JSON, clear_output
from ipywidgets import (
Label, VBox, HBox,
IntSlider, RadioButtons, Checkbox, Output, SelectionSlider, Dropdown, Button,
AppLayout, Layout, interactive_output, FloatRangeSlider, BoundedIntText, BoundedFloatText
)
from SimPEG.electromagnetics.utils import em1d_utils
from discretize import TensorMesh
from discretize.utils import volume_average
import dill
import matplotlib.pyplot as plt
from emrecharge.app import AppState
from scipy.stats import percentileofscore
import json
import dill
from traitlets import HasTraits, Int, Unicode, List, Float, Bool, TraitType, Dict
fname_aem = './aem_system.json'
fname_ttem = './ttem_system.json'
def get_vertical_discretization(n_layer, minimum_dz, geomtric_factor):
"""
Creates a list of vertical discretizations generate from a geometric series.
>>> minimum_dz * geomtric_factor ** np.arange(n_layer)
Parameters
----------
n_layer : int
The number of discretizations
minimum_dz : float
The smallest discretization
geometric_factor : float
The expansion factor of the discretizations
Returns
(n_layer) : numpy.ndarray
The cell widths
"""
hz = minimum_dz * (geomtric_factor) ** np.arange(n_layer)
return hz
class EMSIMULATION:
"""EMSIMULATION"""
def __init__(self):
self.n_layer = 3
self._initialize()
# Widget variables
def _initialize(self):
# AEM
with open(fname_aem, 'r') as infile:
system_aem = json.load(infile)
unit_conversion = 1e-12
area = system_aem['tx_area']
i_start_hm = 10
i_start_lm = 10
i_end_hm = 10 + 27
i_end_lm = 10 + 20 -2
waveform_hm = np.array(system_aem['hm']['waveform'])
waveform_lm = np.array(system_aem['lm']['waveform'])
time_input_currents_HM = waveform_hm[:,0]
input_currents_HM = waveform_hm[:,1]
time_input_currents_LM = waveform_lm[:,0]
input_currents_LM = waveform_lm[:,1]
time_gates = np.array(system_aem['time_gates'])
GateTimeShift=system_aem['lm']['GateTimeShift']
MeaTimeDelay=system_aem['lm']['MeaTimeDelay']
NoGates=int(system_aem['lm']['NoGates'])
t0_lm = waveform_lm[:,0].max()
times_lm = (time_gates[:NoGates,0] + GateTimeShift + MeaTimeDelay)[i_start_lm:i_end_lm]
GateTimeShift=system_aem['hm']['GateTimeShift']
MeaTimeDelay=system_aem['hm']['MeaTimeDelay']
NoGates=int(system_aem['hm']['NoGates'])
t0_hm = waveform_hm[:,0].max()
times_hm = (time_gates[:NoGates,0] + GateTimeShift + MeaTimeDelay)[i_start_hm:i_end_hm]
receiver_location = np.array([13.25, 0., 30.])
receiver_orientation = "z" # "x", "y" or "z"
times = np.logspace(-5, -2, 31) # time channels
source_location = np.array([0., 0., 30.])
source_radius = 10.
current_amplitude = 1.
source_orientation = 'z'
source_list = []
waveform_lm = tdem.sources.PiecewiseLinearWaveform(times=time_input_currents_LM, currents=input_currents_LM)
waveform_hm = tdem.sources.PiecewiseLinearWaveform(times=time_input_currents_HM, currents=input_currents_HM)
# Define receivers at each location.
dbzdt_receiver_lm = tdem.receivers.PointMagneticFluxTimeDerivative(
receiver_location, times_lm, receiver_orientation,
)
dbzdt_receiver_hm = tdem.receivers.PointMagneticFluxTimeDerivative(
receiver_location, times_hm, receiver_orientation,
)
# Must define the transmitter properties and associated receivers
src_0 = tdem.sources.MagDipole(
[dbzdt_receiver_lm],
location=source_location,
waveform=waveform_lm,
orientation=source_orientation,
)
src_1 = tdem.sources.MagDipole(
[dbzdt_receiver_hm],
location=source_location,
waveform=waveform_hm,
orientation=source_orientation,
)
# this needs to be fixed later
source_list.append(src_0)
source_list.append(src_1)
survey = tdem.Survey(source_list)
wire_map = maps.Wires(("sigma", self.n_layer), ("thickness", self.n_layer-1))
self._simulation_aem = tdem.Simulation1DLayered(
survey=survey, thicknessesMap=wire_map.thickness, sigmaMap=wire_map.sigma,
)
# tTEM
with open(fname_ttem, 'r') as infile:
system_ttem = json.load(infile)
area = system_ttem['tx_area']
i_start_hm = int(system_ttem['hm']['RemoveInitialGates'])
i_start_lm = int(system_ttem['lm']['RemoveInitialGates'])
i_end_lm = int(system_ttem['lm']['RemoveGatesFrom'])
waveform_hm = np.array(system_ttem['hm']['waveform'])
waveform_lm = np.array(system_ttem['lm']['waveform'])
time_input_currents_hm = waveform_hm[:,0] - waveform_hm[:,0].max()
input_currents_hm = waveform_hm[:,1]
time_input_currents_lm = waveform_lm[:,0] - waveform_lm[:,0].max()
input_currents_lm = waveform_lm[:,1]
time_gates = np.array(system_ttem['time_gates'])
GateTimeShift=system_ttem['lm']['GateTimeShift']
MeaTimeDelay=system_ttem['lm']['MeaTimeDelay']
NoGates=int(system_ttem['lm']['NoGates'])
t0_lm = waveform_lm[:,0].max()
times_lm = (time_gates[:NoGates,0] + GateTimeShift + MeaTimeDelay)[i_start_lm:i_end_lm]
GateTimeShift=system_ttem['hm']['GateTimeShift']
MeaTimeDelay=system_ttem['hm']['MeaTimeDelay']
NoGates=int(system_ttem['hm']['NoGates'])
t0_hm = waveform_hm[:,0].max()
times_hm = (time_gates[:NoGates,0] + GateTimeShift + MeaTimeDelay)[i_start_hm:]
receiver_location = np.array([9.28, 0., 0.45])
receiver_orientation = "z" # "x", "y" or "z"
source_location = np.array([0., 0., 0.5])
source_orientation = 'z'
source_list = []
waveform_lm = tdem.sources.PiecewiseLinearWaveform(times=time_input_currents_lm, currents=input_currents_lm)
waveform_hm = tdem.sources.PiecewiseLinearWaveform(times=time_input_currents_hm, currents=input_currents_hm)
# Define receivers at each location.
dbzdt_receiver_lm = tdem.receivers.PointMagneticFluxTimeDerivative(
receiver_location, times_lm, receiver_orientation,
)
dbzdt_receiver_hm = tdem.receivers.PointMagneticFluxTimeDerivative(
receiver_location, times_hm, receiver_orientation,
)
# Must define the transmitter properties and associated receivers
src_0 = tdem.sources.MagDipole(
[dbzdt_receiver_lm],
location=source_location,
waveform=waveform_lm,
orientation=source_orientation,
)
src_1 = tdem.sources.MagDipole(
[dbzdt_receiver_hm],
location=source_location,
waveform=waveform_hm,
orientation=source_orientation,
)
# this needs to be fixed later
source_list.append(src_0)
source_list.append(src_1)
survey_ttem = tdem.Survey(source_list)
self._simulation_ttem = tdem.Simulation1DLayered(
survey=survey_ttem, thicknessesMap=wire_map.thickness, sigmaMap=wire_map.sigma,
)
@property
def times_lm(self):
src = self._simulation.survey.source_list[0]
return src.receiver_list[0].times
@property
def times_hm(self):
src = self._simulation.survey.source_list[1]
return src.receiver_list[0].times
def simulate(self, em_system, z, h, rho_background, rho_layer):
thicknesses = np.array([z, h])
rho = np.array([rho_background, rho_layer, rho_background])
sigma = 1./rho
model = np.r_[sigma, thicknesses]
if em_system == 'AEM':
self._simulation = self._simulation_aem
elif em_system == 'tTEM':
self._simulation = self._simulation_ttem
dpred = self._simulation.dpred(model)
return dpred
em_sim = EMSIMULATION()
def get_cmap_range_value(vmin, vmax, vrange):
return [vmin + 0.0*vrange, vmax]
res_min, res_max = 1, 1000
res_range = res_max-res_min
depth_min, depth_max = 0, 300
depth_range = depth_max-depth_min
state = AppState()
w_em_system=widgets.RadioButtons(options=["AEM", "tTEM"], value='AEM', layout={"width":"200px"})
w_connect_profiles=widgets.Checkbox(value=False, description='Use Profile 2')
state.register_stateful_widget(w_connect_profiles, "connect_profiles", Bool(w_connect_profiles.value))
state.register_stateful_widget(w_em_system, "em_system", Unicode(w_em_system.value))
w_rho_range=widgets.IntRangeSlider(
min=res_min, max=res_max,
value=(3, 300),
step=1,
continuous_update=False
)
state.register_stateful_widget(w_rho_range, "rho_range", List(get_cmap_range_value(res_min, res_max, res_range)))
w_depth_range=widgets.IntRangeSlider(
min=depth_min, max=depth_max,
value=(0, 200),
step=1,
continuous_update=False
)
state.register_stateful_widget(w_depth_range, "depth_range", List(get_cmap_range_value(depth_min, depth_max, depth_range)),)
w_z_1=widgets.FloatText(min=1, max=300, step=1, value=20, description='$z$', continuous_update=False, layout={"width":"300px"})
state.register_stateful_widget(w_z_1, "z_1", Float(w_z_1.value))
w_h_1=widgets.FloatText(min=1, max=300, step=1, value=20, description='$h$', continuous_update=False, layout={"width":"300px"})
state.register_stateful_widget(w_h_1, "h_1", Float(w_h_1.value))
w_rho_background_1=widgets.FloatText(min=1, max=1000, step=1, value=20, description='$\\rho_b$', continuous_update=False, layout={"width":"300px"})
state.register_stateful_widget(w_rho_background_1, "rho_background_1", Float(w_rho_background_1.value))
w_rho_layer_1=widgets.FloatText(min=1, max=1000, step=1, value=50, description='$\\rho$', continuous_update=False, layout={"width":"300px"})
state.register_stateful_widget(w_rho_layer_1, "rho_layer_1", Float(w_rho_layer_1.value))
w_z_2=widgets.FloatText(disabled=not w_connect_profiles.value, min=1, max=300, step=1, value=20, description='$z$', continuous_update=False, layout={"width":"300px"})
state.register_stateful_widget(w_z_2, "z_2", Float(w_z_2.value))
w_h_2=widgets.FloatText(disabled=not w_connect_profiles.value, min=1, max=300, step=1, value=20, description='$h$', continuous_update=False, layout={"width":"300px"})
state.register_stateful_widget(w_h_2, "h_2", Float(w_h_2.value))
w_rho_background_2=widgets.FloatText(disabled=not w_connect_profiles.value, min=1, max=1000, step=1, value=20, description='$\\rho_b$', continuous_update=False, layout={"width":"300px"})
state.register_stateful_widget(w_rho_background_2, "rho_background_2", Float(w_rho_background_2.value))
w_rho_layer_2=widgets.FloatText(disabled=not w_connect_profiles.value, min=1, max=1000, step=1, value=20, description='$\\rho$', continuous_update=False, layout={"width":"300px"})
state.register_stateful_widget(w_rho_layer_2, "rho_layer_2", Float(w_rho_layer_2.value))
def apply_change(change):
if change["new"]:
w_z_2.disabled = False
w_h_2.disabled = False
w_rho_background_2.disabled = False
w_rho_layer_2.disabled = False
else:
w_z_2.disabled = True
w_h_2.disabled = True
w_rho_background_2.disabled = True
w_rho_layer_2.disabled = True
state.widgets["connect_profiles"].observe(apply_change, names=['value'])
def apply_change_system(change):
if change["new"]=='AEM':
w_depth_range.value = (0, 300)
else:
w_depth_range.value = (0, 70)
state.widgets["em_system"].observe(apply_change_system, names=['value'])
def plot(**kwargs):
rho_layer_1 = kwargs['rho_layer_1']
rho_background_1 = kwargs['rho_background_1']
z_1 = kwargs['z_1']
h_1 = kwargs['h_1']
em_system = kwargs['em_system']
z_2 = kwargs['z_2']
h_2 = kwargs['h_2']
rho_range = kwargs['rho_range']
depth_range = kwargs['depth_range']
connect_profiles = kwargs['connect_profiles']
if connect_profiles:
rho_layer_2 = kwargs['rho_layer_2']
rho_background_2 = kwargs['rho_background_2']
else:
rho_layer_2 = rho_background_1
rho_background_2 = rho_background_1
thickness_1 = np.r_[z_1, h_1]
thickness_2 = np.r_[z_2, h_2]
rho_1 = np.r_[rho_background_1, rho_layer_1, rho_background_1]
rho_2 = np.r_[rho_background_2, rho_layer_2, rho_background_2]
dpred_1 = em_sim.simulate(em_system, z_1, h_1, rho_background_1, rho_layer_1)
dpred_2 = em_sim.simulate(em_system, z_2, h_2, rho_background_2, rho_layer_2)
difference = (dpred_1-dpred_2) / dpred_2 * 100
detectability = np.sqrt((difference**2).sum()/len(difference))
fig, axs = plt.subplots(1,2, figsize=(10, 5))
ax1, ax2 = axs
if not connect_profiles:
line_style = '--'
label_1 = 'Layer'
label_2 = 'Background'
else:
line_style = '-.'
label_1 = 'Profile 1'
label_2 = 'Profile 2'
utils.plot_1d_layer_model(thickness_2, rho_2, ax=ax1, **{'color':'C0', 'linestyle':line_style, 'label': label_2})
utils.plot_1d_layer_model(thickness_1, rho_1, ax=ax1, **{'color':'C0', 'label': label_1})
ax2.loglog(em_sim.times_lm * 1e3, -dpred_1[:em_sim.times_lm.size], color='C0', label='LM')
ax2.loglog(em_sim.times_lm * 1e3, -dpred_2[:em_sim.times_lm.size], color='C0', linestyle=line_style)
ax2.loglog(em_sim.times_hm * 1e3, -dpred_1[em_sim.times_lm.size:], color='C1', label='HM')
ax2.loglog(em_sim.times_hm * 1e3, -dpred_2[em_sim.times_lm.size:], color='C1', linestyle=line_style)
ax2.set_ylabel("Voltage (V/A-m$^4$)")
ax2.set_title("EM data (Detectability={:.1f}%)".format(detectability))
ax2.set_xlabel("Time (ms)")
ax1.set_xlabel("Resistivity (Ohm-m)")
ax2.grid(which='both', alpha=0.2)
ax1.set_title("Vertical resistivity profile")
ax1.set_xlim(rho_range[0], rho_range[1])
ax1.set_ylim(depth_range[1], depth_range[0])
ax1.legend(fontsize=10)
ax2.legend(fontsize=10)
plt.tight_layout()
plt.show()
return
# Some texts
top_label_1 = Label(value="EM system:", layout={"width":"100px"})
top_label_2 = Label(value="Resistivity range", layout={"width":"100px"})
top_label_3 = Label(value="Depth range", layout={"width":"100px"})
edit_plot_limit_label = Label(value="Edit plot access limit:", layout={"width":"200px"})
top_panel_1 = HBox([top_label_1, w_em_system, w_connect_profiles], layout={"align_items":"center"})
top_panel_2 = HBox([top_label_2, w_rho_range], layout={"align_items":"center"})
top_panel_3 = HBox([top_label_3, w_depth_range], layout={"align_items":"center"})
# top_panel_4 = HBox([run_simulation_label, w_run_simulation])
left_label = Label(value="Resistivity profile 1:")
left_panel = VBox(
[left_label, w_z_1, w_h_1, w_rho_background_1, w_rho_layer_1], layout={"align_items":"center"}
)
right_label = Label(value="Resistivity profile 2:")
right_panel = VBox(
[right_label, w_z_2, w_h_2, w_rho_background_2, w_rho_layer_2], layout={"align_items":"center"}
)
w_output = interactive_output(
plot,
{ k: state.widgets[k] for k in state.widgets.keys()}
);
bottom_panel = HBox([left_panel, right_panel], layout={"align_items":"center"})
panel = VBox([top_panel_1, bottom_panel, w_output, edit_plot_limit_label, top_panel_2, top_panel_3])
display(panel)Loading...
# %matplotlib widget
# class APP:
# def __init__(self, state):
# self.state = state
# self.scatter_facecolor=[0., 0., 0., 1.]
# self.scatter_zorder=1
# self.line_plot_color='k'
# self.line_plot_linestyle='-.'
# self.line_plot_linewidth=0.5
# self.line_plot_zorder=2
# self.dpi=100
# plt.ioff()
# self.fig, self.axs = plt.subplots(1,2, figsize=(10*0.8, 5*0.8))
# self.fig.canvas.header_visible = False
# self.fig.canvas.toolbar_visible = True
# self.fig.canvas.footer_visible = True
# self.fig.canvas.resizable = False
# self.fig.canvas.capture_scroll = True
# plt.ion()
# plt.rc('font', size=10) # controls default text sizes
# plt.rc('axes', titlesize=8) # fontsize of the axes title
# plt.rc('axes', labelsize=8) # fontsize of the x and y labels
# plt.rc('xtick', labelsize=8) # fontsize of the tick labels
# plt.rc('ytick', labelsize=8) # fontsize of the tick labels
# self.output = Output();
# with self.output:
# self.update_plot();
# # items = self.fig.canvas.toolbar.toolitems
# # new_tools = [items[0], items[3], items[4]]
# # self.fig.canvas.toolbar.toolitems = new_tools
# def update_plot(self, *args):
# rho_layer_1 = self.state.rho_layer_1
# rho_background_1 = self.state.rho_background_1
# z_1 = self.state.z_1
# h_1 = self.state.h_1
# em_system = self.state.em_system
# z_2 = self.state.z_2
# h_2 = self.state.h_2
# rho_range = self.state.rho_range
# depth_range = self.state.depth_range
# connect_profiles = self.state.connect_profiles
# if connect_profiles:
# rho_layer_2 = self.state.rho_layer_2
# rho_background_2 = self.state.rho_background_2
# else:
# rho_layer_2 = rho_background_1
# rho_background_2 = rho_background_1
# thickness_1 = np.r_[z_1, h_1]
# thickness_2 = np.r_[z_2, h_2]
# rho_1 = np.r_[rho_background_1, rho_layer_1, rho_background_1]
# rho_2 = np.r_[rho_background_2, rho_layer_2, rho_background_2]
# dpred_1 = em_sim.simulate(em_system, z_1, h_1, rho_background_1, rho_layer_1)
# dpred_2 = em_sim.simulate(em_system, z_2, h_2, rho_background_2, rho_layer_2)
# difference = (dpred_1-dpred_2) / dpred_2 * 100
# detectability = np.sqrt((difference**2).sum()/len(difference))
# plt.figure(self.fig.number)
# ax1, ax2 = self.axs
# if not connect_profiles:
# line_style = '--'
# label_1 = 'Layer'
# label_2 = 'Background'
# else:
# line_style = '-.'
# label_1 = 'Profile 1'
# label_2 = 'Profile 2'
# utils.plot_1d_layer_model(thickness_2, rho_2, ax=ax1, **{'color':'C0', 'linestyle':line_style, 'label': label_2})
# utils.plot_1d_layer_model(thickness_1, rho_1, ax=ax1, **{'color':'C0', 'label': label_1})
# ax2.loglog(em_sim.times_lm * 1e3, -dpred_1[:em_sim.times_lm.size], color='C0', label='LM')
# ax2.loglog(em_sim.times_lm * 1e3, -dpred_2[:em_sim.times_lm.size], color='C0', linestyle=line_style)
# ax2.loglog(em_sim.times_hm * 1e3, -dpred_1[em_sim.times_lm.size:], color='C1', label='HM')
# ax2.loglog(em_sim.times_hm * 1e3, -dpred_2[em_sim.times_lm.size:], color='C1', linestyle=line_style)
# ax2.set_ylabel("Voltage (V/A-m$^4$)")
# ax2.set_title("EM data (Detectability={:.1f}%)".format(detectability))
# ax2.set_xlabel("Time (ms)")
# ax1.set_xlabel("Resistivity (Ohm-m)")
# ax2.grid(which='both', alpha=0.2)
# ax1.set_title("Vertical resistivity profile")
# ax1.set_xlim(rho_range[0], rho_range[1])
# ax1.set_ylim(depth_range[1], depth_range[0])
# ax1.legend()
# ax2.legend()
# plt.tight_layout()
# plt.show(self.fig)
# return
# # w_plot_output = interactive_output(
# # plot,
# # { k: state.widgets[k] for k in state.widgets.keys() if k != 'colocation_distance' }
# # );
# # self.panel = VBox([top_panel_1, bottom_panel, w_plot_output, edit_plot_limit_label, top_panel_2, top_panel_3])
# @property
# def widget(self):
# return self.output# app = APP(state)
# panel = VBox([top_panel_1, bottom_panel, app.widget, edit_plot_limit_label, top_panel_2, top_panel_3])
# display(panel)