Algoritmo para la estimación de la turbidez sobre el Río Paraná

Autor/a
Afiliación

Víctor Gustavo Gómez

Fecha de publicación

12 de julio de 2025

Resumen

Este sitio web contiene información sobre la estimación de la turbidez por teledetección en la cuenca media del Río Paraná. La turbidez es uno parámetros de interés dentro proyecto Estimar indicadores de calidad de agua en la cuenca media del río Paraná para el desarrollo de un algoritmo mediante técnicas de teledetección satelital (MSECRE0008604), desarrollado por el Grupo de Investigación Sobre Temas Ambientales y Químicos (GISTAQ) de la Universidad Tecnológica Nacional Facultad Regional Resistencia (UTN-FRRe).

Se utilizarán imágenes del satélite Sentinel-2 con corrección automática, de las cuales se obtiene la reflectancia de superficie del agua. Se buscará la relación entre la reflectancia y la turbidez por métodos de regresión tradionales y machine learning. Una vez obtenido el algoritmo que relacione ambas propiedades, se desarrollaran mapas de distribución espacial.

Palabras clave

GISTAQ, UTN, FRRe, Algoritmo, Turbidez, Machine learning, Teledetección

1 Turbidez

La turbidez se refiere a la opacidad o falta de claridad en un líquido provocada por la presencia de partículas suspendidas. Este fenómeno es un indicador clave en el monitoreo de la calidad del agua y su influencia en diferentes ecosistemas es significativa.

La turbidez es un indicador de la calidad del agua, reflejando la presencia de partículas en suspensión. Su medición es crucial para garantizar la potabilidad del agua y la salud de los ecosistemas acuáticos. Este fenómeno puede ser resultado de diversas causas, como la erosión del suelo, la actividad biológica y la contaminación. La comprensión de la turbidez y su impacto es esencial para la gestión de recursos hídricos y la protección del medio ambiente.

La turbidez viene determinada por la dispersión de la luz causada por la materia suspendida en el agua, se obtiene normalmente mediante un turbidímetro, que proporciona medidas en Nephelometric Turbidity Unit (NTU) y mide la dispersión de un rayo de luz en el agua a 90º de la luz incidente [1].

Muchas propiedades, como la clorofila-a (Chl-a), sólidos suspendidos totales (SST) y la materia orgánica disuelta coloreada (CDOM), se utilizan a menudo como indicadores del estado del agua. Estos constituyentes del agua a su vez son responsables de la turbidez.

Existe una fuerte correlación entre turbidez y sólidos suspendidos totales, por lo que se puede estimar SST a partir de la turbidez. Por lo general, es una relación directa, a mayor concentración de SST mayor turbidez.

Existe una relación inversa entre la Turbidez y la profundidad del disco de Secchi (a valores bajos de secchi mayor turbidez), por lo que también se puede estimar turbidez a partir de mediciones de disco de secchi.

1.1 Métodos tradicionales

Tabla 1: Características principales de algoritmos tradicionales para la estimación de turbidez.
Ecuación Bandas (nm) Métricas Aguas Plataforma Referencia
\(1.559e^{35.533 \cdot B03} \\ 1.879e^{37.745(B03 \cdot B5)/(B04+B12)}\) B03, B04, B05, B12 \(R^{2}\), RMSE, MAE Lago1 Sentinel-2 [2]
\(2677.2 \cdot B04^{1.856}\) B04 \(R^{2}\), RMSE, Bias Interiores variadas2 Landsat-8 [3]
\(969-1.5468 \cdot R_{1200nm}+2.07 \frac{B8A}{B02}\) B02, B8A, 1200nm IOA, SI, RMSE, MAE Río3 Landsat-8 [4]
\(y=-1.1+5.8 \frac{B02}{B04} \\ y=3.896-4.186 \frac{B02}{B03}\) B02, B03, B04 \(R^{2}\), RMSE Río4 Landsat-8 [5]
\(y=37661 \cdot B8A^{2}+1845 \cdot B8A <br> y=531.5- \frac{B04}{0.88}\) B04, B8A \(R^{2}\), RMSE, MAPE Estuario5 Pléiades [6]

Múltiples modelos (lineal, logaritmos, inversa, cuadrática, exponencial, potencial) y plataformas (Sentinel-2, Landsat-5 y Landsat-8) emplean el cociente de bandas B04/B03 [7].

Modelos de estimación a partir de Sentinel-2 y Landsat-8 utilizan regresiones lineales, cuadráticas y logarítmicas empleando B02, B03, B04, B01 (con menos apariciones) y cocientes entre éstas [8].

1.2 Métodos de aprendizaje automático

El aprendizaje automático es un subconjunto de la inteligencia artificial que permite que un sistema aprenda y mejore de forma autónoma, sin necesidad de una programación explícita, a través del análisis de grandes cantidades de datos. El aprendizaje automático permite que los sistemas informáticos se ajusten y mejoren continuamente a medida que acumulan más “experiencias”. Por lo tanto, el rendimiento de estos sistemas puede mejorar si se proporcionan conjuntos de datos más grandes y variados para su procesamiento.

Cuando se entrenan modelos de machine learning, cada conjunto de datos y cada modelo necesitan un conjunto diferente de “hiperparámetros”. Los hiperparámetros son variables de configuración externa que se utilizan para administrar el entrenamiento de modelos de machine learning. Controlan de forma directa la estructura, funciones y rendimiento de los modelos. Los hiperparámetros son los parámetros de un modelo de aprendizaje automático, que no se aprenden durante el entrenamiento, sino que se establecen antes de que comience.

El “ajuste de hiperparámetros” permite modificar el rendimiento del modelo para lograr resultados óptimos. Este proceso es una parte fundamental del machine learning. El ajuste de hiperparámetros puede ser manual o automático. A pesar de que el ajuste manual es lento y tedioso, permite entender mejor cómo afectan al modelo las ponderaciones de los hiperparámetros. El proceso de ajuste de hiperparámetros es iterativo, y debe probar diferentes combinaciones de parámetros y valores.

En el aprendizaje automático es importante utilizar técnicas de “validación cruzada” , de modo que el modelo no se centre únicamente en una única porción de sus datos. La validación cruzada o cross-validation es una técnica utilizada para evaluar los resultados de un análisis estadístico y garantizar que son independientes de la partición entre datos de entrenamiento y prueba. La idea básica de la validación cruzada es dividir los datos en conjuntos de entrenamiento y validación, y luego entrenar el modelo en el conjunto de entrenamiento y evaluar su rendimiento en el conjunto de validación. Este proceso se repite varias veces, con diferentes subconjuntos de los datos utilizados para el entrenamiento y la validación, y se calcula el rendimiento promedio.

En los procesos de machine learning supervisado se utilizan diversos algoritmos y técnicas de cálculo, generalmente calculados mediante el uso de programas como R o Python.

Dependiendo del tipo de datos que se usen para el entrenamiento, será de modelo de aprendizaje automático que se use. A grandes rasgos, existen tres tipos de modelos que se usan en el aprendizaje automático: aprendizaje supervisado , no supervisado y por refuerzo.

Consultando el trabajo de otros investigadores, se observa que utilizan principalmente el aprendizaje automático supervisado. Este tipo aprendizaje supervisado utiliza un conjunto de entrenamiento para enseñar a los modelos a producir el resultado deseado. Este conjunto de datos de entrenamiento incluye entradas y salidas correctas, que permiten al modelo aprender con el tiempo. El algoritmo mide su precisión a través de la función de pérdida, ajustando hasta que el error se haya minimizado lo suficiente.

Yang Zhe y otros, utilizaron como datos de entrada la reflectancia de superficie y datos de salida la turbidez, utilizaron los modelos SVR (support vector regression), random forest (RF) y eXtreme Gradiente Boostring (XGBoost). Los hiperparámetros de cada modelo se determinaron mediante una búsqueda en cuadrícula de validación cruzada en Scikit-Learn de Python [9].

Ma Yue y otros, utilizaron varios modelos de aprendizaje automático, usaron Python 3.7 tanto para la predicción de la turbidez del agua y la optimización de la los hiperparámetros [2].

Zhao y otros probaron 14 modelos de machine learning en un estanque de peces con un dispositivo de construction propia, de los cuales ETR, Bagging, RFR, and ABR son los que presentaron un mejor desempeño en la estimación de la turbidez. Los algoritmos se implementaron utilizando Python 3.6 y bibliotecas de aprendizaje scikit [10].

Tabla 2: Características principales de algoritmos de aprendizaje automático para la estimación de turbidez.
Modelo de machine learning Cuerpo de agua Métricas Plataforma Referencia
SVR, ELM ,BP ,CART ,GBT ,RF ,KNN Lagos RMSE, \(R^{2}\), MAE Sentinel-MSI [2]
eXtreme Gradient Boosting (XGBoost), support vector regression (SVR), random forest (RF) Lago RMSE, \(R^{2}\), MAPE Sentinel-2A/B y Landsat-8/9 [9]
linear regression (LR), ridge regression (RR), least absolute shrinkage and selection operator regression(LASSO), elastic net regression (ENR), k-nearest neighbor regression (KNN), Gaussian process regression (GPR), decision tree regression (DTR), support vector regression (SVR), multilayer perceptron regression (MLP), adaptive boosting regression (ABR), gradient boosting regression (GBR), bootstrap aggregating regression (Bagging), random forest regression (RFR), and extreme tree regression (ETR) Estanque de peces MAE, MRSE, MAPE, \(R^{2}\), RE, Acc Propia [10]

