data = pd.read_csv('compositos.txt', sep='\t')
litologias = data['lito'].unique()
combinaciones = list(itertools.combinations(litologias, 2))
for lito_1, lito_2 in combinaciones:
litologia_1 = data[data['lito'] == lito_1]
litologia_2 = data[data['lito'] == lito_2]
if litologia_1.empty or litologia_2.empty:
continue
tree_lito_2 = cKDTree(litologia_2[['x', 'y', 'z']])
distancias_minimas_1, _ = tree_lito_2.query(litologia_1[['x', 'y', 'z']])
leyes_correspondientes_1 = litologia_1['cu'].values
resultados_1 = pd.DataFrame({
'distancia': distancias_minimas_1,
'cu': leyes_correspondientes_1
})
resultados_1['intervalo'] = pd.cut(resultados_1['distancia'], bins=np.arange(0, 110, 10), right=False)
ley_media_por_intervalo_1 = resultados_1.groupby('intervalo')['cu'].mean()
tree_lito_1 = cKDTree(litologia_1[['x', 'y', 'z']])
distancias_minimas_2, _ = tree_lito_1.query(litologia_2[['x', 'y', 'z']])
leyes_correspondientes_2 = litologia_2['cu'].values
resultados_2 = pd.DataFrame({
'distancia': distancias_minimas_2,
'cu': leyes_correspondientes_2
})
resultados_2['intervalo'] = pd.cut(resultados_2['distancia'], bins=np.arange(0, 110, 10), right=False)
ley_media_por_intervalo_2 = resultados_2.groupby('intervalo')['cu'].mean()
if ley_media_por_intervalo_1.dropna().empty or ley_media_por_intervalo_2.dropna().empty:
print(f"Combinación Lito {lito_1} vs Lito {lito_2} no tiene suficientes datos en ambas direcciones.")
continue
plt.figure(figsize=(6, 4))
plt.plot(ley_media_por_intervalo_1.index.categories.mid, ley_media_por_intervalo_1.values, marker='o', label=f'Lito {lito_1} respecto a Lito {lito_2}')
plt.plot(-ley_media_por_intervalo_2.index.categories.mid, ley_media_por_intervalo_2.values, marker='x', label=f'Lito {lito_2} respecto a Lito {lito_1}')
plt.xlabel('Intervalos de Distancia (m)')
plt.ylabel('Ley Media de Cu (%)')
plt.title(f'Distancia al contacto (Lito {lito_1} vs Lito {lito_2})')
plt.axvline(0, color='black', linewidth=0.5)
plt.ylim(0, max(ley_media_por_intervalo_1.max(), ley_media_por_intervalo_2.max()) * 1.1)
plt.grid(False)
plt.legend()
plt.show()
