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)