2 Procesamiento de datos

Para el procesamiento de los datos se utilizará la librería pandas de Python.

En el proyecto tenemos dos archivos .csv que contienen los datos:

-base_de_datos_lab.csv → contiene resultados de laboratorio

-base_de_datos_gis.csv → contiene datos espectrales

Importamos la librería pandas para usarla, la nombramos como “pd” para simplificar

import pandas as pd

Leemos los archivos .csv por separado y definimos dos DataFrame

Un DataFrame es basicamente una tabla, donde la información se organiza en filas y columnas. Los datos de la misma columna contienen el mismo tipo de datos, pandas agrega por defecto un “índice” que nos ayuda a identificar una fila en particular.

Con la función pd.read_csv le idicamos a pandas que queremos leer archivos .csv.

df1_lab = pd.read_csv(r"D:\GIT\Turbidez\datos\base_de_datos_lab.csv")
df2_gis = pd.read_csv(r"D:\GIT\Turbidez\datos\base_de_datos_gis_acolite.csv")

df1_lab DataFrame de datos provenientes del laboratorio.

df2_gis DataFrame de datos espectrales provenientes del sensor MSI de Sentinel-2.

Nota: Se debió colocar la “r” delante de la dirección para que lea los archivos.

Video de YouTube ¿Qué es un DataFrame?

Verificamos que los datos se han leído y se crearon correctamente ambos DataFrame por separado, con print

df1_lab.head()

df2_gis.head()
punto banda reflect fecha longitud latitud
0 1 B01 0.041367 2023-05-11 -58.868047 -27.464687
1 1 B02 0.046669 2023-05-11 -58.868047 -27.464687
2 1 B03 0.085630 2023-05-11 -58.868047 -27.464687
3 1 B04 0.128799 2023-05-11 -58.868047 -27.464687
4 1 B05 0.128704 2023-05-11 -58.868047 -27.464687

Nota: La primer columna (donde se ven los valores 0,1,2,3,4, …) es el índice que agrega pandas por defecto al DataFrame, esa columna no forma parte del csv original.

Combinamos ambos DataFrame para tener los datos en una única base de datos.

Lo defimos como df_combinado, para esto utilizamos la función pd.merge para realizar la combinación.

df_combinado = pd.merge(df1_lab, df2_gis, on=['latitud', 'longitud','fecha'], how='inner')

on=[‘latitud’, ‘longitud’,‘fecha’] Especifica las columnas por las cuales se unirán los dos DataFrames. En este caso, por coincidencias exactas en latitud, longitud y fecha.

how=‘inner’ Es el tipo de unión, significa que solo se conservarán las filas que tengan coincidencias en ambos DataFrames en las columnas mencionadas.

Verificamos que la combinacion se haya realizado correctamente

df_combinado.head()
fecha longitud latitud param valor punto banda reflect
0 2023-05-11 -58.868047 -27.464687 ph 6.8 1 B01 0.041367
1 2023-05-11 -58.868047 -27.464687 ph 6.8 1 B02 0.046669
2 2023-05-11 -58.868047 -27.464687 ph 6.8 1 B03 0.085630
3 2023-05-11 -58.868047 -27.464687 ph 6.8 1 B04 0.128799
4 2023-05-11 -58.868047 -27.464687 ph 6.8 1 B05 0.128704

Filtramos la turbidez del DataFrame

Del DataFrame combinado, solo nos interesa las filas que contengan los valores de turbidez, ya que es la propiedad de estudio en este sitio web. Por ello nos quedamos con las filas en donde la columna param sea igual a turb.

Creamos un nuevo DataFrame a partir de df_combinado.

df_turbidez = df_combinado[(df_combinado['param'] == 'turb')]

Verificamos

df_turbidez.head()
fecha longitud latitud param valor punto banda reflect
33 2023-05-11 -58.868047 -27.464687 turb 185.5 1 B01 0.041367
34 2023-05-11 -58.868047 -27.464687 turb 185.5 1 B02 0.046669
35 2023-05-11 -58.868047 -27.464687 turb 185.5 1 B03 0.085630
36 2023-05-11 -58.868047 -27.464687 turb 185.5 1 B04 0.128799
37 2023-05-11 -58.868047 -27.464687 turb 185.5 1 B05 0.128704

Nota: Se han eliminado las filas de ph, cond, sol_sus, etc. Se conserva el índice del DataFrame original.

Eliminamos las columnas que no nos interesan

Creamos un nuevo DataFrame a partir de df_turbidez. Con la función .drop especificamos que columnas queremos eliminar.

df_turbidez_banda = df_turbidez.drop(columns=['longitud','latitud','punto','fecha','param'])

df_turbidez_banda.head()
valor banda reflect
33 185.5 B01 0.041367
34 185.5 B02 0.046669
35 185.5 B03 0.085630
36 185.5 B04 0.128799
37 185.5 B05 0.128704

Nota: Solo nos quedamos con las columnas valor,banda y reflect.

Cambiamos el nombre de la columna “valor” por el de “turbidez” Como los valores de esa columna son de turbidez, directamente cambiamos el nombre de la columna con la función .rename

df_turbidez_banda.rename(columns={'valor': 'turbidez'}, inplace=True)

df_turbidez_banda.head()
turbidez banda reflect
33 185.5 B01 0.041367
34 185.5 B02 0.046669
35 185.5 B03 0.085630
36 185.5 B04 0.128799
37 185.5 B05 0.128704

Creamos la tabla final

Usamos la función .pivot_table para que las columnas sean la turbidez y las distintas bandas de Sentinel-2 (B01, B02 , B03…)

df_final = df_turbidez_banda.pivot_table(
    index='turbidez',
    columns='banda',
    values='reflect',  
    )

¿Qué significa cada término en el argumento de la función pivot_table?

index=‘turbidez’ → Las filas serán los valores únicos de ‘turbidez’

columns=‘banda’ → Las columnas serán los valores únicos de ‘banda’ (B01, B02…)

values=‘reflect’ → El contenido de la tabla será lo que haya en la columna ‘reflect’

Creamos archivo final con los valores de turbidez y reflectancia de cada banda

-Forma 1 Creamos un .csv con los datos que nos interesan con la función .to_csv

 df_final.to_csv('Turbidez_FINAL.csv', index=True)

IMPORTANTE: index=True , debe ser así para que el índice que definimos al pivotar sea una columna visible en el archivo csv.

Recordemos que la columna index es solo visible solo por la librería pandas, Si guardamos con index=False, se omite y no se guarda en el csv.

-Forma 2 Si decidimos poner index=False tenemos que usar una función adicional antes de exportar, debido a que la turbidez está en el índice, no como una columna.

Por lo tanto, luego de hacer el pivot se debe agregar una línea de código.

df_final = df_turbidez_banda.pivot_table(
    index='turbidez',
    columns='banda',
    values='reflect',  
    )

df_final = df_final.reset_index() #Línea de código que se debe agregar

reset_index() convierte el índice definido previamente como “turbidez” en una columna normal y reemplaza el índice por uno numérico estándar (0, 1, 2…). que es el predeterminado por pandas.

Finalmente exportamos como archivo csv.

df_final.to_csv('Turbidez_FINAL2.csv', index=False)

Verificación tabla final

import pandas as pd

df = pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')

df.head()
turbidez B01 B02 B03 B04 B05 B06 B07 B08 B11 B12 B8A
0 2.0 0.030264 0.035694 0.052761 0.038619 0.028929 0.017878 0.018247 0.014811 0.008830 0.006162 0.015677
1 3.0 0.018744 0.027450 0.044980 0.028040 0.018781 0.007693 0.007571 0.007623 0.001773 0.000112 0.003305
2 4.0 0.029141 0.036459 0.054631 0.039793 0.029767 0.022395 0.022966 0.019097 0.009430 0.007976 0.020290
3 4.5 0.022306 0.027934 0.045472 0.028727 0.020835 0.014462 0.017479 0.013988 0.000088 0.000361 0.013232
4 5.0 0.029902 0.030043 0.047235 0.033111 0.024060 0.011116 0.012519 0.007006 0.001982 0.000181 0.009074

3 Prueba de modelo lineal

Regresión lineal con método de mínimos cuadrados

Reemplamos el código del tutorial de sklearn por nuestros datos.

Debemos importar las funciones necesarias para:

1- Leer el archivo .csv creado previamente , que contiene los datos de turbidez y los valores de reflectancia para cada banda.

2- Para usar el método de mínimos cuadrados

Importamos las funciones

import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt

pandas → para leer los datos

train_test_split → para dividir los datos en entrenamiento y validación

