supporting notebook

Resistivity to FCD Mapping

Resistivity to FCD mapping App

App parameters

  • ρFD\rho_{FD}: resistivity of the fine-dominated
  • ρCD\rho_{CD}: resistivity of the coarse-dominated
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import dill
from emrecharge.colocation import find_closest_locations, find_locations_in_distance, compute_fraction_for_aem_layer
from emrecharge.rockphysics import (
    compute_interval,
    compute_integral,
    from_sigma_to_fraction,
    compute_rho_to_cf_mappings_with_rho
)
from scipy.stats import percentileofscore
from ipywidgets import Output, BoundedIntText, BoundedFloatText, Label, HBox, VBox, widgets, interactive_output
from IPython.display import display, clear_output
from traitlets import Int, Float
from traitlets import HasTraits, Int, Unicode, List, Float, Bool, TraitType, Dict
from emrecharge.app import AppState
import matplotlib as mpl
df_lith_collar = pd.read_csv("collar.csv")
df_lith_interval = pd.read_csv("interval.csv")
# Code 0: fine-dominated; Code 1: coarse-dominated
df_lith_interval['Code'] = df_lith_interval['Code']-1.
xy_lithology = df_lith_collar[['X', 'Y']].values
df_lith_collar = df_lith_collar.set_index('Hole_ID')
df_lith_interval = df_lith_interval.set_index('Hole_ID')
lith_index = np.array(df_lith_collar.index.to_list())
results = dill.load(open("rho_1d_inversion.pik", "rb"))
xy_em = results['topography'][:,:2]
n_sounding = xy_em.shape[0]
hz = results['hz']
resistivity = results['rho'].reshape((n_sounding, len(hz)))
n_layer = len(hz)

from matplotlib.colors import LogNorm
import matplotlib
matplotlib.rcParams['font.size'] = 12
matplotlib.rcParams['pcolor.shading'] = 'auto'
from matplotlib.ticker import FormatStrFormatter

state = AppState()
style = {'description_width': 'initial'}
w_rho_fd=widgets.FloatSlider(
    min=10, max=25., 
    value=16,
    step=0.1,
    continuous_update=False, 
    disabled=False,
    description='rho_FD',
    style=style,
    layout={"width":"500px", "left": "20px"}
)

w_rho_cd=widgets.FloatSlider(
    min=10, max=25., 
    value=20,
    step=0.1,
    continuous_update=False, 
    disabled=False,
    description='rho_CD',
    style=style,
    layout={"width":"500px", "left": "20px"}
)
state.register_stateful_widget(w_rho_cd, "rho_cd", Float(w_rho_cd.value))
state.register_stateful_widget(w_rho_fd, "rho_fd", Float(w_rho_fd.value))

