import numpy as np
import pandas as pd
import dill
from emrecharge.rockphysics import (
from_sigma_to_fraction,
compute_rho_to_cf_mappings_with_rho,
)
import matplotlib.pyplot as plt
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
results_metrics = dill.load(open("results.pik", "rb"))
df_rock_physics = pd.read_csv("rock_physics.csv")
thresholds = [10, 20, 30, 40, 50]
rho_fine_above = df_rock_physics['rho_fine_dominated_above'].median()
rho_coarse_above = df_rock_physics['rho_coarse_dominated_above'].median()
rho_min, rho_max = 13, 33
nx = results_metrics['nx']
ny = results_metrics['ny']
mapping = compute_rho_to_cf_mappings_with_rho(rho_fine_above, rho_coarse_above, rho_min, rho_max)
import matplotlib
matplotlib.rcParams['font.size'] = 12
matplotlib.rcParams['pcolor.shading'] = 'auto'
state = AppState()
METRIC_PATH_LENGTH = 'Path length'
METRIC_NORM_PATH_LENGTH = 'Normalized path length'
METRIC_DEPTH_TO_NO_FLOW = 'Depth to no-flow'
optins_metrics = [
METRIC_PATH_LENGTH,
METRIC_NORM_PATH_LENGTH,
METRIC_DEPTH_TO_NO_FLOW
]
w_metric=widgets.RadioButtons(
options=optins_metrics,
description="Select Recharge Metric",
layout={"width":"500px"},
style = {'description_width': 'initial'},
)
state.register_stateful_widget(w_metric, "metric", Unicode(METRIC_PATH_LENGTH))
THRESHOLD_MIN = 10
THRESHOLD_MAX = 50
THRESHOLD_STEP = 10
THRESHOLD_DEFAULT = 20
THRESHOLDS = list(range(THRESHOLD_MIN, THRESHOLD_MAX+1, THRESHOLD_STEP))
w_threshold=widgets.IntSlider(
min=THRESHOLD_MIN, max=THRESHOLD_MAX, step=THRESHOLD_STEP,
description="FCD Threshold",
value=THRESHOLD_DEFAULT,
disabled=False,
style = {'description_width': 'initial'},
layout={"width":"500px"}
)
state.register_stateful_widget(w_threshold, "threshold", Int(THRESHOLD_DEFAULT))
max_path_length = 500 / 0.3048 * 4
max_water_depth = results_metrics['depth_to_wt'].max()
cmap_values = {}
cmap_values[METRIC_PATH_LENGTH] = (0, max_path_length)
cmap_values[METRIC_NORM_PATH_LENGTH] = (1, 5.)
cmap_values[METRIC_DEPTH_TO_NO_FLOW] = (0, max_water_depth)
cmap_limits = {}
cmap_limits[METRIC_PATH_LENGTH] = (0, 5e3)
cmap_limits[METRIC_NORM_PATH_LENGTH] = (1, 10)
cmap_limits[METRIC_DEPTH_TO_NO_FLOW] = (0, max_water_depth)
w_cmap_range=widgets.FloatRangeSlider(
min=cmap_limits[METRIC_PATH_LENGTH][0],
max=cmap_limits[METRIC_PATH_LENGTH][1],
step=0.1,
value=cmap_values[METRIC_PATH_LENGTH],
continuous_update=False,
description="", readout_format=".1f")
w_cmap_label=widgets.Label(value="Colormap Range:")
box_cmap_controls=widgets.VBox([w_cmap_label, w_cmap_range])
state.register_stateful_widget(w_cmap_range, "cmap_range", List(cmap_values[METRIC_PATH_LENGTH]))
def update_cmap_range(evt):
state.widgets["cmap_range"].min = cmap_limits[evt.owner.metric][0]
state.widgets["cmap_range"].max = cmap_limits[evt.owner.metric][1]
state.widgets["cmap_range"].value = cmap_values[evt.owner.metric]
state.observe(update_cmap_range, ["metric", "threshold"])
def plot_mapping(threshold):
ii = THRESHOLDS.index(threshold)
fig, mapping_ax = plt.subplots(1,1, figsize=(10, 4))
mapping_ax.plot(np.log10(mapping["rho"]), mapping["func_cf_above"](np.log10(mapping["rho"])) * 100, 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("Fraction coarse dominated (%)")
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')
ylim = mapping_ax.get_ylim()
mapping_ax.fill_between((np.log10(rho_min), np.log10(rho_max)), threshold, 100, color='green', alpha=0.2, label='flow')
mapping_ax.fill_between((np.log10(rho_min), np.log10(rho_max)), 0, threshold, color='gold', alpha=0.2, label='no-flow')
mapping_ax.set_ylim(ylim)
mapping_ax.set_xlim(np.log10(rho_min), np.log10(rho_max))
legend = mapping_ax.legend(fontsize=10)
plt.tight_layout()
def plot(**kwargs):
threshold = kwargs['threshold']
metric = kwargs['metric']
vmin, vmax = kwargs['cmap_range']
ii = THRESHOLDS.index(threshold)
if metric == METRIC_PATH_LENGTH:
cmap = 'Greens_r'
values = results_metrics['ds'][0,ii,:].copy()
elif metric == METRIC_NORM_PATH_LENGTH:
cmap = 'Greens_r'
values = results_metrics['ds'][0,ii,:].copy() / results_metrics['depth_to_wt']
elif metric == METRIC_DEPTH_TO_NO_FLOW:
cmap = 'Greens'
values = results_metrics['d_to_no_flows'][0,ii,:].copy()
values[~results_metrics['mask_inds']] = np.nan
fig, ax = plt.subplots(figsize=(10, 8))
out = ax.pcolormesh(
results_metrics['x_node'], results_metrics['y_node'], values.reshape((nx, ny), order='F').T, cmap=cmap, vmin=vmin, vmax=vmax
)
plt.colorbar(out, ax=ax, fraction=0.02)
ax.set_aspect(1)
ax.set_facecolor('lightgrey')
ax.set_xlabel("Easting (ft)")
ax.set_ylabel("Northing (ft)")
title = "{:s}, FCD threshold={:.0f}%".format(metric, threshold)
ax.set_title(title)
w_output = interactive_output(
plot,
{ k: state.widgets[k] for k in state.widgets.keys()}
);
w_output_map = interactive_output(
plot_mapping,
{'threshold':w_threshold}
);
panel = VBox([w_output_map, w_threshold, w_metric, w_output, box_cmap_controls])
display(panel)