LinearRegression → para crear el modelo de regresión lineal

mean_squared_error y r2_score → Para medir el desempeño del modelo (RMSE y R2)

matplotlib.pyplot → para realizar gráficos

Leemos los datos de interés y los dividimos en entrenamiento y validación.

Lectura de datos, lo hacemos con pd.read_csv

df = pd.read_csv(r"D:\GIT\Turbidez\Turbidez_FINAL.csv")

X = df[['B05']]  
y = df['turbidez']  

Definimos:

df, un DataFrame que contiene los datos de turbidez y reflectancia;

X para que sea los valores de reflectancia de una banda en particular;

y para que sea la turbidez.

Dividimos en datos para el entrenamiento y validación

X_train, X_test, y_train, y_test = train_test_split(
  X, y, test_size=0.25, shuffle=True, random_state=42 
)

Conjunto de datos para el entrenamiento X_train y X_test

Conjunto de datos para la validación y_train y y_test

Creamos el modelo de regresión

regressor = LinearRegression().fit(X_train, y_train)

LinearRegression() Crea el modelo de regresión lineal llamado “regressor”. Este modelo, se entrena con la función .fit a partir de los datos de entrenamiento (X_train, y_train)

Evaluamos el modelo generado a partir de las métricas de desempeño.

y_pred = regressor.predict(X_test)
p_rmse = mean_squared_error(y_test, y_pred)
p_r2 = r2_score(y_test, y_pred)

Con la función print visualizamos los valores de las métricas de desempeño

print("RMSE", p_rmse) 
print("R2", p_r2)
RMSE 2006.5290988949998
R2 0.6851071928788597

Para ver la ecuación del modelo de regresión

# Coeficientes (pendientes)
coef = regressor.coef_

# Intercepto (ordenada al origen)
intercept = regressor.intercept_

print(f"La ecuación es: y = {coef[0]:.3f}x + {intercept:.3f}")
La ecuación es: y = 1325.285x + -48.720

Visualizamos los resultados comparando el conjunto de entrenamiento y validación.

Código
fig, ax = plt.subplots(ncols=2, figsize=(10, 5), sharex=True, sharey=True)

#Gráfico de entrenamiento
ax[0].plot(
    X_train,
    regressor.predict(X_train),
    linewidth=3,
    color="#17A77E",
    label="Modelo",
)

ax[0].scatter(X_train, y_train, label="Entrenamiento", color = "#9D50A6", alpha = .6)
ax[0].set(xlabel="Banda 5", ylabel="Turbidez", title="Conjunto de entrenamiento")
ax[0].legend()

#Gráfico de validación
ax[1].plot(X_test, y_pred, linewidth=3, color="#17A77E", label="Modelo")
ax[1].scatter(X_test, y_test, label="Validación", color = "#9D50A6", alpha = .6)
ax[1].set(xlabel="Banda 5", ylabel="Turbidez", title="Conjunto de validación")
ax[1].legend()

#Ecuación de la recta
coef = regressor.coef_[0]
intercept = regressor.intercept_
equation = f"y = {coef:.2f}x + {intercept:.2f}"
# Mostrar la ecuación en ambos subgráficos (opcionalmente, puedes usar solo uno)
for a in ax:
    a.text(0.05, 0.95, equation, transform=a.transAxes,
           fontsize=10, verticalalignment='top',
           bbox=dict(boxstyle="round", facecolor="white", alpha=0.7))

fig.suptitle("Regresión lineal")

plt.show()

4 Pruebas de correlación

Para ser más rigurosos, agregamos más etapas durante el entrenamiento de nuestro modelo lineal.

En esta etapa nos interesa estimar la turbidez por medio de modelos lineales sencillos, para encontrar mejores modelos podemos aplicar transformaciones logarítmicas con el objetivo principal de linealizar relaciones no lineales.

Primeramente se probará un modelo sin aplicar ninguna transformaciones en las variables.

En las siguientes secciones, se aplicará logaritmo en la turbidez, en las bandas espectrales y en en ambas, con el fin de evaluar en rendimiento de los modelos y seleccionar los mejores, según sus métricas de desempeño.

Se utilizará el coeficiente de correlación lineal r, para evaluar la relación lineal entre la turdidez y las distintas bandas (con y sin transformación logarítmica).

Esta medida indica cuanto se relacionan dos variables, puede tomar valores desde -1 a +1:

• +1 correlación lineal perfecta positiva.

• 0 sin correlación.

• -1 correlación lineal perfecta negativa.

4.1 Prueba de correlación C1: (sin aplicar logaritmo)

Importamos pandas para leer los datos

Código
import pandas as pd

Datos= pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')
#print (Datos.head())

Calculamos coeficiente de correlacion lineal “r” entre la turbidez y cada banda con la función .corr de pandas.

Código
bandas = [col for col in Datos.columns if col.startswith('B')]

correlaciones = {}

for banda in bandas:
    r = Datos['turbidez'].corr(Datos[banda])
    correlaciones[banda] = r

#Creamos un Data Frame.
df_correlaciones1 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones1 = df_correlaciones1.sort_values(by='r', ascending=False).reset_index(drop=True)

df_correlaciones1['Banda'] = "turb vs " + df_correlaciones1['Banda'].astype(str)
df_correlaciones1.rename(columns={'Banda': 'Combinación 1'}, inplace=True)

df_correlaciones1.to_csv(r'D:\GIT\Turbidez\Datos creados\Correlaciones\C1_turb_vs_banda.csv', index=False)


#Creo archivo csv para usarlo para entrenar modelos
Datos.to_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv", index=False)

df_correlaciones1
Combinación 1 r
0 turb vs B05 0.884580
1 turb vs B06 0.866118
2 turb vs B08 0.857019
3 turb vs B07 0.856416
4 turb vs B04 0.812887
5 turb vs B8A 0.764714
6 turb vs B03 0.643901
7 turb vs B02 0.577170
8 turb vs B01 0.571294
9 turb vs B12 0.395981
10 turb vs B11 0.372982

Gráfica Turbidez en funcion de las distintas bandas

Código
import pandas as pd
import matplotlib.pyplot as plt

# Cargar los datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")

# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'turbidez']

# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'turbidez']

# Cantidad de bandas
n = len(bandas)

# Calcular las filas y columnas para el grid de subplots
filas = (n + 2) // 3  # hasta 3 columnas por fila
columnas = 3

# Crear figura y ejes
fig, axes = plt.subplots(filas, columnas, figsize=(12, filas * 3))
axes = axes.flatten()

# Graficar cada banda
for i, banda in enumerate(bandas):
    ax = axes[i]
    ax.scatter(df[banda], df['turbidez'], alpha=0.6, s=10)
    ax.set_title(f'Turbidez vs {banda}', fontsize=10)
    ax.set_xlabel(banda)
    ax.set_ylabel('Turbidez')
    ax.grid(True)

# Ocultar los subplots vacíos si hay
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])

# Ajustar espacio
plt.tight_layout()
plt.show()

Las mejores bandas para predecir turbidez (por correlación lineal) son:

B05 > B06 > B08 > B07 > B04

4.2 Prueba de correlación C2: logaritmo en la turbidez

Importamos nunpy para operar con funciones matemáticas

Creamos un nuevo DataFrame para aplicarle el logaritmo a la colomna turbidez.

Código
import numpy as np

Datos_turb_log = pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')
Datos_turb_log['turbidez'] = np.log(Datos_turb_log['turbidez'])

#Cambio el nombre la columna "turbidez" luego de aplicar el logaritmo
Datos_turb_log = Datos_turb_log.rename(columns={'turbidez': 'ln_turbidez'})

#Creo archivo csv para usarlo para entrenar modelos
Datos_turb_log.to_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C2_ln_turb_banda.csv", index=False)
#print(Datos_turb_log.head())

Calculamos r entre el ln_turb y cada banda, con la función .corr de pandas.

Código
bandas = [col for col in Datos_turb_log.columns if col.startswith('B')]

correlaciones = {}

for banda in bandas:
    r = Datos_turb_log['ln_turbidez'].corr(Datos_turb_log[banda])
    correlaciones[banda] = r

#Creamos un Data Frame.
df_correlaciones2 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones2 = df_correlaciones2.sort_values(by='r', ascending=False).reset_index(drop=True)

df_correlaciones2['Banda'] = "ln_turb vs " + df_correlaciones2['Banda'].astype(str)
df_correlaciones2.rename(columns={'Banda': 'Combinación 2'}, inplace=True)

df_correlaciones2.to_csv(r'D:\GIT\Turbidez\Datos creados\Correlaciones\C2_ln_turb_vs_banda.csv', index=False)

df_correlaciones2
Combinación 2 r
0 ln_turb vs B05 0.802144
1 ln_turb vs B04 0.731431
2 ln_turb vs B06 0.687748
3 ln_turb vs B07 0.680786
4 ln_turb vs B08 0.674117
5 ln_turb vs B8A 0.563360
6 ln_turb vs B03 0.483412
7 ln_turb vs B02 0.385967
8 ln_turb vs B01 0.338065
9 ln_turb vs B12 0.194032
10 ln_turb vs B11 0.174042