def plot(**kwargs):
    rho_coarse_above = kwargs['rho_cd']
    rho_fine_above = kwargs['rho_fd']
    topo = results['topography']
    resistivity = results['rho'].reshape((n_sounding, len(hz))).T
    m_to_ft_conversion = 1./0.3048
    depth = np.cumsum(np.r_[0., hz]) * m_to_ft_conversion
    x = topo[:,0]*m_to_ft_conversion
    x_node = np.r_[x[0]-25.*m_to_ft_conversion, x+25.*m_to_ft_conversion]
    rho_min, rho_max = results['rho'].min(), results['rho'].max()

    mapping = compute_rho_to_cf_mappings_with_rho(rho_fine_above, rho_coarse_above, rho_min, rho_max)
    fcd_em = mapping['func_cf_above'](np.log10(results['rho']))
    fcd = fcd_em.reshape((n_sounding, len(hz))).T * 100

    fig, axs = plt.subplots(2,1, figsize=(8, 5))
    ax1, ax2 = axs
    rho_min, rho_max = 10, 40
    # _,_ = model.plot_section(i_line=0, x_axis='x', aspect=5, cmap='turbo', clim=(8, 80), dx=25, physical_property=rho_true, ax=ax2, show_colorbar=False)
    out_rho = ax1.pcolormesh(x_node, depth, resistivity, shading='auto', norm=LogNorm(vmin=rho_min, vmax=rho_max), cmap='turbo')
    out_fcd = ax2.pcolormesh(x_node, depth, fcd, shading='auto', vmin=0, vmax=100, cmap='RdBu_r')
    # _, _= model.plot_section(i_line=0, x_axis='x', aspect=5, cmap='turbo', clim=(8, 80), dx=25, alpha=1, ax=ax2, show_colorbar=False)
    titles = ["(a) Resistivity", "(b) Fraction coarse dominated"]
    for ii, ax in enumerate(axs):
        ax.set_xlim(-2500*m_to_ft_conversion, 2500*m_to_ft_conversion)
        ax.set_ylim(70*m_to_ft_conversion,  0)
        ax.set_aspect(20)
        if ii == 0:
            ax.set_xlabel("")
        else:
            ax.set_xlabel("Easting (ft)")
        ax.set_ylabel("Depth (ft)")
        ax.set_title(titles[ii])

    cb_rho=plt.colorbar(out_rho, ax=ax1, orientation='vertical', fraction=0.02) 
    cb_rho.set_label("Resistivity ($\Omega$m)")
    cb_rho.ax.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
    cb_fcd=plt.colorbar(out_fcd, ax=ax2, orientation='vertical', fraction=0.02) 
    cb_fcd.set_label("FCD (%)")
    cb_fcd.ax.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
    plt.tight_layout()

def plot_map(**kwargs):
    rho_coarse_above = kwargs['rho_cd']
    rho_fine_above = kwargs['rho_fd']
    topo = results['topography']
    resistivity = results['rho'].reshape((n_sounding, len(hz))).T
    m_to_ft_conversion = 1./0.3048
    depth = np.cumsum(np.r_[0., hz]) * m_to_ft_conversion
    x = topo[:,0]*m_to_ft_conversion
    x_node = np.r_[x[0]-25.*m_to_ft_conversion, x+25.*m_to_ft_conversion]
    rho_min, rho_max = results['rho'].min(), results['rho'].max()

    mapping = compute_rho_to_cf_mappings_with_rho(rho_fine_above, rho_coarse_above, rho_min, rho_max)
    fcd_em = mapping['func_cf_above'](np.log10(results['rho']))
    fcd = fcd_em.reshape((n_sounding, len(hz))).T * 100
    # print (fcd.mean())
    fig, mapping_ax = plt.subplots(1,1, figsize=(8, 3))
    mapping_ax.plot(np.log10(mapping["rho"]), mapping["func_cf_above"](np.log10(mapping["rho"])), color='k', alpha=0.5, lw=2)
    mapping_ax.set_xlabel("Resistivity ($\Omega$m)")
    mapping_ax.set_title("Resistivity to FCD mapping")
    mapping_ax.set_ylabel("Coarse fraction")

    xticks = mapping_ax.get_xticks()

    mapping_ax.set_xticks(xticks)
    mapping_ax.set_xticklabels(["{:.0f}".format(10**x) for x in xticks])   
    mapping_ax_1 = mapping_ax.twinx()
    mapping_ax_1.hist(np.log10(results['rho']), bins=np.linspace(np.log10(rho_min), np.log10(rho_max), 40), alpha=0.2, color=None, edgecolor='k', label='EM\nresistivity')
    mapping_ax.set_xlim(np.log10(rho_min), np.log10(rho_max))
    legend = mapping_ax_1.legend(fontsize=10)
    plt.tight_layout()


w_output_map = interactive_output(
    plot_map,
    { k: state.widgets[k] for k in state.widgets.keys()}
);

w_output = interactive_output(
    plot,
    { k: state.widgets[k] for k in state.widgets.keys()}
);

panel = VBox([w_output_map, w_rho_fd, w_rho_cd, w_output])
display(panel)