import pandas as pd
Algoritmo para la estimación de la turbidez sobre el Río Paraná
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.
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
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].
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
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.
= pd.read_csv(r"D:\GIT\Turbidez\datos\base_de_datos_lab.csv")
df1_lab = pd.read_csv(r"D:\GIT\Turbidez\datos\base_de_datos_gis_acolite.csv") df2_gis
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.
= pd.merge(df1_lab, df2_gis, on=['latitud', 'longitud','fecha'], how='inner') df_combinado
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_combinado[(df_combinado['param'] == 'turb')] df_turbidez
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.drop(columns=['longitud','latitud','punto','fecha','param'])
df_turbidez_banda
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
={'valor': 'turbidez'}, inplace=True)
df_turbidez_banda.rename(columns
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_turbidez_banda.pivot_table(
df_final ='turbidez',
index='banda',
columns='reflect',
values )
¿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
'Turbidez_FINAL.csv', index=True) df_final.to_csv(
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_turbidez_banda.pivot_table(
df_final ='turbidez',
index='banda',
columns='reflect',
values
)
= df_final.reset_index() #Línea de código que se debe agregar df_final
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.
'Turbidez_FINAL2.csv', index=False) df_final.to_csv(
Verificación tabla final
import pandas as pd
= pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')
df
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
= pd.read_csv(r"D:\GIT\Turbidez\Turbidez_FINAL.csv")
df
= df[['B05']]
X = df['turbidez'] y
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
= train_test_split(
X_train, X_test, y_train, y_test =0.25, shuffle=True, random_state=42
X, y, test_size )
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
= LinearRegression().fit(X_train, y_train) regressor
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.
= regressor.predict(X_test)
y_pred = mean_squared_error(y_test, y_pred)
p_rmse = r2_score(y_test, y_pred) p_r2
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)
= regressor.coef_
coef
# Intercepto (ordenada al origen)
= regressor.intercept_
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
= plt.subplots(ncols=2, figsize=(10, 5), sharex=True, sharey=True)
fig, ax
#Gráfico de entrenamiento
0].plot(
ax[
X_train,
regressor.predict(X_train),=3,
linewidth="#17A77E",
color="Modelo",
label
)
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()
ax[
#Gráfico de validación
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()
ax[
#Ecuación de la recta
= regressor.coef_[0]
coef = regressor.intercept_
intercept = f"y = {coef:.2f}x + {intercept:.2f}"
equation # Mostrar la ecuación en ambos subgráficos (opcionalmente, puedes usar solo uno)
for a in ax:
0.05, 0.95, equation, transform=a.transAxes,
a.text(=10, verticalalignment='top',
fontsize=dict(boxstyle="round", facecolor="white", alpha=0.7))
bbox
"Regresión lineal")
fig.suptitle(
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
= pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')
Datos#print (Datos.head())
Calculamos coeficiente de correlacion lineal “r” entre la turbidez y cada banda con la función .corr de pandas.
Código
= [col for col in Datos.columns if col.startswith('B')]
bandas
= {}
correlaciones
for banda in bandas:
= Datos['turbidez'].corr(Datos[banda])
r = r
correlaciones[banda]
#Creamos un Data Frame.
= pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
df_correlaciones1 #Ordenamos por mayor a menos valor de r.
= df_correlaciones1.sort_values(by='r', ascending=False).reset_index(drop=True)
df_correlaciones1
'Banda'] = "turb vs " + df_correlaciones1['Banda'].astype(str)
df_correlaciones1[={'Banda': 'Combinación 1'}, inplace=True)
df_correlaciones1.rename(columns
r'D:\GIT\Turbidez\Datos creados\Correlaciones\C1_turb_vs_banda.csv', index=False)
df_correlaciones1.to_csv(
#Creo archivo csv para usarlo para entrenar modelos
r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv", index=False)
Datos.to_csv(
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
= pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")
df
# Seleccionar las columnas de bandas (excepto turbidez)
= [col for col in df.columns if col != 'turbidez']
bandas
# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
= [col for col in df.columns if col != 'turbidez']
bandas
# Cantidad de bandas
= len(bandas)
n
# Calcular las filas y columnas para el grid de subplots
= (n + 2) // 3 # hasta 3 columnas por fila
filas = 3
columnas
# Crear figura y ejes
= plt.subplots(filas, columnas, figsize=(12, filas * 3))
fig, axes = axes.flatten()
axes
# Graficar cada banda
for i, banda in enumerate(bandas):
= axes[i]
ax 'turbidez'], alpha=0.6, s=10)
ax.scatter(df[banda], df[f'Turbidez vs {banda}', fontsize=10)
ax.set_title(
ax.set_xlabel(banda)'Turbidez')
ax.set_ylabel(True)
ax.grid(
# 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
= pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')
Datos_turb_log 'turbidez'] = np.log(Datos_turb_log['turbidez'])
Datos_turb_log[
#Cambio el nombre la columna "turbidez" luego de aplicar el logaritmo
= Datos_turb_log.rename(columns={'turbidez': 'ln_turbidez'})
Datos_turb_log
#Creo archivo csv para usarlo para entrenar modelos
r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C2_ln_turb_banda.csv", index=False)
Datos_turb_log.to_csv(#print(Datos_turb_log.head())
Calculamos r entre el ln_turb y cada banda, con la función .corr de pandas.
Código
= [col for col in Datos_turb_log.columns if col.startswith('B')]
bandas
= {}
correlaciones
for banda in bandas:
= Datos_turb_log['ln_turbidez'].corr(Datos_turb_log[banda])
r = r
correlaciones[banda]
#Creamos un Data Frame.
= pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
df_correlaciones2 #Ordenamos por mayor a menos valor de r.
= 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[={'Banda': 'Combinación 2'}, inplace=True)
df_correlaciones2.rename(columns
r'D:\GIT\Turbidez\Datos creados\Correlaciones\C2_ln_turb_vs_banda.csv', index=False)
df_correlaciones2.to_csv(
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
= pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C2_ln_turb_banda.csv")
df
# Seleccionar las columnas de bandas (excepto turbidez)
= [col for col in df.columns if col != 'ln_turbidez']
bandas
# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
= [col for col in df.columns if col != 'ln_turbidez']
bandas
# Cantidad de bandas
= len(bandas)
n
# Calcular las filas y columnas para el grid de subplots
= (n + 2) // 3 # hasta 3 columnas por fila
filas = 3
columnas
# Crear figura y ejes
= plt.subplots(filas, columnas, figsize=(12, filas * 3))
fig, axes = axes.flatten()
axes
# Graficar cada banda
for i, banda in enumerate(bandas):
= axes[i]
ax 'ln_turbidez'], alpha=0.6, s=10)
ax.scatter(df[banda], df[f'ln_turbidez vs {banda}', fontsize=10)
ax.set_title(
ax.set_xlabel(banda)'ln_turbidez')
ax.set_ylabel(True)
ax.grid(
# 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
= pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')
Datos
# Guardar cantidad original de filas
= Datos.shape[0]
filas_originales
# Filtrar filas: nos quedamos solo con aquellas donde **todos los valores sean mayores a 0**
= Datos[(Datos > 0).all(axis=1)]
Datos_filtrado
# Calcular cuántas filas se eliminaron
= Datos_filtrado.shape[0]
filas_finales = filas_originales - filas_finales
filas_eliminadas
# 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:
= Datos_filtrado
DatosC3_C4 #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_C4.copy()
DatosC3
= [col for col in DatosC3.columns if col.startswith('B')]
col
= np.log(DatosC3[col])
DatosC3[col]
#Cambiamos en nombre las columnas, agremamos ln_ a cada columna
= ['ln_' + col for col in DatosC3.columns]
DatosC3.columns
={'ln_turbidez': 'turbidez'}, inplace=True)
DatosC3.rename(columns
#Creo archivo csv para usarlo para entrenar modelos
r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C3_turb_vs_ln_banda.csv", index=False)
DatosC3.to_csv(
#DatosC3
Calculamos r entre turb y ln de cada banda, con la función .corr de pandas.
Código
= [col for col in DatosC3.columns if col.startswith('ln_B')]
bandas_ln
= {}
correlaciones
for banda in bandas_ln:
= DatosC3['turbidez'].corr(DatosC3[banda])
r = r
correlaciones[banda]
#Creamos un Data Frame.
= pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
df_correlaciones3 #Ordenamos por mayor a menos valor de r.
= df_correlaciones3.sort_values(by='r', ascending=False).reset_index(drop=True)
df_correlaciones3
'Banda'] = "turb vs " + df_correlaciones3['Banda'].astype(str)
df_correlaciones3[={'Banda': 'Combinación 3'}, inplace=True)
df_correlaciones3.rename(columns
r'D:\GIT\Turbidez\Datos creados\Correlaciones\C3_turb_vs_ln_banda.csv', index=False)
df_correlaciones3.to_csv(
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
= pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C3_turb_vs_ln_banda.csv")
df
# Seleccionar las columnas de bandas (excepto turbidez)
= [col for col in df.columns if col.startswith('ln_B')]
bandas
# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
= [col for col in df.columns if col.startswith('ln_B')]
bandas
# Cantidad de bandas
= len(bandas)
n
# Calcular las filas y columnas para el grid de subplots
= (n + 2) // 3 # hasta 3 columnas por fila
filas = 3
columnas
# Crear figura y ejes
= plt.subplots(filas, columnas, figsize=(12, filas * 3))
fig, axes = axes.flatten()
axes
# Graficar cada banda
for i, banda in enumerate(bandas):
= axes[i]
ax 'turbidez'], alpha=0.6, s=10)
ax.scatter(df[banda], df[f'Turbidez vs {banda}', fontsize=10)
ax.set_title(
ax.set_xlabel(banda)'Turbidez')
ax.set_ylabel(True)
ax.grid(
# 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
= pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')
Datos
# Guardar cantidad original de filas
= Datos.shape[0]
filas_originales
# Filtrar filas: nos quedamos solo con aquellas donde **todos los valores sean mayores a 0**
= Datos[(Datos > 0).all(axis=1)]
Datos_filtrado
# Calcular cuántas filas se eliminaron
= Datos_filtrado.shape[0]
filas_finales = filas_originales - filas_finales
filas_eliminadas
# 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:
= Datos_filtrado
DatosC3_C4 #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
= np.log(DatosC3_C4)
DatosC4 #Cambiamos en nombre las columnas, agremamos ln_ a cada columna
= ['ln_' + col for col in DatosC4.columns]
DatosC4.columns
#Creo archivo csv para usarlo para entrenar modelos
r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv", index=False)
DatosC4.to_csv(
#DatosC4.head()
Calculamos r entre el ln_turb y ln de cada banda, con la función .corr de pandas.
Código
= [col for col in DatosC4.columns if col.startswith('ln_B')]
bandas_ln
= {}
correlaciones
for banda in bandas_ln:
= DatosC4['ln_turbidez'].corr(DatosC4[banda])
r = r
correlaciones[banda]
#Creamos un Data Frame.
= pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
df_correlaciones4 #Ordenamos por mayor a menos valor de r.
= 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[={'Banda': 'Combinación 4'}, inplace=True)
df_correlaciones4.rename(columns
r'D:\GIT\Turbidez\Datos creados\Correlaciones\C4_ln_turb_vs_ln_banda.csv', index=False)
df_correlaciones4.to_csv(
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
= pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")
df
# Seleccionar las columnas de bandas (excepto turbidez)
= [col for col in df.columns if col != 'ln_turbidez']
bandas
# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
= [col for col in df.columns if col != 'ln_turbidez']
bandas
# Cantidad de bandas
= len(bandas)
n
# Calcular las filas y columnas para el grid de subplots
= (n + 2) // 3 # hasta 3 columnas por fila
filas = 3
columnas
# Crear figura y ejes
= plt.subplots(filas, columnas, figsize=(12, filas * 3))
fig, axes = axes.flatten()
axes
# Graficar cada banda
for i, banda in enumerate(bandas):
= axes[i]
ax 'ln_turbidez'], alpha=0.6, s=10)
ax.scatter(df[banda], df[f'ln_turbidez vs {banda}', fontsize=10)
ax.set_title(
ax.set_xlabel(banda)'ln_turbidez')
ax.set_ylabel(True)
ax.grid(
# 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
= pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C1_turb_vs_banda.csv")
C1 = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C2_ln_turb_vs_banda.csv")
C2= pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C3_turb_vs_ln_banda.csv")
C3 = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C4_ln_turb_vs_ln_banda.csv")
C4
= pd.concat([C1, C2 ,C3, C4], axis=1)
C_comb
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).
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
= pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")
Datos = Datos['turbidez']
y
# Variables
= ['B05']
variables_fijas = ['B06', 'B08', 'B07', 'B04', 'B8A', 'B03', 'B02', 'B01', 'B12', 'B11']
variables_a_agregar
# Resultados
= []
resultados
# Configuración de bootstrapping
= 1000
n_iteraciones_bootstrap
# Entrenamiento incremental
for i in range(len(variables_a_agregar) + 1):
= variables_fijas + variables_a_agregar[:i]
variables_usadas = Datos[variables_usadas].values
X
# División entrenamiento/test
= train_test_split(
X_train, X_test, y_train, y_test =0.25, shuffle=True, random_state=42
X, y, test_size
)
# Bootstrapping: obtener coeficientes promedio y métricas de entrenamiento
= run_bootstrap_linear_regression_analysis(
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot =n_iteraciones_bootstrap)
X_train, y_train.to_numpy(), n_iteraciones
# Modelo final con coeficientes promedio
= LinearRegression()
modelo_final = coef_prom
modelo_final.coef_ = intercept_prom
modelo_final.intercept_
# Predicción sobre testeo
= modelo_final.predict(X_test)
y_pred = r2_score(y_test, y_pred)
r2 = np.sqrt(mean_squared_error(y_test, y_pred))
rmse
# R² ajustado
= len(y_test)
n_obs = X_test.shape[1]
n_vars = 1 - (1 - r2) * (n_obs - 1) / (n_obs - n_vars - 1)
r2_ajustado
# AIC
= y_test - y_pred
residuals = np.sum(residuals ** 2)
rss = X_test.shape[1] + 1 # +1 por el intercepto
k = n_obs * np.log(rss / n_obs) + 2 * k
aic
# 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
= pd.DataFrame(resultados)
df_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
= [r["AIC"] for r in resultados]
aic_scores = [r["num_variables"] for r in resultados]
n_vars
='o', color='purple')
plt.plot(n_vars, aic_scores, marker'AIC vs Número de Variables')
plt.title('Número de variables')
plt.xlabel('AIC')
plt.ylabel(True)
plt.grid( 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
= pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")
Datos = Datos['ln_turbidez']
y
# Variables
= ['ln_B05']
variables_fijas = ['ln_B04', 'ln_B06', 'ln_B07', 'ln_B08', 'ln_B8A', 'ln_B03', 'ln_B02', 'ln_B01','ln_B11' ,'ln_B12']
variables_a_agregar
# Resultados
= []
resultados
# Configuración de bootstrapping
= 1000
n_iteraciones_bootstrap
# Entrenamiento incremental
for i in range(len(variables_a_agregar) + 1):
= variables_fijas + variables_a_agregar[:i]
variables_usadas = Datos[variables_usadas].values
X
# División entrenamiento/test
= train_test_split(
X_train, X_test, y_train, y_test =0.25, shuffle=True, random_state=42
X, y, test_size
)
# Bootstrapping: obtener coeficientes promedio y métricas de entrenamiento
= run_bootstrap_linear_regression_analysis(
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot =n_iteraciones_bootstrap)
X_train, y_train.to_numpy(), n_iteraciones
# Modelo final con coeficientes promedio
= LinearRegression()
modelo_final = coef_prom
modelo_final.coef_ = intercept_prom
modelo_final.intercept_
# Predicción sobre testeo
= modelo_final.predict(X_test)
y_pred = r2_score(y_test, y_pred)
r2 = np.sqrt(mean_squared_error(y_test, y_pred))
rmse
# Cálculo R² ajustado
= len(y_test)
n_obs = X_test.shape[1]
n_vars
if n_obs - n_vars - 1 != 0:
= 1 - (1 - r2) * (n_obs - 1) / (n_obs - n_vars - 1)
r2_ajustado else:
= np.nan # o podrías usar 0.0 si preferís
r2_ajustado
# AIC
= y_test - y_pred
residuals = np.sum(residuals ** 2)
rss = X_test.shape[1] + 1 # +1 por el intercepto
k = n_obs * np.log(rss / n_obs) + 2 * k
aic
# 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
= pd.DataFrame(resultados)
df_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
= [r["AIC"] for r in resultados]
aic_scores = [r["num_variables"] for r in resultados]
n_vars
='o', color='purple')
plt.plot(n_vars, aic_scores, marker'AIC vs Número de Variables')
plt.title('Número de variables')
plt.xlabel('AIC')
plt.ylabel(True)
plt.grid( 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
= pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")
df
# Variables predictoras
= ['B05']
variables = df[variables].values
X_completo = df['turbidez'].values
y_completo
# Separar entrenamiento y testeo
= train_test_split(
X_train, X_test, y_train, y_test =0.25, random_state=42
X_completo, y_completo, test_size
)
# Ejecutar bootstrapping
= 1000
n_iteraciones_config = \
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot =n_iteraciones_config)
run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones
# Modelo final con coeficientes promedio
= LinearRegression()
modelo_final_promedio = coef_prom
modelo_final_promedio.coef_ = intercept_prom
modelo_final_promedio.intercept_
# Predicción sobre entrenamiento y test
= modelo_final_promedio.predict(X_train)
y_train_pred = modelo_final_promedio.predict(X_test)
y_test_pred
# Métricas en testeo
= r2_score(y_test, y_test_pred)
r2_test = np.sqrt(mean_squared_error(y_test, y_test_pred))
rmse_test
# ------------------ Gráfico con estilo personalizado ------------------ #
=(5, 5))
plt.figure(figsize
# Entrenamiento
="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')
plt.scatter(y_train, y_train_pred, color
# Test
="red", alpha=0.7, label="Datos de testeo", marker='^')
plt.scatter(y_test, y_test_pred, color
# Línea ideal
= min(np.min(y_train), np.min(y_test))
min_val = max(np.max(y_train), np.max(y_test))
max_val '--', color="#17A77E", lw=2, label="Línea ideal")
plt.plot([min_val, max_val], [min_val, max_val],
# Línea de tendencia sobre test
= np.polyfit(y_test, y_test_pred, 1) # [pendiente, intercepto]
coef_linea = np.linspace(min_val, max_val, 100)
x_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
y_tendencia '-', color='black', lw=2, label="Línea de tendencia")
plt.plot(x_tendencia, y_tendencia,
"Turbidez real (NTU)")
plt.xlabel("Turbidez estimada (NTU)")
plt.ylabel(
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()True, linestyle='--', alpha=0.5)
plt.grid(
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
= " + ".join(terminos_latex)
ecuacion_latexB05 = f"turbidez = {ecuacion_latexB05} + {intercept_prom:.4f}"
ecuacion_latexB05
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
= pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")
df
# Variables predictoras
= ['B05', 'B06', 'B08', 'B07', 'B04', 'B8A']
variables = df[variables].values
X_completo = df['turbidez'].values
y_completo
# Separar entrenamiento y testeo
= train_test_split(
X_train, X_test, y_train, y_test =0.25, random_state=42
X_completo, y_completo, test_size
)
# Ejecutar bootstrapping
= 1000
n_iteraciones_config = \
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot =n_iteraciones_config)
run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones
# Modelo final con coeficientes promedio
= LinearRegression()
modelo_final_promedio = coef_prom
modelo_final_promedio.coef_ = intercept_prom
modelo_final_promedio.intercept_
# Predicción sobre entrenamiento y test
= modelo_final_promedio.predict(X_train)
y_train_pred = modelo_final_promedio.predict(X_test)
y_test_pred
# Métricas en testeo
= r2_score(y_test, y_test_pred)
r2_test = np.sqrt(mean_squared_error(y_test, y_test_pred))
rmse_test
# ------------------ Gráfico con estilo personalizado ------------------ #
=(5, 5))
plt.figure(figsize
# Entrenamiento
="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')
plt.scatter(y_train, y_train_pred, color
# Test
="red", alpha=0.7, label="Datos de testeo", marker='^')
plt.scatter(y_test, y_test_pred, color
# Línea ideal
= min(np.min(y_train), np.min(y_test))
min_val = max(np.max(y_train), np.max(y_test))
max_val '--', color="#17A77E", lw=2, label="Línea ideal")
plt.plot([min_val, max_val], [min_val, max_val],
# Línea de tendencia sobre test
= np.polyfit(y_test, y_test_pred, 1) # [pendiente, intercepto]
coef_linea = np.linspace(min_val, max_val, 100)
x_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
y_tendencia '-', color='black', lw=2, label="Línea de tendencia")
plt.plot(x_tendencia, y_tendencia,
"Turbidez real (NTU)")
plt.xlabel("Turbidez estimada (NTU)")
plt.ylabel(
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()True, linestyle='--', alpha=0.5)
plt.grid(
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
= " + ".join(terminos_latex)
ecuacion_latex = f"turbidez = {ecuacion_latex} + {intercept_prom:.4f}"
ecuacion_latex
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
= pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")
df
# Variables predictoras
= ['ln_B05']
variables = df[variables].values
X_completo = df['ln_turbidez'].values
y_completo
# Separar entrenamiento y testeo
= train_test_split(
X_train, X_test, y_train, y_test =0.25, random_state=42
X_completo, y_completo, test_size
)
# Ejecutar bootstrapping
= 1000
n_iteraciones_config = \
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot =n_iteraciones_config)
run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones
# Modelo final con coeficientes promedio
= LinearRegression()
modelo_final_promedio = coef_prom
modelo_final_promedio.coef_ = intercept_prom
modelo_final_promedio.intercept_
# Predicción sobre entrenamiento y test
= modelo_final_promedio.predict(X_train)
y_train_pred = modelo_final_promedio.predict(X_test)
y_test_pred
# Métricas en testeo
= r2_score(y_test, y_test_pred)
r2_test = np.sqrt(mean_squared_error(y_test, y_test_pred))
rmse_test
# ------------------ Gráfico con estilo personalizado ------------------ #
=(5, 5))
plt.figure(figsize
# Entrenamiento
="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')
plt.scatter(y_train, y_train_pred, color
# Test
="red", alpha=0.7, label="Datos de testeo", marker='^')
plt.scatter(y_test, y_test_pred, color
# Línea ideal
= min(np.min(y_train), np.min(y_test))
min_val = max(np.max(y_train), np.max(y_test))
max_val '--', color="#17A77E", lw=2, label="Línea ideal")
plt.plot([min_val, max_val], [min_val, max_val],
# Línea de tendencia sobre test
= np.polyfit(y_test, y_test_pred, 1) # [pendiente, intercepto]
coef_linea = np.linspace(min_val, max_val, 100)
x_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
y_tendencia '-', color='black', lw=2, label="Línea de tendencia")
plt.plot(x_tendencia, y_tendencia,
"Turbidez real (NTU)")
plt.xlabel("Turbidez estimada (NTU)")
plt.ylabel(
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()True, linestyle='--', alpha=0.5)
plt.grid(
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
= " + ".join(terminos_latex)
ecuacion_latex = f"\\ln(\\text{{turbidez}}) = {ecuacion_latex} + {intercept_prom:.4f}"
ecuacion_latex
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
= np.exp(y_train_pred)
y_train_pred_real = np.exp(y_test_pred)
y_test_pred_real
= np.exp(y_train)
y_train_real = np.exp(y_test)
y_test_real
= np.sqrt(mean_squared_error(y_test_real, y_test_pred_real))
rmse_test_real = r2_score(y_test_real, y_test_pred_real)
r2_test_real
#Grafico
=(5, 5))
plt.figure(figsize
# Entrenamiento (en escala real)
="#9D50A6", alpha=0.5, label="Entrenamiento")
plt.scatter(y_train_real, y_train_pred_real, color
# Test (en escala real)
="red", alpha=0.7, label="Test")
plt.scatter(y_test_real, y_test_pred_real, color
# Línea ideal
= min(y_train_real.min(), y_test_real.min())
min_val = max(y_train_real.max(), y_test_real.max())
max_val '--', color="#17A77E", lw=2, label="Línea ideal")
plt.plot([min_val, max_val], [min_val, max_val],
# Línea de tendencia sobre test
= np.polyfit(y_test, y_test_pred, 1) # [pendiente, intercepto]
coef_linea = np.linspace(min_val, max_val, 100)
x_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
y_tendencia '-', color='black', lw=2, label="Línea de tendencia")
plt.plot(x_tendencia, y_tendencia,
"Turbidez real (NTU)")
plt.xlabel("Turbidez estimada (NTU)")
plt.ylabel(
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()True, linestyle='--', alpha=0.5)
plt.grid(
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
= pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")
df
# Variables predictoras
= ['ln_B05','ln_B04','ln_B06']
variables = df[variables].values
X_completo = df['ln_turbidez'].values
y_completo
# Separar entrenamiento y testeo
= train_test_split(
X_train, X_test, y_train, y_test =0.25, random_state=42
X_completo, y_completo, test_size
)
# Ejecutar bootstrapping
= 1000
n_iteraciones_config = \
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot =n_iteraciones_config)
run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones
# Modelo final con coeficientes promedio
= LinearRegression()
modelo_final_promedio = coef_prom
modelo_final_promedio.coef_ = intercept_prom
modelo_final_promedio.intercept_
# Predicción sobre entrenamiento y test
= modelo_final_promedio.predict(X_train)
y_train_pred = modelo_final_promedio.predict(X_test)
y_test_pred
# Métricas en testeo
= r2_score(y_test, y_test_pred)
r2_test = np.sqrt(mean_squared_error(y_test, y_test_pred))
rmse_test
# ------------------ Gráfico con estilo personalizado ------------------ #
=(5, 5))
plt.figure(figsize
# Entrenamiento
="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')
plt.scatter(y_train, y_train_pred, color
# Test
="red", alpha=0.7, label="Datos de testeo", marker='^')
plt.scatter(y_test, y_test_pred, color
# Línea ideal
= min(np.min(y_train), np.min(y_test))
min_val = max(np.max(y_train), np.max(y_test))
max_val '--', color="#17A77E", lw=2, label="Línea ideal")
plt.plot([min_val, max_val], [min_val, max_val],
# Línea de tendencia sobre test
= np.polyfit(y_test, y_test_pred, 1) # [pendiente, intercepto]
coef_linea = np.linspace(min_val, max_val, 100)
x_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
y_tendencia '-', color='black', lw=2, label="Línea de tendencia")
plt.plot(x_tendencia, y_tendencia,
"Turbidez real (NTU)")
plt.xlabel("Turbidez estimada (NTU)")
plt.ylabel(
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()True, linestyle='--', alpha=0.5)
plt.grid(
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
= " + ".join(terminos_latex)
ecuacion_latex = f"\\ln(\\text{{turbidez}}) = {ecuacion_latex} + {intercept_prom:.4f}"
ecuacion_latex
from IPython.display import display, Math
display(Math(ecuacion_latex))
# Ecuacion real guardada como string
# Exponenciar el intercepto
= np.exp(intercept_prom)
intercepto_escala_real
# 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
= f"\\text{{turbidez}} = {intercepto_escala_real:.4f} \\cdot " + " \\cdot ".join(terminos_real)
ecuacion_real_latex
# (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
= f"turbidez = {intercepto_escala_real:.4f} * " + " * ".join([
ecuacion_turb 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
= np.exp(y_train_pred)
y_train_pred_real = np.exp(y_test_pred)
y_test_pred_real
= np.exp(y_train)
y_train_real = np.exp(y_test)
y_test_real
= np.sqrt(mean_squared_error(y_test_real, y_test_pred_real))
rmse_test_real = r2_score(y_test_real, y_test_pred_real)
r2_test_real
#Grafico
=(5, 5))
plt.figure(figsize
# Entrenamiento (en escala real)
="#9D50A6", alpha=0.5, label="Entrenamiento")
plt.scatter(y_train_real, y_train_pred_real, color
# Test (en escala real)
="red", alpha=0.7, label="Test")
plt.scatter(y_test_real, y_test_pred_real, color
# Línea ideal
= min(y_train_real.min(), y_test_real.min())
min_val = max(y_train_real.max(), y_test_real.max())
max_val '--', color="#17A77E", lw=2, label="Línea ideal")
plt.plot([min_val, max_val], [min_val, max_val],
# Línea de tendencia sobre test
= np.polyfit(y_test, y_test_pred, 1) # [pendiente, intercepto]
coef_linea = np.linspace(min_val, max_val, 100)
x_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
y_tendencia '-', color='black', lw=2, label="Línea de tendencia")
plt.plot(x_tendencia, y_tendencia,
"Turbidez real (NTU)")
plt.xlabel("Turbidez estimada (NTU)")
plt.ylabel(
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()True, linestyle='--', alpha=0.5)
plt.grid(
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
= r"D:\GIT\Turbidez\recorte_acolite\2025-06-09.tif"
ruta
# Abrir el archivo
with rs.open(ruta) as raster:
#Definimos las bandas
#Bandas RGB para ver área de estudio
= raster.read(2) #Azul
B02 = raster.read(3) #Verde
B03 = raster.read(4) #Rojo
B04
#Para calcular NDWI con B03 y B11
= raster.read(10)
B11
#para el modelo
= raster.read(5)
B05 = raster.read(6)
B06 = raster.read(7)
B07
# Composión de imagen
= np.stack([B04, B03, B02], axis=-1)
rgb
# Mejorar constraste
def stretch(banda):
= np.percentile(banda, (2, 98))
p2, p98 return np.clip((banda - p2) / (p98 - p2 + 1e-6), 0, 1)
= np.stack([stretch(B04), stretch(B03), stretch(B02)], axis=-1)
rgb_stretched
# Gráfico
plt.imshow(rgb_stretched)"Río Paraná")
plt.title( 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
= (B03 - B11) / (B03 + B11 + 1e-6)
ndwi
#+ 1e-6 para evitar división por cero
# Binarización: agua = 1, no agua = 0
= np.where(ndwi >= 0.6, 1, 0)
ndwi_bin
# Gráfico máscara de agua
='Greys')
plt.imshow(ndwi_bin, cmap"NDWI")
plt.title(="Agua (1) / No Agua (0)")
plt.colorbar(label 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 ---
= "turbidez = 1290.3786 * B05 + 10.5"
ecuacion_latexB05
# --- Extraer coeficiente e intercepto ---
= re.match(r"turbidez\s*=\s*([-\d.]+)\s*\*\s*B05\s*([\+\-]\s*[\d.]+)", ecuacion_latexB05)
match = float(match.group(1))
coef_B05 = float(match.group(2).replace(" ", "")) # Quitar espacios en el intercepto
intercepto
# --- Calcular turbidez directamente ---
= coef_B05 * B05 + intercepto
turbidez
# --- Aplicar máscara NDWI (solo donde hay agua) ---
= np.where(ndwi_bin == 1, turbidez, np.nan)
turbidez_agua
# --- Calcular percentiles para ajustar visualización ---
= np.nanpercentile(turbidez_agua, [5, 95]) p5, p95
Mapa de turbidez
Código
=(5.5, 5.5))
plt.figure(figsize= plt.imshow(turbidez_agua, cmap="rainbow", vmin=p5, vmax=p95)
img "Turbidez 2025-06-09")
plt.title(="Turbidez (NTU)")
plt.colorbar(img, label"off")
plt.axis(
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 ---
= "turbidez = 822.8952 * B05^7.6903 * B04^-5.3331 * B06^-1.1797"
ecuacion_turb
# --- Extraer intercepto y coeficientes ---
# Formato esperado: turbidez = intercepto * Bxx^coef * ...
= re.match(r"turbidez\s*=\s*([0-9.]+)\s*\*\s*(.+)", ecuacion_turb)
match = float(match.group(1))
intercepto = match.group(2)
terminos
# Obtener variables y coeficientes
= r"(B\d{2})\^(-?\d+\.\d+)"
patron = zip(*re.findall(patron, terminos))
variables, coeficientes = np.array([float(c) for c in coeficientes])
coeficientes = np.log(intercepto)
intercept_ln
# --- 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 ---
= intercept_ln
ln_turbidez for var, coef in zip(variables, coeficientes):
+= coef * np.log(bandas[var])
ln_turbidez
= np.exp(ln_turbidez)
turbidez
# --- Aplicar máscara NDWI ---
= np.where(ndwi_bin == 1, turbidez, np.nan)
turbidez_agua
# --- Visualizar con escala entre percentil 5 y 95 ---
= np.nanpercentile(turbidez_agua, [5, 95])
p5, p95 #print(f"Visualización ajustada entre P5: {p5:.2f} y P95: {p95:.2f} NTU")
Mapa de turbidez
Código
=(5.5, 5.5))
plt.figure(figsize= plt.imshow(turbidez_agua, cmap="rainbow", vmin=p5, vmax=p95)
img "Turbidez 2025-06-09")
plt.title(="Turbidez (NTU)")
plt.colorbar(img, label"off")
plt.axis(
plt.tight_layout() plt.show()
8 Métodos de aprendizaje automático
Proximamente…