Gráfica

Código
import pandas as pd
import matplotlib.pyplot as plt

# Cargar los datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C2_ln_turb_banda.csv")

# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'ln_turbidez']

# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'ln_turbidez']

# Cantidad de bandas
n = len(bandas)

# Calcular las filas y columnas para el grid de subplots
filas = (n + 2) // 3  # hasta 3 columnas por fila
columnas = 3

# Crear figura y ejes
fig, axes = plt.subplots(filas, columnas, figsize=(12, filas * 3))
axes = axes.flatten()

# Graficar cada banda
for i, banda in enumerate(bandas):
    ax = axes[i]
    ax.scatter(df[banda], df['ln_turbidez'], alpha=0.6, s=10)
    ax.set_title(f'ln_turbidez vs {banda}', fontsize=10)
    ax.set_xlabel(banda)
    ax.set_ylabel('ln_turbidez')
    ax.grid(True)

# Ocultar los subplots vacíos si hay
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])

# Ajustar espacio
plt.tight_layout()
plt.show()

El orden general se mantiene, pero los valores de r bajan un poco al aplicar el logaritmo a la turbidez. Las mejores bandas siguen siendo: B05 > B04 > B06 > B07 > B08

Aunque r haya bajado ligeramente, esas bandas siguen siendo las más prometedoras predictoras lineales.

4.3 Prueba de correlación C3: logaritmo en las bandas

Anteriormente se estaba trabajando con una base de datos de reflectancias procesadas con Sen2Cor. Actualmente trabajamos con una nueva base de datos con las reflectancias corregidas por ACOLITE.

Acolite y Sen2Cor son dos herramientas de software utilizadas para la corrección atmosférica de imágenes satelitales. Acolite se destaca en ambientes acuáticos eutróficos, mientras que Sen2Cor es un procesador más general para la creación de productos Sentinel-2 de nivel 2A.

IMPORTANTE: Este cambio provocó errores en el código a partir de esta sección en adelante, el error se debe a que hay valores de reflectancia que son cero y negativos en la banda 12.

En las pruebas de correlación C3 y C4 aplicamos logaritmo a las bandas. Por defincion de logaritmo no ponemos tomar valores cero ni negativos. Para salvar este error, se procedió a filtrar filas con esos valores.

Creaamos un nuevo DataFrame para aplicarle el logaritmo a las bandas.

Código
import pandas as pd

# Leer los datos
Datos = pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')

# Guardar cantidad original de filas
filas_originales = Datos.shape[0]

# Filtrar filas: nos quedamos solo con aquellas donde **todos los valores sean mayores a 0**
Datos_filtrado = Datos[(Datos > 0).all(axis=1)]

# Calcular cuántas filas se eliminaron
filas_finales = Datos_filtrado.shape[0]
filas_eliminadas = filas_originales - filas_finales

# Mostrar resultado
print(f"Se eliminaron {filas_eliminadas} filas que contenían ceros o valores negativos.")

# Si querés seguir trabajando con el DataFrame filtrado:
DatosC3_C4 = Datos_filtrado
#DatosC3_C4
Se eliminaron 23 filas que contenían ceros o valores negativos.

Aplicamos logaritmo a las bandas

Código
import numpy as np

DatosC3 = DatosC3_C4.copy()

col = [col for col in DatosC3.columns if col.startswith('B')]

DatosC3[col] = np.log(DatosC3[col])

#Cambiamos en nombre las columnas, agremamos ln_ a cada columna

DatosC3.columns = ['ln_' + col for col in DatosC3.columns]

DatosC3.rename(columns={'ln_turbidez': 'turbidez'}, inplace=True)

#Creo archivo csv para usarlo para entrenar modelos
DatosC3.to_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C3_turb_vs_ln_banda.csv", index=False)

#DatosC3

Calculamos r entre turb y ln de cada banda, con la función .corr de pandas.

Código
bandas_ln = [col for col in DatosC3.columns if col.startswith('ln_B')]

correlaciones = {}

for banda in bandas_ln:
    r = DatosC3['turbidez'].corr(DatosC3[banda])
    correlaciones[banda] = r

#Creamos un Data Frame.
df_correlaciones3 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones3 = df_correlaciones3.sort_values(by='r', ascending=False).reset_index(drop=True)

df_correlaciones3['Banda'] = "turb vs " + df_correlaciones3['Banda'].astype(str)
df_correlaciones3.rename(columns={'Banda': 'Combinación 3'}, inplace=True)


df_correlaciones3.to_csv(r'D:\GIT\Turbidez\Datos creados\Correlaciones\C3_turb_vs_ln_banda.csv', index=False)

df_correlaciones3
Combinación 3 r
0 turb vs ln_B05 0.723516
1 turb vs ln_B06 0.712018
2 turb vs ln_B07 0.703632
3 turb vs ln_B08 0.703631
4 turb vs ln_B04 0.683685
5 turb vs ln_B8A 0.631773
6 turb vs ln_B03 0.572489
7 turb vs ln_B02 0.504570
8 turb vs ln_B01 0.485318
9 turb vs ln_B12 0.308112
10 turb vs ln_B11 0.303731

ln_B05 , ln_B06 , ln_B07 , ln_B08

Estas cuatro variables transformadas logarítmicamente se correlacionan de forma fuerte con la turbidez.

Gráfica

Código
import pandas as pd
import matplotlib.pyplot as plt

# Cargar los datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C3_turb_vs_ln_banda.csv")

# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col.startswith('ln_B')]

# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col.startswith('ln_B')]

# Cantidad de bandas
n = len(bandas)

# Calcular las filas y columnas para el grid de subplots
filas = (n + 2) // 3  # hasta 3 columnas por fila
columnas = 3

# Crear figura y ejes
fig, axes = plt.subplots(filas, columnas, figsize=(12, filas * 3))
axes = axes.flatten()

# Graficar cada banda
for i, banda in enumerate(bandas):
    ax = axes[i]
    ax.scatter(df[banda], df['turbidez'], alpha=0.6, s=10)
    ax.set_title(f'Turbidez vs {banda}', fontsize=10)
    ax.set_xlabel(banda)
    ax.set_ylabel('Turbidez')
    ax.grid(True)

# Ocultar los subplots vacíos si hay
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])

# Ajustar espacio
plt.tight_layout()
plt.show()

4.4 Prueba de correlación C4: logaritmo en turbidez y en bandas

Primeramente verificamos que no hayan valores negativos y cero. Ya que por defincion de logaritno no los podremos utilizar.

Código
import pandas as pd

# Leer los datos
Datos = pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')

# Guardar cantidad original de filas
filas_originales = Datos.shape[0]

# Filtrar filas: nos quedamos solo con aquellas donde **todos los valores sean mayores a 0**
Datos_filtrado = Datos[(Datos > 0).all(axis=1)]

# Calcular cuántas filas se eliminaron
filas_finales = Datos_filtrado.shape[0]
filas_eliminadas = filas_originales - filas_finales

# Mostrar resultado
print(f"Se eliminaron {filas_eliminadas} filas que contenían ceros o valores negativos.")

# Si querés seguir trabajando con el DataFrame filtrado:
DatosC3_C4 = Datos_filtrado
#DatosC3_C4
Se eliminaron 23 filas que contenían ceros o valores negativos.

Creaamos un nuevo DataFrame para aplicarle el logaritmo a todas las columnas

Importamos nunpy para operar con funciones matemáticas

Código
import numpy as np

DatosC4 = np.log(DatosC3_C4)
#Cambiamos en nombre las columnas, agremamos ln_ a cada columna
DatosC4.columns = ['ln_' + col for col in DatosC4.columns]

#Creo archivo csv para usarlo para entrenar modelos
DatosC4.to_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv", index=False)

#DatosC4.head()

Calculamos r entre el ln_turb y ln de cada banda, con la función .corr de pandas.

Código
bandas_ln = [col for col in DatosC4.columns if col.startswith('ln_B')]

correlaciones = {}

for banda in bandas_ln:
    r = DatosC4['ln_turbidez'].corr(DatosC4[banda])
    correlaciones[banda] = r

#Creamos un Data Frame.
df_correlaciones4 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones4 = df_correlaciones4.sort_values(by='r', ascending=False).reset_index(drop=True)

df_correlaciones4['Banda'] = "ln_turb vs " + df_correlaciones4['Banda'].astype(str)
df_correlaciones4.rename(columns={'Banda': 'Combinación 4'}, inplace=True)

df_correlaciones4.to_csv(r'D:\GIT\Turbidez\Datos creados\Correlaciones\C4_ln_turb_vs_ln_banda.csv', index=False)

df_correlaciones4
Combinación 4 r
0 ln_turb vs ln_B05 0.871536
1 ln_turb vs ln_B04 0.808946
2 ln_turb vs ln_B06 0.806073
3 ln_turb vs ln_B07 0.797959
4 ln_turb vs ln_B08 0.786710
5 ln_turb vs ln_B8A 0.695210
6 ln_turb vs ln_B03 0.546315
7 ln_turb vs ln_B02 0.444453
8 ln_turb vs ln_B01 0.378123
9 ln_turb vs ln_B11 0.200915
10 ln_turb vs ln_B12 0.187022

Gráfica

Código
import pandas as pd
import matplotlib.pyplot as plt

# Cargar los datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")

# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'ln_turbidez']

# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'ln_turbidez']

# Cantidad de bandas
n = len(bandas)

# Calcular las filas y columnas para el grid de subplots
filas = (n + 2) // 3  # hasta 3 columnas por fila
columnas = 3

# Crear figura y ejes
fig, axes = plt.subplots(filas, columnas, figsize=(12, filas * 3))
axes = axes.flatten()

# Graficar cada banda
for i, banda in enumerate(bandas):
    ax = axes[i]
    ax.scatter(df[banda], df['ln_turbidez'], alpha=0.6, s=10)
    ax.set_title(f'ln_turbidez vs {banda}', fontsize=10)
    ax.set_xlabel(banda)
    ax.set_ylabel('ln_turbidez')
    ax.grid(True)

# Ocultar los subplots vacíos si hay
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])

# Ajustar espacio
plt.tight_layout()
plt.show()

B05, B04, B06, B07, B08 son los mejores candidatos como predictores en un modelo de regresión. Tienen una correlación fuerte con la turbidez, por encima de 0.78.

Comparación con combinaciones anteriores:

Combinación Variables transformadas Mejor r
1 turbidez vs Banda B05 (0.88)
2 ln(turb) vs Banda B05 (0.80)
3 turbidez vs ln(Banda) B05 (0.72)
4 ln(turb) vs ln(Banda) B05 (0.87)

5 Seleción de variables

Se resume en una tabla las pruebas de correlación realizadas.

Código
import pandas as pd 

C1 = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C1_turb_vs_banda.csv")
C2= pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C2_ln_turb_vs_banda.csv")
C3 = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C3_turb_vs_ln_banda.csv")
C4 = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C4_ln_turb_vs_ln_banda.csv")

C_comb = pd.concat([C1, C2 ,C3, C4], axis=1)

C_comb
Combinación 1 r Combinación 2 r Combinación 3 r Combinación 4 r
0 turb vs B05 0.884580 ln_turb vs B05 0.802144 turb vs ln_B05 0.723516 ln_turb vs ln_B05 0.871536
1 turb vs B06 0.866118 ln_turb vs B04 0.731431 turb vs ln_B06 0.712018 ln_turb vs ln_B04 0.808946
2 turb vs B08 0.857019 ln_turb vs B06 0.687748 turb vs ln_B07 0.703632 ln_turb vs ln_B06 0.806073
3 turb vs B07 0.856416 ln_turb vs B07 0.680786 turb vs ln_B08 0.703631 ln_turb vs ln_B07 0.797959
4 turb vs B04 0.812887 ln_turb vs B08 0.674117 turb vs ln_B04 0.683685 ln_turb vs ln_B08 0.786710
5 turb vs B8A 0.764714 ln_turb vs B8A 0.563360 turb vs ln_B8A 0.631773 ln_turb vs ln_B8A 0.695210
6 turb vs B03 0.643901 ln_turb vs B03 0.483412 turb vs ln_B03 0.572489 ln_turb vs ln_B03 0.546315
7 turb vs B02 0.577170 ln_turb vs B02 0.385967 turb vs ln_B02 0.504570 ln_turb vs ln_B02 0.444453
8 turb vs B01 0.571294 ln_turb vs B01 0.338065 turb vs ln_B01 0.485318 ln_turb vs ln_B01 0.378123
9 turb vs B12 0.395981 ln_turb vs B12 0.194032 turb vs ln_B12 0.308112 ln_turb vs ln_B11 0.200915
10 turb vs B11 0.372982 ln_turb vs B11 0.174042 turb vs ln_B11 0.303731 ln_turb vs ln_B12 0.187022

A partir del análisis con las correlaciones entre turbidez (y su logaritmo) y distintas bandas (y sus logaritmos). Se puede observar que en las combinaciones 1 y 4 las correlaciones son mas altas respecto a C2 y C3. Por que que se seguirá el análisis con estas combinaciones:

• Combinación 1: turb vs bandas

• Combinación 4: ln_turb vs ln_bandas

Se entrenarán modelos con ambas combinaciones por separado y se evalurá su desempeño.

Todos los modelos serán sometidos a bootstrapping con enfoque out-of-bag (OOB):

→ Donde se tomará el 75% de los datos con reemplazo para entrenamiento; → Luego se usa el 25% restante (no muestreados) como conjunto de prueba (test).

Boostrapping

Proponemos avanzar hacia un modelo multivariable, utilizando el AIC [11] como criterio de seleccion de variables, este criterio es muy conveniente para encontrar el modelo de mejor balance entre ajuste y complejidad.

¿Qué hace el AIC?

→ Penaliza modelos con demasiadas variables para evitar sobreajuste.

→ Se busca minimizar el AIC: cuanto más bajo, mejor el modelo.

→ No requiere que las variables tengan alta correlación individual, pero sí que mejoren el modelo conjunto.

AIC=n⋅ln(RMSE²)+2k

5.1 Combinación 1: turb vs bandas

Código
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from utils import run_bootstrap_linear_regression_analysis  # Asegurate de que esta función devuelve lo necesario

# Leer los datos
Datos = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")
y = Datos['turbidez']

# Variables
variables_fijas = ['B05']
variables_a_agregar = ['B06', 'B08', 'B07', 'B04', 'B8A', 'B03', 'B02', 'B01', 'B12', 'B11']

# Resultados
resultados = []

# Configuración de bootstrapping
n_iteraciones_bootstrap = 1000

# Entrenamiento incremental
for i in range(len(variables_a_agregar) + 1):
    variables_usadas = variables_fijas + variables_a_agregar[:i]
    X = Datos[variables_usadas].values

    # División entrenamiento/test
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.25, shuffle=True, random_state=42
    )

    # Bootstrapping: obtener coeficientes promedio y métricas de entrenamiento
    coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = run_bootstrap_linear_regression_analysis(
    X_train, y_train.to_numpy(), n_iteraciones=n_iteraciones_bootstrap)


    # Modelo final con coeficientes promedio
    modelo_final = LinearRegression()
    modelo_final.coef_ = coef_prom
    modelo_final.intercept_ = intercept_prom

    # Predicción sobre testeo
    y_pred = modelo_final.predict(X_test)
    r2 = r2_score(y_test, y_pred)
    rmse = np.sqrt(mean_squared_error(y_test, y_pred))

    # R² ajustado
    n_obs = len(y_test)
    n_vars = X_test.shape[1]
    r2_ajustado = 1 - (1 - r2) * (n_obs - 1) / (n_obs - n_vars - 1)

    # AIC
    residuals = y_test - y_pred
    rss = np.sum(residuals ** 2)
    k = X_test.shape[1] + 1  # +1 por el intercepto
    aic = n_obs * np.log(rss / n_obs) + 2 * k

    # Guardar resultados
    resultados.append({
        "variables": ", ".join(variables_usadas),
        "num_variables": len(variables_usadas),
        "R²_train (bootstrap)": r2_train_boot,
        "RMSE_train (bootstrap)": rmse_train_boot,
        "R²_test": r2,
        "R²-ajustado": r2_ajustado,
        "RMSE_test": rmse,
        "AIC": aic
    })

# Convertir a DataFrame
df_resultados = pd.DataFrame(resultados)
df_resultados
variables num_variables R²_train (bootstrap) RMSE_train (bootstrap) R²_test R²-ajustado RMSE_test AIC
0 B05 1 0.812785 35.553960 0.684411 0.658112 44.843791 110.489184
1 B05, B06 2 0.814920 35.350663 0.680995 0.622995 45.085807 112.639890
2 B05, B06, B08 3 0.821109 34.754564 0.678413 0.581937 45.267941 114.752774
3 B05, B06, B08, B07 4 0.841697 32.693595 0.699094 0.565358 43.788175 115.822186
4 B05, B06, B08, B07, B04 5 0.911450 24.451943 0.885156 0.813378 27.051801 104.337100
5 B05, B06, B08, B07, B04, B8A 6 0.945052 19.261674 0.916569 0.845057 23.057077 101.863237
6 B05, B06, B08, B07, B04, B8A, B03 7 0.944640 19.333827 0.919647 0.825902 22.627765 103.336975
7 B05, B06, B08, B07, B04, B8A, B03, B02 8 0.969925 14.250175 0.949158 0.867811 17.999165 98.929111
8 B05, B06, B08, B07, B04, B8A, B03, B02, B01 9 0.975921 12.750825 0.963266 0.880614 15.299489 96.378945
9 B05, B06, B08, B07, B04, B8A, B03, B02, B01, B12 10 0.981370 11.215662 0.976452 0.897960 12.249413 92.153385
10 B05, B06, B08, B07, B04, B8A, B03, B02, B01, B... 11 0.981549 11.161500 0.975182 0.838680 12.575594 94.889222

Gráfico de AIC

Código
aic_scores = [r["AIC"] for r in resultados]
n_vars = [r["num_variables"] for r in resultados]

plt.plot(n_vars, aic_scores, marker='o', color='purple')
plt.title('AIC vs Número de Variables')
plt.xlabel('Número de variables')
plt.ylabel('AIC')
plt.grid(True)
plt.show()

¿Por qué el AIC baja al agregar más variables, si se supone que penaliza la complejidad?

AIC=n⋅ln(RMSE²)+2k

donde:

n = número de observaciones k = número de parámetros (incluye el intercepto) ln(RMSE²) representa el ajuste del modelo

AIC no busca modelos simples en sí, sino la mejor compensación entre ajuste y complejidad.

A pesar de agregar más variables, el error (RMSE) está disminuyendo mucho, más rápido que lo que penaliza el AIC con +2k

5.2 Combinación 4: ln_turb vs ln_bandas

Código
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from utils import run_bootstrap_linear_regression_analysis  # Asegurate de que esta función devuelve lo necesario

# Leer los datos
Datos = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")
y = Datos['ln_turbidez']

# Variables
variables_fijas = ['ln_B05']
variables_a_agregar = ['ln_B04', 'ln_B06', 'ln_B07', 'ln_B08', 'ln_B8A', 'ln_B03', 'ln_B02', 'ln_B01','ln_B11' ,'ln_B12']

# Resultados
resultados = []

# Configuración de bootstrapping
n_iteraciones_bootstrap = 1000

# Entrenamiento incremental
for i in range(len(variables_a_agregar) + 1):
    variables_usadas = variables_fijas + variables_a_agregar[:i]
    X = Datos[variables_usadas].values

    # División entrenamiento/test
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.25, shuffle=True, random_state=42
    )

    # Bootstrapping: obtener coeficientes promedio y métricas de entrenamiento
    coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = run_bootstrap_linear_regression_analysis(
    X_train, y_train.to_numpy(), n_iteraciones=n_iteraciones_bootstrap)


    # Modelo final con coeficientes promedio
    modelo_final = LinearRegression()
    modelo_final.coef_ = coef_prom
    modelo_final.intercept_ = intercept_prom

    # Predicción sobre testeo
    y_pred = modelo_final.predict(X_test)
    r2 = r2_score(y_test, y_pred)
    rmse = np.sqrt(mean_squared_error(y_test, y_pred))

    # Cálculo R² ajustado
    n_obs = len(y_test)
    n_vars = X_test.shape[1]

    if n_obs - n_vars - 1 != 0:
     r2_ajustado = 1 - (1 - r2) * (n_obs - 1) / (n_obs - n_vars - 1)
    else:
     r2_ajustado = np.nan  # o podrías usar 0.0 si preferís



    # AIC
    residuals = y_test - y_pred
    rss = np.sum(residuals ** 2)
    k = X_test.shape[1] + 1  # +1 por el intercepto
    aic = n_obs * np.log(rss / n_obs) + 2 * k

    # Guardar resultados
    resultados.append({
        "variables": ", ".join(variables_usadas),
        "num_variables": len(variables_usadas),
        "R²_train (bootstrap)": r2_train_boot,
        "RMSE_train (bootstrap)": rmse_train_boot,
        "R²_test": r2,
        "R²-ajustado": r2_ajustado,
        "RMSE_test": rmse,
        "AIC": aic
    })

# Convertir a DataFrame
df_resultados = pd.DataFrame(resultados)
df_resultados
variables num_variables R²_train (bootstrap) RMSE_train (bootstrap) R²_test R²-ajustado RMSE_test AIC
0 ln_B05 1 0.817553 0.616178 0.314326 0.200047 0.900744 2.327460
1 ln_B05, ln_B04 2 0.900090 0.455976 0.646828 0.505559 0.646452 -0.980108
2 ln_B05, ln_B04, ln_B06 3 0.929909 0.381918 0.755352 0.571866 0.538039 -1.917199
3 ln_B05, ln_B04, ln_B06, ln_B07 4 0.931418 0.377784 0.741808 0.397551 0.552732 0.513886
4 ln_B05, ln_B04, ln_B06, ln_B07, ln_B08 5 0.930868 0.379296 0.731719 0.061016 0.563427 2.820532
5 ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A 6 0.963140 0.276959 0.682023 -1.225842 0.613396 6.180090
6 ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A... 7 0.964138 0.273184 0.645954 NaN 0.647251 9.039664
7 ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A... 8 0.965335 0.268586 0.684545 3.208184 0.610958 10.116370
8 ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A... 9 0.967875 0.258560 0.705021 2.032425 0.590797 11.579466
9 ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A... 10 0.967400 0.260462 0.720219 1.652821 0.575376 13.156290
10 ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A... 11 0.959762 0.289372 0.454362 1.954866 0.803517 20.499890

Gráfico de AIC

Código
aic_scores = [r["AIC"] for r in resultados]
n_vars = [r["num_variables"] for r in resultados]

plt.plot(n_vars, aic_scores, marker='o', color='purple')
plt.title('AIC vs Número de Variables')
plt.xlabel('Número de variables')
plt.ylabel('AIC')
plt.grid(True)
plt.show()

6 Entrenamiento de Modelos

6.1 Combinación 1: Turbidez vs Bandas

6.1.1 Modelo con B05

Código
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from utils import run_bootstrap_linear_regression_analysis

# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")

# Variables predictoras
variables = ['B05']
X_completo = df[variables].values
y_completo = df['turbidez'].values

# Separar entrenamiento y testeo
X_train, X_test, y_train, y_test = train_test_split(
    X_completo, y_completo, test_size=0.25, random_state=42
)

# Ejecutar bootstrapping
n_iteraciones_config = 1000
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = \
    run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones=n_iteraciones_config)

# Modelo final con coeficientes promedio
modelo_final_promedio = LinearRegression()
modelo_final_promedio.coef_ = coef_prom
modelo_final_promedio.intercept_ = intercept_prom

# Predicción sobre entrenamiento y test
y_train_pred = modelo_final_promedio.predict(X_train)
y_test_pred = modelo_final_promedio.predict(X_test)

# Métricas en testeo
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))

# ------------------ Gráfico con estilo personalizado ------------------ #
plt.figure(figsize=(5, 5))

# Entrenamiento
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')

# Test
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7, label="Datos de testeo", marker='^')

# Línea ideal
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")

# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")

plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
    f"Regresión lineal con Bootstrapping ({n_iteraciones_config} modelos)\n"
    f"R² entrenamiento: {r2_train_boot:.4f}, RMSE: {rmse_train_boot:.4f} | "
    f"R² testeo: {r2_test:.4f}, RMSE: {rmse_test:.4f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

# ------------------ Coeficientes ------------------ #

# Diccionario para renombrar las variables en notación matemática
nombre_vars_latex = {
    'B05': 'B05',
    'B04': 'B04',
    'B06': 'B06'
}

# Construir términos con nombres renombrados
terminos_latex = [
    f"{coef:.4f} \\cdot {nombre_vars_latex[var]}"
    for var, coef in zip(variables, coef_prom)
]

# Construir la ecuación completa en LaTeX
ecuacion_latexB05 = " + ".join(terminos_latex)
ecuacion_latexB05 = f"turbidez = {ecuacion_latexB05} + {intercept_prom:.4f}"

from IPython.display import display, Math
display(Math(ecuacion_latexB05))

\(\displaystyle turbidez = 1291.6163 \cdot B05 + -46.2433\)

6.1.2 Modelo con 6 bandas

Código
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from utils import run_bootstrap_linear_regression_analysis

# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")

# Variables predictoras
variables = ['B05', 'B06', 'B08', 'B07', 'B04', 'B8A']
X_completo = df[variables].values
y_completo = df['turbidez'].values

# Separar entrenamiento y testeo
X_train, X_test, y_train, y_test = train_test_split(
    X_completo, y_completo, test_size=0.25, random_state=42
)

# Ejecutar bootstrapping
n_iteraciones_config = 1000
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = \
    run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones=n_iteraciones_config)

# Modelo final con coeficientes promedio
modelo_final_promedio = LinearRegression()
modelo_final_promedio.coef_ = coef_prom
modelo_final_promedio.intercept_ = intercept_prom

# Predicción sobre entrenamiento y test
y_train_pred = modelo_final_promedio.predict(X_train)
y_test_pred = modelo_final_promedio.predict(X_test)

# Métricas en testeo
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))

# ------------------ Gráfico con estilo personalizado ------------------ #
plt.figure(figsize=(5, 5))

# Entrenamiento
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')

# Test
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7, label="Datos de testeo", marker='^')

# Línea ideal
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")

# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")

plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
    f"Regresión lineal con Bootstrapping ({n_iteraciones_config} modelos)\n"
    f"R² entrenamiento: {r2_train_boot:.4f}, RMSE: {rmse_train_boot:.4f} | "
    f"R² testeo: {r2_test:.4f}, RMSE: {rmse_test:.4f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

# ------------------ Coeficientes ------------------ #

# Diccionario para renombrar las variables en notación matemática
nombre_vars_latex = {
    'B05': 'B05',
    'B04': 'B04',
    'B06': 'B06',
    'B07': 'B07',
    'B8A': 'B8A',
    'B08': 'B08'
}

# Construir términos con nombres renombrados
terminos_latex = [
    f"{coef:.4f} \\cdot {nombre_vars_latex[var]}"
    for var, coef in zip(variables, coef_prom)
]

# Construir la ecuación completa en LaTeX
ecuacion_latex = " + ".join(terminos_latex)
ecuacion_latex = f"turbidez = {ecuacion_latex} + {intercept_prom:.4f}"

from IPython.display import display, Math
display(Math(ecuacion_latex))

\(\displaystyle turbidez = -261.1831 \cdot B05 + 4596.6799 \cdot B06 + -153.5938 \cdot B08 + 1699.0469 \cdot B07 + -652.0188 \cdot B04 + -4532.9188 \cdot B8A + -5.0671\)

6.2 Modelos con logaritmo en turbidez y bandas

6.2.1 Modelo con 1 variable (ln_B05)

Escala logarítmica

Código
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from utils import run_bootstrap_linear_regression_analysis

# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")

# Variables predictoras
variables = ['ln_B05']
X_completo = df[variables].values
y_completo = df['ln_turbidez'].values

# Separar entrenamiento y testeo
X_train, X_test, y_train, y_test = train_test_split(
    X_completo, y_completo, test_size=0.25, random_state=42
)

# Ejecutar bootstrapping
n_iteraciones_config = 1000
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = \
    run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones=n_iteraciones_config)

# Modelo final con coeficientes promedio
modelo_final_promedio = LinearRegression()
modelo_final_promedio.coef_ = coef_prom
modelo_final_promedio.intercept_ = intercept_prom

# Predicción sobre entrenamiento y test
y_train_pred = modelo_final_promedio.predict(X_train)
y_test_pred = modelo_final_promedio.predict(X_test)

# Métricas en testeo
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))

# ------------------ Gráfico con estilo personalizado ------------------ #
plt.figure(figsize=(5, 5))

# Entrenamiento
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')

# Test
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7, label="Datos de testeo", marker='^')

# Línea ideal
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")

# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")

plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
    f"Regresión lineal con Bootstrapping ({n_iteraciones_config} modelos)\n"
    f"R² entrenamiento: {r2_train_boot:.4f}, RMSE: {rmse_train_boot:.4f} | "
    f"R² testeo: {r2_test:.4f}, RMSE: {rmse_test:.4f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

# ------------------ Coeficientes ------------------ #

# Diccionario para renombrar las variables en notación matemática
nombre_vars_latex = {
    'ln_B05': 'ln(B05)',
    'ln_B04': 'ln(B04)',
    'ln_B06': 'ln(B06)'
}

# Construir términos con nombres renombrados
terminos_latex = [
    f"{coef:.4f} \\cdot {nombre_vars_latex[var]}"
    for var, coef in zip(variables, coef_prom)
]

# Construir la ecuación completa en LaTeX
ecuacion_latex = " + ".join(terminos_latex)
ecuacion_latex = f"\\ln(\\text{{turbidez}}) = {ecuacion_latex} + {intercept_prom:.4f}"

from IPython.display import display, Math
display(Math(ecuacion_latex))

\(\displaystyle \ln(\text{turbidez}) = 1.7728 \cdot ln(B05) + 8.0888\)

Escala real

Código
# Transformo a escala real

y_train_pred_real = np.exp(y_train_pred)
y_test_pred_real = np.exp(y_test_pred)

y_train_real = np.exp(y_train)
y_test_real = np.exp(y_test)

rmse_test_real = np.sqrt(mean_squared_error(y_test_real, y_test_pred_real))
r2_test_real = r2_score(y_test_real, y_test_pred_real)

#Grafico

plt.figure(figsize=(5, 5))

# Entrenamiento (en escala real)
plt.scatter(y_train_real, y_train_pred_real, color="#9D50A6", alpha=0.5, label="Entrenamiento")

# Test (en escala real)
plt.scatter(y_test_real, y_test_pred_real, color="red", alpha=0.7, label="Test")

# Línea ideal
min_val = min(y_train_real.min(), y_test_real.min())
max_val = max(y_train_real.max(), y_test_real.max())
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")

# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")

plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
    f"Regresión lineal en escala real\n"
    f"R² testeo: {r2_test_real:.4f}, RMSE: {rmse_test_real:.4f} NTU"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

6.2.2 Modelo con 3 variables

Escala Logarítmica

Código
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from utils import run_bootstrap_linear_regression_analysis

# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")

# Variables predictoras
variables = ['ln_B05','ln_B04','ln_B06']
X_completo = df[variables].values
y_completo = df['ln_turbidez'].values

# Separar entrenamiento y testeo
X_train, X_test, y_train, y_test = train_test_split(
    X_completo, y_completo, test_size=0.25, random_state=42
)

# Ejecutar bootstrapping
n_iteraciones_config = 1000
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = \
    run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones=n_iteraciones_config)

# Modelo final con coeficientes promedio
modelo_final_promedio = LinearRegression()
modelo_final_promedio.coef_ = coef_prom
modelo_final_promedio.intercept_ = intercept_prom

# Predicción sobre entrenamiento y test
y_train_pred = modelo_final_promedio.predict(X_train)
y_test_pred = modelo_final_promedio.predict(X_test)

# Métricas en testeo
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))

# ------------------ Gráfico con estilo personalizado ------------------ #
plt.figure(figsize=(5, 5))

# Entrenamiento
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')

# Test
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7, label="Datos de testeo", marker='^')

# Línea ideal
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")

# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")

plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
    f"Regresión lineal con Bootstrapping ({n_iteraciones_config} modelos)\n"
    f"R² entrenamiento: {r2_train_boot:.4f}, RMSE: {rmse_train_boot:.4f} | "
    f"R² testeo: {r2_test:.4f}, RMSE: {rmse_test:.4f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

# ------------------ Coeficientes ------------------ #

# Diccionario para renombrar las variables en notación matemática
nombre_vars_latex = {
    'ln_B05': 'ln(B05)',
    'ln_B04': 'ln(B04)',
    'ln_B06': 'ln(B06)'
}

# Construir términos con nombres renombrados
terminos_latex = [
    f"{coef:.4f} \\cdot {nombre_vars_latex[var]}"
    for var, coef in zip(variables, coef_prom)
]

# Construir la ecuación completa en LaTeX
ecuacion_latex = " + ".join(terminos_latex)
ecuacion_latex = f"\\ln(\\text{{turbidez}}) = {ecuacion_latex} + {intercept_prom:.4f}"

from IPython.display import display, Math
display(Math(ecuacion_latex))

# Ecuacion real guardada como string
# Exponenciar el intercepto
intercepto_escala_real = np.exp(intercept_prom)

# Armar términos con potencias
terminos_real = [
    f"{var[3:]}^{{{coef:.4f}}}"  # Elimina "ln_" del nombre de la banda
    for var, coef in zip(variables, coef_prom)
]

# Construir la ecuación en LaTeX como string reutilizable
ecuacion_real_latex = f"\\text{{turbidez}} = {intercepto_escala_real:.4f} \\cdot " + " \\cdot ".join(terminos_real)

# (Opcional) Mostrar
from IPython.display import Math, display
display(Math(ecuacion_real_latex))

# Guardar en otra variable si querés reutilizar luego como string sin formato LaTeX
ecuacion_turb = f"turbidez = {intercepto_escala_real:.4f} * " + " * ".join([
    f"{var[3:]}^{coef:.4f}" for var, coef in zip(variables, coef_prom)
])

\(\displaystyle \ln(\text{turbidez}) = 7.7188 \cdot ln(B05) + -5.3687 \cdot ln(B04) + -1.1799 \cdot ln(B06) + 6.6959\)

\(\displaystyle \text{turbidez} = 809.0508 \cdot B05^{7.7188} \cdot B04^{-5.3687} \cdot B06^{-1.1799}\)

Escala Real

Código
# Transformo a escala real

y_train_pred_real = np.exp(y_train_pred)
y_test_pred_real = np.exp(y_test_pred)

y_train_real = np.exp(y_train)
y_test_real = np.exp(y_test)

rmse_test_real = np.sqrt(mean_squared_error(y_test_real, y_test_pred_real))
r2_test_real = r2_score(y_test_real, y_test_pred_real)

#Grafico

plt.figure(figsize=(5, 5))

# Entrenamiento (en escala real)
plt.scatter(y_train_real, y_train_pred_real, color="#9D50A6", alpha=0.5, label="Entrenamiento")

# Test (en escala real)
plt.scatter(y_test_real, y_test_pred_real, color="red", alpha=0.7, label="Test")

# Línea ideal
min_val = min(y_train_real.min(), y_test_real.min())
max_val = max(y_train_real.max(), y_test_real.max())
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")

# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")

plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
    f"Regresión lineal en escala real\n"
    f"R² testeo: {r2_test_real:.4f}, RMSE: {rmse_test_real:.4f} NTU"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

7 Procesamiento de productos satelitales

Para el procesamiento se utilizará la librería rasterio en python.

Lectura de imágen raster y visualización del área de estudio.

Código
import rasterio as rs
import matplotlib.pyplot as plt
import numpy as np

ruta = r"D:\GIT\Turbidez\recorte_acolite\2025-06-09.tif"

# Abrir el archivo
with rs.open(ruta) as raster:

    #Definimos las bandas 
    #Bandas RGB para ver área de estudio
    B02 = raster.read(2) #Azul
    B03 = raster.read(3) #Verde
    B04 = raster.read(4) #Rojo

    #Para calcular NDWI con B03 y B11
    B11 = raster.read(10)

    #para el modelo
    B05 = raster.read(5)
    B06 = raster.read(6)
    B07 = raster.read(7)

# Composión de imagen 
rgb = np.stack([B04, B03, B02], axis=-1)

# Mejorar constraste
def stretch(banda):
    p2, p98 = np.percentile(banda, (2, 98))
    return np.clip((banda - p2) / (p98 - p2 + 1e-6), 0, 1)

rgb_stretched = np.stack([stretch(B04), stretch(B03), stretch(B02)], axis=-1)

# Gráfico
plt.imshow(rgb_stretched)
plt.title("Río Paraná")
plt.show()

Se utilizrá el índice NDWI para recortar la imágen satelital, y poder trabajar únicamente con los píxeles de agua.

\(NDWI=(B03-B11)/(B03+B11)\)

Código
import numpy as np

ndwi = (B03 - B11) / (B03 + B11 + 1e-6)

#+ 1e-6 para evitar división por cero

# Binarización: agua = 1, no agua = 0
ndwi_bin = np.where(ndwi >= 0.6, 1, 0)

# Gráfico máscara de agua
plt.imshow(ndwi_bin, cmap='Greys')
plt.title("NDWI")
plt.colorbar(label="Agua (1) / No Agua (0)")
plt.show()

7.1 Modelo sin logaritmo con B05

Aplicamos la ecuación hallada para ver la distribución de la turbidez en el río.

Código
import numpy as np
import matplotlib.pyplot as plt
import re

# --- Ecuación como string ---
ecuacion_latexB05 = "turbidez = 1290.3786 * B05 + 10.5"

# --- Extraer coeficiente e intercepto ---
match = re.match(r"turbidez\s*=\s*([-\d.]+)\s*\*\s*B05\s*([\+\-]\s*[\d.]+)", ecuacion_latexB05)
coef_B05 = float(match.group(1))
intercepto = float(match.group(2).replace(" ", ""))  # Quitar espacios en el intercepto

# --- Calcular turbidez directamente ---
turbidez = coef_B05 * B05 + intercepto

# --- Aplicar máscara NDWI (solo donde hay agua) ---
turbidez_agua = np.where(ndwi_bin == 1, turbidez, np.nan)

# --- Calcular percentiles para ajustar visualización ---
p5, p95 = np.nanpercentile(turbidez_agua, [5, 95])

Mapa de turbidez

Código
plt.figure(figsize=(5.5, 5.5))
img = plt.imshow(turbidez_agua, cmap="rainbow", vmin=p5, vmax=p95)
plt.title("Turbidez 2025-06-09")
plt.colorbar(img, label="Turbidez (NTU)")
plt.axis("off")
plt.tight_layout()
plt.show()

7.2 Modelo con logaritmo (B04, B05 ,B06)

Aplicamos las ecuaciones halladas para ver la distribución de la turbidez en el río.

Código
import numpy as np
import matplotlib.pyplot as plt
import rasterio as rs
import re

# --- Ecuación dinámica como string ---
ecuacion_turb = "turbidez = 822.8952 * B05^7.6903 * B04^-5.3331 * B06^-1.1797"

# --- Extraer intercepto y coeficientes ---
# Formato esperado: turbidez = intercepto * Bxx^coef * ...
match = re.match(r"turbidez\s*=\s*([0-9.]+)\s*\*\s*(.+)", ecuacion_turb)
intercepto = float(match.group(1))
terminos = match.group(2)

# Obtener variables y coeficientes
patron = r"(B\d{2})\^(-?\d+\.\d+)"
variables, coeficientes = zip(*re.findall(patron, terminos))
coeficientes = np.array([float(c) for c in coeficientes])
intercept_ln = np.log(intercepto)

# --- Preparar diccionario de bandas ---
bandas = {
    "B04": np.clip(B04, 1e-3, None),
    "B05": np.clip(B05, 1e-3, None),
    "B06": np.clip(B06, 1e-3, None)
}

# --- Calcular ln(turbidez) con log(banda) * coef ---
ln_turbidez = intercept_ln
for var, coef in zip(variables, coeficientes):
    ln_turbidez += coef * np.log(bandas[var])

turbidez = np.exp(ln_turbidez)

# --- Aplicar máscara NDWI ---
turbidez_agua = np.where(ndwi_bin == 1, turbidez, np.nan)

# --- Visualizar con escala entre percentil 5 y 95 ---
p5, p95 = np.nanpercentile(turbidez_agua, [5, 95])
#print(f"Visualización ajustada entre P5: {p5:.2f} y P95: {p95:.2f} NTU")

Mapa de turbidez

Código
plt.figure(figsize=(5.5, 5.5))
img = plt.imshow(turbidez_agua, cmap="rainbow", vmin=p5, vmax=p95)
plt.title("Turbidez 2025-06-09")
plt.colorbar(img, label="Turbidez (NTU)")
plt.axis("off")
plt.tight_layout()
plt.show()

8 Métodos de aprendizaje automático

Proximamente…

Referencias

[1]
J. cc Delegido et al., «Turbidity and secchi disc depth with sentinel-2 in different trophic status reservoirs at the comunidad valenciana», Revista de Teledeteccion, vol. 2019, pp. 15-24, 2019, doi: 10.4995/raet.2019.12603.
[2]
Y. Ma et al., «Remote sensing of turbidity for lakes in Northeast China using sentinel-2 images with machine learning algorithms», IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 14, pp. 9132-9146, 2021, doi: 10.1109/JSTARS.2021.3109292.
[3]
A. K. M. A. Hossain, C. Mathias, y R. Blanton, «Remote sensing of turbidity in the tennessee river using landsat 8 satellite», Remote Sensing, vol. 13, sep. 2021, doi: 10.3390/rs13183785.
[4]
M. Najafzadeh y S. Basirian, «Evaluation of River Water Quality Index Using Remote Sensing and Artificial Intelligence Models», Remote Sensing, vol. 15, may 2023, doi: 10.3390/rs15092359.
[5]
M. Allam, M. Y. A. Khan, y Q. Meng, «Retrieval of turbidity on a spatio-temporal scale using landsat 8 SR: A case study of the Ramganga River in the Ganges Basin, India», Applied Sciences (Switzerland), vol. 10, jun. 2020, doi: 10.3390/app10113702.
[6]
Y. Luo, D. Doxaran, y Q. Vanhellemont, «Retrieval and validation of water turbidity at metre-scale using Pleiades satellite data: A case study in the Gironde Estuary», Remote Sensing, vol. 12, mar. 2020, doi: 10.3390/rs12060946.
[7]
M. Shen, S. Wang, Y. Li, M. Tang, y Y. Ma, «Pattern of turbidity change in the middle reaches of the yarlung zangbo river, southern tibetan plateau, from 2007 to 2017», Remote Sensing, vol. 13, pp. 1-25, ene. 2021, doi: 10.3390/rs13020182.
[8]
Y. O. Ouma, K. Noor, y K. Herbert, «Modelling Reservoir Chlorophyll- a, TSS, and Turbidity Using Sentinel-2A MSI and Landsat-8 OLI Satellite Sensors with Empirical Multivariate Regression», Journal of Sensors, vol. 2020, 2020, doi: 10.1155/2020/8858408.
[9]
Z. Yang et al., «Combined Retrievals of Turbidity from Sentinel-2A/B and Landsat-8/9 in the Taihu Lake through Machine Learning», Remote Sensing, vol. 15, sep. 2023, doi: 10.3390/rs15174333.
[10]
Y. Zhao et al., «Retrieval of Water Quality Parameters Based on Near-Surface Remote Sensing and Machine Learning Algorithm», Remote Sensing, vol. 14, nov. 2022, doi: 10.3390/rs14215305.
[11]
S. Chatterjee y A. S. Hadi, Regression analysis by example. John Wiley & Sons, 2015.

Notas

  1. 0,83 - 112,26 NTU.↩︎

  2. 2,3 - 107,02 NTU.↩︎

  3. IOA = index of agreement
    SI = scatter index.↩︎

  4. 20,6 - 112 NTU
    2,3 - 15,4 NTU.↩︎

  5. MAPE = Mean Absolute Percentage Error
    0 - 1300 NTU
    0 - 80 NTU.↩︎