Análisis de Sólidos Suspendidos

Autor/a

Afiliación

Ariadna Malena Seba

Fecha de publicación

26 de junio de 2025

Resumen

Este sitio web aborda el análisis de sólidos suspendidos totales (SST) en cuerpos de agua mediante el uso de imágenes satelitales y técnicas de aprendizaje automático. Presenta una introducción teórica sobre la importancia de los SST como indicador ambiental y luego desarrolla una parte práctica con Python, donde se integran datos espectrales y mediciones reales para aplicar modelos de regresión. Se evalúa el desempeño de estos modelos y se exploran relaciones entre bandas espectrales y SST, proponiendo mejoras metodológicas para estudios futuros.

Palabras clave

GISTAQ, UTN, FRRe, Quarto

1 Sólidos suspendidos totales, por Vera Geneyer (https://github.com/VeraGeneyer)

Los sólidos suspendidos totales (TSM): es la cantidad de materia en suspensión en el agua, que incluye plancton, minerales, arena, y microorganismos. Se determinan como el residuo no filtrable de una muestra de agua. Niveles altos (TSM) pueden reducir la transparencia del agua, limitar la luz y y transportar sustancias tóxicas, afectando la vida acuática y la calidad del agua. Este parámetro, medido mediante sensores remotos, nos da información sobre el estado físico del cuerpo de agua y están relacionados con factores como la humedad, temperatura y entre otros, que es vital para detectar riesgos al ecosistema y cumplir con las normas ambientales.

1.1 Métodos tradicionales

Tabla 1: Características principales de algoritmos tradicionales para la estimación de sólidos suspendidos.

Ecuación	Bandas (nm)	Métricas	Aguas	Plataforma	Referencia
\(-229.34 \left( \frac{B03}{B08} \right)^{3}+1001.65 \left( \frac{B03}{B08} \right)^{2}-1422.7 \left( \frac{B03}{B08} \right)+665.17\)	B03, B08	\(R^{2}\)	Embalse¹	Landsat-8	[1]
\(-244.83+40.21 \cdot B01-3.67 \cdot NDWI\)	B01, NDWI (B03, B08)	\(R^{2}\), RMSE, d	Río²	GeoEye	[2]

De acuerdo a un estudio que analizó 48 cuerpos de agua, la estimación de TSM se hizo en su mayoría por modelos lineales, siendo la banda B8A la más frecuente [3].

1.2 Métodos de aprendizaje automático

El aprendizaje automático (ML) es una rama de la inteligencia artificial cuyo objetivo es desarrollar algoritmos capaces de resolver problemas mediante el análisis de datos y la creación de funciones que describen el comportamiento de fenómenos monitoreados [4]. Los modelos de aprendizaje automático más utilizados y mencionados por los investigadores para predecir la concentración de SST son:

Bosque Aleatorio (RF) y Refuerzo Adaptativo (AdB), modelos que se destacan por su robustez ante datos complejos y ruidosos. Estos algoritmos construyen múltiples árboles de decisión que analizan las relaciones entre características como el uso del suelo o el volumen de escorrentía y los niveles de SST [5].
Redes Neuronales Artificiales (ANN), copian las redes neuronales biológicas y aprenden patrones complejos en grandes volúmenes de datos, como los niveles de SST en distintas condiciones ambientales [5],
k-Nearest Neighbors (kNN), en sus variantes de ponderación uniforme y variable, que estima el SST en función de la cercanía en características de nuevos puntos de muestreo con datos históricos [5].

El aprendizaje automático es esencial para mejorar la precisión y rapidez en el análisis de la calidad del agua, proporcionando un monitoreo más eficiente y menos costoso en comparación con los métodos tradicionales, especialmente en áreas de difícil acceso o con datos limitados.

Tabla 2: Características principales de algoritmos de aprendizaje automático para la estimación de sólidos suspendidos.

Modelo de machine learning	Software	Agua	Datos	Métricas	Referencias
Bagging y Random Forest	Programa R	Bahía	Muestreo	Prueba de normalidad multivalente Mardia-tests y Royston	[4]
Regresión lineal, LASSO, regresión de vectores de soporte (SVR), K vecinos más cercanos (KNN), bosque aleatorio (RF) y redes neuronales artificiales (ANN).	-	Lago y embalse	Sentinel-2 y UAV	\(R^{2}\)	[6]
Regresión lineal, regresión de vectores de soporte (SVR), K vecinos más cercanos (KNN), bosque aleatorio (RF) y redes neuronales artificiales (ANN).	Programa Python	Lagos	Estación de monitoreo (Sensores para cada parámetro)	\(R^{2}\), NSE y RMSE	[5]

2 Desarrollo del modelo con Python

2.1 Procesamiento de datos (Sen2Cor)

Se detalla el procedimiento técnico que implementé para procesar información ambiental georreferenciada con el objetivo de analizar el comportamiento del parámetro sólidos suspendidos (sol_sus) en una región específica (pixel 3x3). Para esto, utilicé el lenguaje Python y la biblioteca pandas, que resulta particularmente eficiente para el manejo de estructuras tabulares.

2.1.1 Carga de datos

Primero importo la biblioteca pandas, una herramienta en Python que se utiliza para manejar datos en formato tabular (como hojas de cálculo o CSVs). Se le da el alias pd por convención, para simplificar el código.

Luego cargo dos archivos CSV con la función pd.read_csv(), la cual convierte dichos archivos en objetos del tipo DataFrame, que representan tablas en memoria, que son estructuras de datos similares a tablas (parecida a una hoja de Excel). Los conjuntos de datos cargados fueron:

gis_df: contiene información geográfica (latitud, longitud, pixel, etc.).
lab_df: contiene datos de laboratorio, incluyendo el parámetro de interés sol_sus.

Verifico la carga correcta mostrando las primeras filas con la función .head(). Es útil para ver rápidamente cómo es la estructura del archivo: qué columnas hay, qué tipo de datos, si se cargó bien.

import pandas as pd  # pandas es la biblioteca para manejar datos tabulares

# Cargar los archivos de datos
gis_df = pd.read_csv('datos/base_de_datos_gis.csv')
lab_df = pd.read_csv('datos/base_de_datos_lab.csv')

# Ver las primeras filas para asegurarse de que se cargaron bien
gis_df.head(), lab_df.head()

print("Primeras filas de gis_df:")
display(gis_df.head())

print("\nPrimeras filas de lab_df:")
display(lab_df.head())

Primeras filas de gis_df:

	fecha	punto	pixel	banda	reflect	longitud	latitud
0	2023-05-11	1	1x1	B01	0.161744	-58.868047	-27.464687
1	2023-05-11	1	1x1	B02	0.170400	-58.868047	-27.464687
2	2023-05-11	1	1x1	B03	0.208800	-58.868047	-27.464687
3	2023-05-11	1	1x1	B04	0.251000	-58.868047	-27.464687
4	2023-05-11	1	1x1	B05	0.257656	-58.868047	-27.464687


Primeras filas de lab_df:

	fecha	longitud	latitud	param	valor
0	2023-05-11	-58.868047	-27.464687	ph	6.8
1	2023-05-11	-58.868047	-27.464687	cond	141.1
2	2023-05-11	-58.868047	-27.464687	sol_sus	198.0
3	2023-05-11	-58.868047	-27.464687	turb	185.5
4	2023-05-11	-58.866729	-27.466303	ph	6.8

2.1.2 Filtrar el parámetro ‘sol_sus’

En el conjunto de datos del laboratorio lab_df, hay múltiples parámetros medidos (como pH, turbidez, etc.). En este caso, me interesa trabajar solamente con los datos de sólidos suspendidos, identificado como "sol_sus" en la columna param. Este filtrado selectivo lo realicé para limitar el análisis al fenómeno físico-químico de interés.

Filtré el DataFrame para quedarme solo con esas filas, y renombré la columna valor como sol_sus para que sea más claro en los siguientes pasos.

# Filtrar solo las filas donde el parámetro es 'sol_sus'
sol_sus_df = lab_df[lab_df["param"] == "sol_sus"]

# Renombrar la columna 'valor' a 'sol_sus' para que tenga sentido en el merge
sol_sus_df = sol_sus_df.rename(columns={"valor": "sol_sus"})

# Mostrar para confirmar
sol_sus_df.head()
print("Primeras filas de sol_sus_df:")
display(sol_sus_df.head())

Primeras filas de sol_sus_df:

	fecha	longitud	latitud	param	sol_sus
2	2023-05-11	-58.868047	-27.464687	sol_sus	198.0
6	2023-05-11	-58.866729	-27.466303	sol_sus	150.0
10	2023-05-11	-58.864889	-27.468056	sol_sus	101.0
14	2023-05-11	-58.863367	-27.469240	sol_sus	95.0
18	2023-05-11	-58.862454	-27.470561	sol_sus	69.0

2.1.3 Transformar la columna banda en columnas individuales

En este paso, convierto los valores únicos de la columna banda (como B01, B02, etc.) en nombres de columnas. Cada nueva columna contendrá los valores del parámetro reflect correspondientes a esa banda en particular. Esta operación se realiza antes de unir con los valores de sol_sus, ya que el valor de reflectancia depende de la banda, mientras que sol_sus es un dato independiente que se asignará luego por punto, fecha y ubicación.

# Pivotear la tabla para que cada banda sea una columna
gis_pivot = gis_df.pivot_table(
    index=['fecha', 'punto', 'pixel', 'latitud', 'longitud'],
    columns='banda',
    values='reflect'
).reset_index()

# Eliminar el nombre del índice de columnas generado por el pivot
gis_pivot.columns.name = None

print("Primeras filas de gis_pivot:")
display(gis_pivot.head())

Primeras filas de gis_pivot:

	fecha	punto	pixel	latitud	longitud	B01	B02	B03	B04	B05	B06	B07	B08	B11	B12	B8A
0	2023-05-11	1	1x1	-27.464687	-58.868047	0.161744	0.170400	0.208800	0.251000	0.257656	0.202175	0.207500	0.194000	0.118625	0.114738	0.173069
1	2023-05-11	1	3x3	-27.464687	-58.868047	0.161290	0.171256	0.206644	0.250911	0.258058	0.202170	0.207490	0.193722	0.118892	0.115462	0.173137
2	2023-05-11	2	1x1	-27.466303	-58.866729	0.149925	0.170600	0.202100	0.243200	0.247169	0.186631	0.187512	0.175800	0.114100	0.111031	0.159263
3	2023-05-11	2	3x3	-27.466303	-58.866729	0.149906	0.170378	0.202278	0.242689	0.247115	0.186739	0.188027	0.175544	0.113992	0.111210	0.159373
4	2023-05-11	3	1x1	-27.468056	-58.864889	0.149725	0.166200	0.197100	0.223800	0.218175	0.159169	0.163344	0.150700	0.114263	0.113006	0.139662

2.1.4 Combinar datos geoespaciales y de laboratorio

Una vez que las bandas han sido transformadas en columnas, combino esta tabla con los valores de sólidos suspendidos (sol_sus) provenientes del laboratorio. La combinación se hace usando las columnas fecha, latitud y longitud, que permiten identificar los datos correspondientes a un mismo punto geográfico y temporal.

# Realizar el merge por ubicación y fecha
df_merged = pd.merge(
    gis_pivot,
    sol_sus_df[['fecha', 'latitud', 'longitud', 'sol_sus']],
    on=['fecha', 'latitud', 'longitud'],
    how='inner'
)

print("Primeras filas del DataFrame combinado:")
display(df_merged.head())

Primeras filas del DataFrame combinado:

	fecha	punto	pixel	latitud	longitud	B01	B02	B03	B04	B05	B06	B07	B08	B11	B12	B8A	sol_sus
0	2023-05-11	1	1x1	-27.464687	-58.868047	0.161744	0.170400	0.208800	0.251000	0.257656	0.202175	0.207500	0.194000	0.118625	0.114738	0.173069	198.0
1	2023-05-11	1	3x3	-27.464687	-58.868047	0.161290	0.171256	0.206644	0.250911	0.258058	0.202170	0.207490	0.193722	0.118892	0.115462	0.173137	198.0
2	2023-05-11	2	1x1	-27.466303	-58.866729	0.149925	0.170600	0.202100	0.243200	0.247169	0.186631	0.187512	0.175800	0.114100	0.111031	0.159263	150.0
3	2023-05-11	2	3x3	-27.466303	-58.866729	0.149906	0.170378	0.202278	0.242689	0.247115	0.186739	0.188027	0.175544	0.113992	0.111210	0.159373	150.0
4	2023-05-11	3	1x1	-27.468056	-58.864889	0.149725	0.166200	0.197100	0.223800	0.218175	0.159169	0.163344	0.150700	0.114263	0.113006	0.139662	101.0

2.1.5 Filtrado espacial por pixel

Luego de combinar los datos, aplico un filtrado adicional al DataFrame sobre la columna pixel para conservar únicamente las filas correspondientes al área geográfica designada como "3x3". Este paso reduce el dominio de análisis y permite concentrarse en una región de estudio concreta.

# Filtrar solo los datos del pixel 3x3
df_pixel_3x3 = df_merged[df_merged["pixel"] == "3x3"]

print("Primeras filas del pixel 3x3:")
display(df_pixel_3x3.head())

Primeras filas del pixel 3x3:

	fecha	punto	pixel	latitud	longitud	B01	B02	B03	B04	B05	B06	B07	B08	B11	B12	B8A	sol_sus
1	2023-05-11	1	3x3	-27.464687	-58.868047	0.161290	0.171256	0.206644	0.250911	0.258058	0.202170	0.207490	0.193722	0.118892	0.115462	0.173137	198.0
3	2023-05-11	2	3x3	-27.466303	-58.866729	0.149906	0.170378	0.202278	0.242689	0.247115	0.186739	0.188027	0.175544	0.113992	0.111210	0.159373	150.0
5	2023-05-11	3	3x3	-27.468056	-58.864889	0.149718	0.167256	0.197322	0.223511	0.218608	0.159292	0.163228	0.151544	0.114241	0.112931	0.139783	101.0
7	2023-05-11	4	3x3	-27.469240	-58.863367	0.148441	0.163289	0.188733	0.196167	0.183854	0.139388	0.140017	0.135400	0.115602	0.113097	0.128222	95.0
9	2023-05-11	5	3x3	-27.470561	-58.862454	0.145838	0.157733	0.178767	0.173911	0.162600	0.132140	0.133229	0.126656	0.117069	0.113992	0.124548	69.0

2.1.6 Guardar el archivo final

Finalmente, guardo el resultado como un nuevo archivo .csv dentro de la carpeta datos.

Por último, exporto el resultado a un nuevo archivo en formato .csv, mediante la función to_csv() de pandas, con el parámetro index=False para evitar que la columna de índice se incluya en el archivo de salida que pandas crea por defecto. Esto me permite utilizarlo después para visualización o análisis posterior.

# Guardar el archivo CSV dentro de la carpeta "datos"
df_pixel_3x3.to_csv('datos/datos_sol_sus_pixel_3x3.csv', index=False)

2.2 Procesamiento de datos (ACOLITE)

A continuación teniendo la misma lógica de procesamiento que aplique antes, pero usando el nuevo archivo base_de_datos_gis_acolite.csv.

import pandas as pd

# 1. Cargar archivos CSV 
gis_acol_df = pd.read_csv('datos/base_de_datos_gis_acolite.csv')
lab_df = pd.read_csv('datos/base_de_datos_lab.csv')

# 2. Filtrar solo el parámetro 'sol_sus' 
sol_sus_df = lab_df[lab_df["param"] == "sol_sus"]
sol_sus_df = sol_sus_df.rename(columns={"valor": "sol_sus"})

#  3. Pivotear bandas en columnas 
gis_pivot = gis_acol_df.pivot_table(
    index=['fecha', 'punto', 'latitud', 'longitud'],
    columns='banda',
    values='reflect'
).reset_index()

gis_pivot.columns.name = None  # Quitar la etiqueta 'banda'

# 4. Combinar reflectancias con sol_sus 
df_merged = pd.merge(
    gis_pivot,
    sol_sus_df[['fecha', 'latitud', 'longitud', 'sol_sus']],
    on=['fecha', 'latitud', 'longitud'],
    how='inner'
)

# 5. Guardar resultado final
df_merged.to_csv('datos/datos_sol_sus_acolite.csv', index=False)

2.3 Análisis de Regresión Lineal

En este análisis aplico un modelo de regresión lineal simple para estudiar la relación entre la reflectancia y los sólidos suspendidos, utilizando datos experimentales. La regresión lineal es una técnica fundamental del aprendizaje automático supervisado que nos permite predecir un valor continuo basado en una o más variables independientes. A lo largo de este documento, se explican paso a paso las acciones realizadas y los conceptos clave para comprender y replicar este análisis.

2.3.1 Importar librerías

En este paso, cargo las bibliotecas necesarias para procesar datos, ajustar modelos de regresión, evaluar su desempeño y visualizar los resultados.

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

2.3.2 Cargar datos desde un CSV

Importo el archivo .csv con los datos experimentales. Se visualizan las primeras filas para verificar que los datos se han cargado correctamente.

# Cargar el CSV
datos = pd.read_csv('datos/datos_sol_sus_acolite.csv')

# Mostrar las primeras filas para verificar
datos.head()
print("Primeras filas de datos:")
display(datos.head())

Primeras filas de datos:

	fecha	punto	latitud	longitud	B01	B02	B03	B04	B05	B06	B07	B08	B11	B12	B8A	sol_sus
0	2023-05-11	1	-27.464687	-58.868047	0.041367	0.046669	0.085630	0.128799	0.128704	0.091748	0.098365	0.083253	0.002521	0.000000	0.066612	198.0
1	2023-05-11	2	-27.466303	-58.866729	0.032235	0.047916	0.084191	0.123898	0.122481	0.079922	0.082403	0.069021	0.001669	0.000000	0.055924	150.0
2	2023-05-11	3	-27.468056	-58.864889	0.030956	0.044084	0.078444	0.104498	0.095737	0.052162	0.057475	0.044354	0.000209	0.000174	0.034767	101.0
3	2023-05-11	4	-27.469240	-58.863367	0.029646	0.040407	0.070384	0.078357	0.063937	0.032452	0.033994	0.028441	0.001279	0.000000	0.022597	95.0
4	2023-05-11	5	-27.470561	-58.862454	0.026737	0.034968	0.060701	0.056602	0.044067	0.024809	0.026684	0.019209	0.001879	0.000000	0.017930	69.0

2.3.3 Seleccionar variables y dividir en conjuntos

Selecciono las variables relevantes: B01 como variable independiente y sol_sus como variable dependiente. Luego divido el conjunto en dos subconjuntos: uno para entrenar el modelo y otro para probarlo, lo cual sirve para evaluar su capacidad de generalización.

# Selección de variables
X = datos[["B01"]]
y = datos["sol_sus"]

# División en entrenamiento y validación
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, shuffle=False)

2.3.4 Entrenar modelo de regresión lineal

En este paso se entrena un modelo de regresión lineal usando los datos de entrenamiento. El modelo aprende la relación matemática entre la reflectancia y los sólidos suspendidos.

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

2.3.5 Evaluar desempeño del modelo

Una vez entrenado el modelo, evaluo su desempeño usando métricas estadísticas. Estas nos permiten cuantificar qué tan bien el modelo predice los valores de sólidos suspendidos a partir de la reflectancia en los datos de prueba.

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

Métricas de desempeño

El error cuadrático medio es: 1846.79
El coeficiente de determinación (R²) es: -0.252

2.3.6 Visualizar el modelo

Finalmente, se visualiza gráficamente la relación entre reflectancia y sólidos suspendidos, tanto en el conjunto de entrenamiento como en el de prueba. Esto ayuda a interpretar de forma visual cómo se ajusta el modelo a los datos reales.

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

# Gráfico 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=0.6)
ax[0].set(xlabel="Reflectancia", ylabel="Sol_Sus", title="Conjunto de entrenamiento")
ax[0].legend()

# Gráfico 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=0.6)
ax[1].set(xlabel="Reflectancia", ylabel="Sol_Sus", title="Conjunto de validación")
ax[1].legend()

fig.suptitle("Regresión lineal")

plt.show()

2.4 Análisis de regresión por banda

Con el objetivo de profundizar el análisis, se evalúa la relación entre los sólidos suspendidos y cada una de las bandas espectrales disponibles de forma individual. Para ello, se entrena un modelo de regresión lineal simple por cada banda, utilizando los mismos datos experimentales. Este enfoque permite comparar el desempeño predictivo de cada banda por separado mediante métricas como el coeficiente de determinación R², su versión ajustada y el error cuadrático medio (RMSE).

Código

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
import math

# Cargar datos
datos = pd.read_csv('datos/datos_sol_sus_acolite.csv')

# Detectar columnas de bandas
bandas = [col for col in datos.columns if col.startswith("B")]

# Lista para guardar resultados
resultados = []

# Parámetros para organización de gráficos
n_bandas = len(bandas)
ncols = 3  # Número de columnas de la grilla
nrows = math.ceil(n_bandas / ncols)  # Calculamos cuántas filas se necesitan

# Crear figura
fig, axs = plt.subplots(nrows=nrows, ncols=ncols, figsize=(12, 4 * nrows))
axs = axs.flatten()  # Asegura que podamos indexarlos como una lista

for i, banda in enumerate(bandas):
    # Variables
    X = datos[[banda]]
    y = datos["sol_sus"]

    # División entrenamiento/prueba
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, shuffle=False)

    # Ajuste del modelo
    modelo = LinearRegression().fit(X_train, y_train)
    y_train_pred = modelo.predict(X_train)

    # Métricas
    r2 = modelo.score(X_train, y_train)
    n = len(y_train)
    p = X_train.shape[1]
    r2_adj = 1 - (1 - r2) * (n - 1) / (n - p - 1)
    rmse = np.sqrt(mean_squared_error(y_train, y_train_pred))

    resultados.append({
        "Banda": banda,
        "R²": round(r2, 4),
        "R²_ajustado": round(r2_adj, 4),
        "RMSE": round(rmse, 4)
    })

    # Gráfico de entrenamiento
    ax = axs[i]
    ax.scatter(X_train, y_train, alpha=0.6, color="#9D50A6", label="Entrenamiento")
    ax.plot(X_train, y_train_pred, color="#17A77E", linewidth=1.8, label="Modelo")
    ax.set_title(f'{banda}\nR²={r2:.2f}', fontsize=10)
    ax.set_xlabel('Reflectancia', fontsize=8)
    ax.set_ylabel('Sol_Sus', fontsize=8)
    ax.tick_params(axis='both', which='major', labelsize=8)
    ax.legend(fontsize=7)
    ax.grid(True)

# Eliminar ejes sobrantes
for j in range(i + 1, len(axs)):
    fig.delaxes(axs[j])

plt.suptitle("Regresiones lineales por banda (entrenamiento)", fontsize=14)
plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

# Tabla resumen
df_resultados = pd.DataFrame(resultados).sort_values("R²", ascending=False).reset_index(drop=True)
print("Tabla resumen de métricas por banda:")
display(df_resultados.style.hide(axis="index"))

Tabla resumen de métricas por banda:

Banda	R²	R²_ajustado	RMSE
B05	0.182800	0.167100	33.457100
B06	0.179600	0.163800	33.523000
B07	0.173300	0.157400	33.650800
B08	0.146500	0.130100	34.191400
B8A	0.134600	0.118000	34.429300
B04	0.093900	0.076400	35.230300
B12	0.037900	0.019400	36.302600
B11	0.033300	0.014700	36.389200
B01	0.007800	-0.011200	36.864800
B03	0.005000	-0.014100	36.916700
B02	0.001400	-0.017800	36.984700

Conclusiones del análisis por banda

Al aplicar regresión lineal simple a cada banda espectral por separado, se observó que ninguna logra predecir con precisión los sólidos suspendidos. Las bandas B05, B06 y B07 fueron las que mejor se desempeñaron, pero con un R² apenas superior al 18 %, lo que indica que explican muy poca variabilidad. Además, los errores (RMSE) siguen siendo altos en relación con los valores reales.

Esto sugiere que la relación entre reflectancia y sólidos suspendidos no es lineal en escala natural, y que sería necesario probar modelos más complejos. En próximos pasos, pienso aplicar transformaciones logarítmicas y combinaciones de bandas para mejorar el ajuste.

2.5 Análisis de regresión por banda aplicando logaritmo a las variables

En esta etapa del análisis, aplico una transformación logarítmica natural a las variables de reflectancia y sólidos suspendidos antes de ajustar los modelos de regresión lineal. Esta transformación es útil para:

Estabilizar la varianza y reducir la heterocedasticidad.
Linealizar relaciones no lineales entre variables.
Evitar que valores extremos influyan excesivamente en el modelo.

El procedimiento es similar al análisis anterior, pero antes de entrenar el modelo se aplica log(x) a las columnas correspondientes. Para evitar problemas con ceros en los datos, estos se reemplazan por NaN y se eliminan del análisis.

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
import math

# 1. Cargar datos 
datos = pd.read_csv("datos/datos_sol_sus_acolite.csv")
bandas = [c for c in datos.columns if c.startswith("B")]

# 2.  Filtrar positivos (requisito para log) 
mask_pos = (datos["sol_sus"] > 0) & (datos[bandas] > 0).all(axis=1)
datos_log = datos.loc[mask_pos].copy()

# 3.Transformar a log-natural (para entrenar) 
datos_log["log_sol_sus"] = np.log(datos_log["sol_sus"])
for b in bandas:
    datos_log[f"log_{b}"] = np.log(datos_log[b])

log_bandas = [f"log_{b}" for b in bandas]

# 4. Split train / test (test se descarta aquí) 
X_total = datos_log[log_bandas]
y_total = datos_log["log_sol_sus"]
X_train, _, y_train, _ = train_test_split(
    X_total, y_total, test_size=0.2, shuffle=False
)

# 5. Funciones de métricas 
def metrics_log(model, X, y):
    """Métricas sobre los datos en log (las que ya tenías)."""
    y_pred = model.predict(X)
    rmse = np.sqrt(mean_squared_error(y, y_pred))
    r2   = r2_score(y, y_pred)
    n, p = X.shape
    r2_adj = 1 - ((1 - r2) * (n - 1)) / (n - p - 1)
    return rmse, r2, r2_adj

def metrics_original(model, X, y):
    """
    Idem, pero después de volver a la escala original:
       sol_sus_pred  = exp(ŷ_log)
       sol_sus_true  = exp(y_log)
    """
    y_pred_log  = model.predict(X).ravel()
    y_pred_orig = np.exp(y_pred_log)
    y_true_orig = np.exp(y.values)
    rmse = np.sqrt(mean_squared_error(y_true_orig, y_pred_orig))
    r2   = r2_score(y_true_orig, y_pred_orig)
    n, p = X.shape
    r2_adj = 1 - ((1 - r2) * (n - 1)) / (n - p - 1)
    return rmse, r2, r2_adj

# 6.  Entrenar un modelo por banda y registrar métricas 
resultados = []           # para log
resultados_nolog = []     # para escala original
lineas_log = {}           # (x, ŷ_log) para graficar en log-log
lineas_nolog = {}         # (x_original, ŷ_original) para graficar sin log

for b, lb in zip(bandas, log_bandas):
    Xtr = X_train[[lb]]
    reg = LinearRegression().fit(Xtr, y_train)

    # ── Métricas en log
    rmse_l, r2_l, r2a_l = metrics_log(reg, Xtr, y_train)
    resultados.append({"Banda": b, "RMSE_log": rmse_l,
                       "R²_log": r2_l, "R²aj_log": r2a_l})

    # ── Métricas en original
    rmse_o, r2_o, r2a_o = metrics_original(reg, Xtr, y_train)
    resultados_nolog.append({"Banda": b, "RMSE": rmse_o,
                             "R²": r2_o, "R²aj": r2a_o})

    # ── Recta ordenada para ambas escalas
    order = np.argsort(Xtr.values.ravel())
    x_log_ord   = Xtr.values.ravel()[order]             # log(B)
    y_pred_log  = reg.predict(Xtr).ravel()[order]       # log(sol_suŝ)
    lineas_log[b] = (x_log_ord, y_pred_log)

    # Para volver a la escala original:
    x_orig_ord  = np.exp(x_log_ord)                     # B
    y_pred_orig = np.exp(y_pred_log)                    # sol_suŝ
    lineas_nolog[b] = (x_orig_ord, y_pred_orig)

# 7. DataFrames con métricas
res_df_log    = (pd.DataFrame(resultados)
                   .round(3)
                   .sort_values("RMSE_log")
                   .reset_index(drop=True))

res_df_nolog  = (pd.DataFrame(resultados_nolog)
                   .round(3)
                   .sort_values("RMSE")
                   .reset_index(drop=True))

# 8.Gráficos *sin* log (escala original)
n_cols = 3
n_rows = math.ceil(len(bandas) / n_cols)
fig, axes = plt.subplots(n_rows, n_cols, figsize=(15, 4 * n_rows))
axes = axes.flatten()

for idx, b in enumerate(bandas):
    ax = axes[idx]
    # Datos originales de entrenamiento (sin log):
    x_orig = np.exp(X_train[f"log_{b}"])
    y_orig = np.exp(y_train)
    ax.scatter(x_orig, y_orig, alpha=0.6, color="#9D50A6", label="Entrenamiento")
    # Curva modelo en original:
    ax.plot(*lineas_nolog[b], lw=3, color="#17A77E", label="Modelo")
    ax.set(xlabel=b, ylabel="sol_sus", title=f"Banda {b} (escala original)")
    ax.legend()
    ax.grid(True)

# Si sobran ejes, quitarlos
for j in range(len(bandas), len(axes)):
    fig.delaxes(axes[j])

plt.tight_layout()
plt.show()

# 9.  Mostrar tablas de métricas
print("\nMétricas en **escala log-log** (ordenadas por RMSE_log):\n")
print(res_df_log.to_string(index=False))

print("\nMétricas en **escala original** (ordenadas por RMSE):\n")
print(res_df_nolog.to_string(index=False))


Métricas en **escala log-log** (ordenadas por RMSE_log):

Banda  RMSE_log  R²_log  R²aj_log
  B06     0.318   0.304     0.281
  B07     0.319   0.303     0.280
  B05     0.326   0.272     0.248
  B08     0.327   0.267     0.243
  B8A     0.328   0.262     0.237
  B04     0.349   0.166     0.138
  B12     0.360   0.111     0.082
  B11     0.369   0.065     0.034
  B03     0.378   0.017    -0.016
  B01     0.379   0.015    -0.018
  B02     0.380   0.007    -0.026

Métricas en **escala original** (ordenadas por RMSE):

Banda   RMSE     R²   R²aj
  B06 31.661  0.198  0.171
  B07 31.897  0.186  0.159
  B05 32.049  0.178  0.150
  B8A 32.472  0.156  0.128
  B08 32.567  0.151  0.123
  B04 33.962  0.077  0.046
  B12 33.986  0.076  0.045
  B11 34.576  0.043  0.011
  B01 35.843 -0.028 -0.063
  B03 35.940 -0.034 -0.068
  B02 35.969 -0.035 -0.070

Conclusión del ajuste log–log por banda

Aplicar la transformación logarítmica mejoró notablemente el ajuste respecto al análisis en escala natural. Las bandas B06, B07 y B05 mostraron los mejores resultados, con R²_log cercanos a 0.30 y una reducción del error (RMSE) tanto en escala logarítmica como en la original.

Aun así, el desempeño sigue siendo limitado si se usa cada banda por separado. Esto sugiere que para mejorar las predicciones es necesario combinar bandas, aplicar modelos más complejos o probar con índices espectrales derivados.

2.6 Selección de bandas con AIC y desempeño individual

Para evaluar qué bandas individuales son más relevantes para predecir los sólidos suspendidos, se ajustó un modelo de regresión lineal simple en escala logarítmica para cada banda. Además de calcular el coeficiente de determinación (R²) y el error cuadrático medio (RMSE), se calcula ahora el AIC (Criterio de Información de Akaike), que penaliza la complejidad del modelo. Esto ayuda a decidir si la incorporación de una banda es realmente valiosa o si se está agregando complejidad innecesaria al modelo.

Código

# 0. Librerías
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score

# 1. Cargar datos y detectar bandas
datos = pd.read_csv("datos/datos_sol_sus_acolite.csv")
bandas = [c for c in datos.columns if c.startswith("B")]

# 2. Filtrar casos con valores positivos (requisito de log)
mask_pos = (datos["sol_sus"] > 0) & (datos[bandas] > 0).all(axis=1)
datos_log = datos.loc[mask_pos].copy()

# 3. Transformar a log-natural
datos_log["log_sol_sus"] = np.log(datos_log["sol_sus"])
for b in bandas:
    datos_log[f"log_{b}"] = np.log(datos_log[b])

log_bandas = [f"log_{b}" for b in bandas]

# 4. Partición entrenamiento / test 
X_total = datos_log[log_bandas]
y_total = datos_log["log_sol_sus"]

X_train, _, y_train, _ = train_test_split(
    X_total, y_total, test_size=0.20, shuffle=False
)

# 5. Funciones auxiliares (AIC + métricas)
def aic_gauss(model, X, y):
    """
    AIC para regresión lineal con ruido ~ N(0, σ²)
    AIC = n · ln(RSS/n) + 2k
    """
    n = len(y)
    k = X.shape[1] + 1  # predictores + intercepto
    rss = np.sum((y - model.predict(X)) ** 2)
    return n * np.log(rss / n) + 2 * k

def metrics_log(model, X, y):
    y_pred = model.predict(X)
    rmse = np.sqrt(mean_squared_error(y, y_pred))
    r2   = r2_score(y, y_pred)
    n, p = X.shape
    r2_adj = 1 - ((1 - r2) * (n - 1)) / (n - p - 1)
    return rmse, r2, r2_adj

# 6. Ajuste univariable banda-a-banda + AIC
resultados = []

for b in bandas:
    lb = f"log_{b}"
    Xtr = X_train[[lb]]
    reg = LinearRegression().fit(Xtr, y_train)

    rmse_log, r2_log, r2a_log = metrics_log(reg, Xtr, y_train)
    aic_val = aic_gauss(reg, Xtr, y_train)

    resultados.append({
        "Banda"    : b,
        "RMSE_log" : rmse_log,
        "R²_log"   : r2_log,
        "R²aj_log" : r2a_log,
        "AIC"      : aic_val
    })

# 7. Tabla ordenada por AIC
df_log  = (pd.DataFrame(resultados)
             .round(3)
             .sort_values("AIC")
             .reset_index(drop=True))

print("\nAIC y métricas en escala log-log (ordenadas por AIC):\n")
print(df_log.to_string(index=False))


AIC y métricas en escala log-log (ordenadas por AIC):

Banda  RMSE_log  R²_log  R²aj_log     AIC
  B06     0.318   0.304     0.281 -69.242
  B07     0.319   0.303     0.280 -69.191
  B05     0.326   0.272     0.248 -67.798
  B08     0.327   0.267     0.243 -67.584
  B8A     0.328   0.262     0.237 -67.355
  B04     0.349   0.166     0.138 -63.448
  B12     0.360   0.111     0.082 -61.415
  B11     0.369   0.065     0.034 -59.781
  B03     0.378   0.017    -0.016 -58.184
  B01     0.379   0.015    -0.018 -58.112
  B02     0.380   0.007    -0.026 -57.878

Conclusión sobre la selección de bandas con AIC

Los resultados del AIC coinciden con las métricas de desempeño: las bandas B06 y B07 son las más relevantes para predecir sólidos suspendidos, mostrando los valores más bajos de AIC (-69.2 y -69.1 respectivamente), así como los mejores valores de R²_log (~0.30) y RMSE_log (~0.32).

Las demás bandas presentan AIC menos favorables, lo que indica que, individualmente, aportan menos información útil y no justifican la complejidad del modelo.

Por lo tanto, para modelos simples basados en una sola banda, B06 y B07 son las candidatas preferentes. Sin embargo, para un mejor desempeño se recomienda evaluar modelos que combinen varias bandas.

2.6.1 Incorporar validación cruzada en la selección de modelo

En esta etapa incorporo una validación cruzada (CV) para evaluar la robustez y generalización del modelo univariable basado en cada banda espectral. La validación cruzada consiste en dividir el conjunto de entrenamiento en varios “folds” o particiones, entrenar el modelo en todas menos una, y validar en la partición restante, repitiendo este proceso rotando las particiones. Esto ayuda a obtener métricas de desempeño más confiables y evitar sobreajuste al conjunto de entrenamiento.

El flujo completo que implemento es: - Dividir el conjunto completo en un 80 % para entrenamiento y un 20 % para prueba final (hold-out). - En el 80 % de entrenamiento, aplicar validación cruzada K-Fold (5 folds) para evaluar RMSE y coeficiente de determinación ajustado (R² ajustado) de cada banda. - Calcular también el AIC para cada modelo univariable, utilizando el conjunto completo de entrenamiento (80 %) para penalizar complejidad. - Ordenar las bandas por AIC y, en caso de empate, desempatar usando el promedio de RMSE de la validación cruzada. - Finalmente, con la banda ganadora entrenar el modelo sobre todo el 80 % y evaluar su desempeño en el conjunto de prueba (hold-out 20 %).

Esta estrategia permite seleccionar la banda más relevante no solo por la bondad del ajuste, sino también considerando la estabilidad y capacidad predictiva en datos no vistos, ofreciendo una selección más robusta y realista para modelos futuros.

Código

# ─────────────────────────────────────────────────────────────────────
# 0. LIBRERÍAS
# ─────────────────────────────────────────────────────────────────────
import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split, KFold
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.utils import resample

# ─────────────────────────────────────────────────────────────────────
# 1. CARGAR DATOS Y DETECTAR BANDAS
# ─────────────────────────────────────────────────────────────────────
datos  = pd.read_csv("datos/datos_sol_sus_acolite.csv")
bandas = [c for c in datos.columns if c.startswith("B")]

# ─────────────────────────────────────────────────────────────────────
# 2. FILTRAR CASOS > 0  (requisito de log)
# ─────────────────────────────────────────────────────────────────────
mask_pos = (datos["sol_sus"] > 0) & (datos[bandas] > 0).all(axis=1)
datos_log = datos.loc[mask_pos].copy()

# ─────────────────────────────────────────────────────────────────────
# 3. TRANSFORMACIÓN LOGARÍTMICA
# ─────────────────────────────────────────────────────────────────────
datos_log["log_sol_sus"] = np.log(datos_log["sol_sus"])
for b in bandas:
    datos_log[f"log_{b}"] = np.log(datos_log[b])

log_bandas = [f"log_{b}" for b in bandas]

X_full = datos_log[log_bandas].values
y_full = datos_log["log_sol_sus"].values

# ─────────────────────────────────────────────────────────────────────
# 4. SPLIT HOLD-OUT 80/20  (shuffle=False para respetar orden)
# ─────────────────────────────────────────────────────────────────────
X_train_full, X_test_full, y_train_full, y_test_full = train_test_split(
    X_full, y_full, test_size=0.20, shuffle=False
)

# ─────────────────────────────────────────────────────────────────────
# 5. FUNCIONES AUXILIARES
# ─────────────────────────────────────────────────────────────────────
def aic_gauss(model, X, y):
    """AIC con supuesto de ruido gaussiano homocedástico."""
    n = len(y)
    k = X.shape[1] + 1    # +1 intercepto
    rss = np.sum((y - model.predict(X)) ** 2)
    return n * np.log(rss / n) + 2 * k

def r2_adj(r2, n, p):
    return 1 - (1 - r2) * (n - 1) / (n - p - 1)

def cv_metrics(X_mat, y_vec, cv):
    """Devuelve medias y sd del RMSE y medias de R² y R²aj."""
    model = LinearRegression()
    rmse_folds, r2_folds, r2a_folds = [], [], []

    for tr_idx, val_idx in cv.split(X_mat):
        X_tr, X_val = X_mat[tr_idx], X_mat[val_idx]
        y_tr, y_val = y_vec[tr_idx], y_vec[val_idx]

        model.fit(X_tr, y_tr)
        y_pred = model.predict(X_val)

        rmse = np.sqrt(mean_squared_error(y_val, y_pred))
        r2   = r2_score(y_val, y_pred)
        r2a  = r2_adj(r2, len(y_val), X_mat.shape[1])

        rmse_folds.append(rmse)
        r2_folds.append(r2)
        r2a_folds.append(r2a)

    return {
        "rmse_mean": np.mean(rmse_folds),
        "rmse_std":  np.std(rmse_folds, ddof=1),
        "r2_mean":   np.mean(r2_folds),
        "r2a_mean":  np.mean(r2a_folds)
    }

# Configuración de la validación cruzada interior
kf = KFold(n_splits=5, shuffle=True, random_state=42)

# ─────────────────────────────────────────────────────────────────────
# 6. EVALUACIÓN UNIVARIABLE  (CV + AIC)
# ─────────────────────────────────────────────────────────────────────
resultados = []

for idx_band, b in enumerate(bandas):
    X_train_band = X_train_full[:, idx_band].reshape(-1, 1)

    # 6.1  CV
    stats_cv = cv_metrics(X_train_band, y_train_full, kf)

    # 6.2  AIC en todo el 80 %
    modelo_tmp = LinearRegression().fit(X_train_band, y_train_full)
    aic_val    = aic_gauss(modelo_tmp, X_train_band, y_train_full)

    resultados.append({
        "Banda":     b,
        "AIC_train": aic_val,
        "RMSE_CV":   stats_cv["rmse_mean"],
        "CV_std":    stats_cv["rmse_std"],
        "R2_CV":     stats_cv["r2_mean"],
        "R2aj_CV":   stats_cv["r2a_mean"]
    })

df_res = (
    pd.DataFrame(resultados)
      #   PRIMER CRITERIO: AIC  ▸  SEGUNDO: RMSE_CV
      .sort_values(["AIC_train", "RMSE_CV"])
      .reset_index(drop=True)
      .round(4)
)

print("\nComparación banda-a-banda (ordenada por AIC, luego RMSE_CV):\n")
print(df_res.to_string(index=False))

# ─────────────────────────────────────────────────────────────────────
# 7. BANDA GANADORA  →  EVALUACIÓN HOLD-OUT
# ─────────────────────────────────────────────────────────────────────
best_band = df_res.loc[0, "Banda"]
idx_best  = bandas.index(best_band)

X_train_best = X_train_full[:, idx_best].reshape(-1, 1)
X_test_best  = X_test_full[:,  idx_best].reshape(-1, 1)

best_model = LinearRegression().fit(X_train_best, y_train_full)
y_pred_test = best_model.predict(X_test_best)

rmse_test = np.sqrt(mean_squared_error(y_test_full, y_pred_test))
r2_test   = r2_score(y_test_full, y_pred_test)

print(f"\n────────  HOLD-OUT (20 %)  ─────────")
print(f"Banda ganadora: {best_band}")
print(f"RMSE_test  = {rmse_test:.3f}")
print(f"R²_test    = {r2_test:.3f}")


Comparación banda-a-banda (ordenada por AIC, luego RMSE_CV):

Banda  AIC_train  RMSE_CV  CV_std   R2_CV  R2aj_CV
  B06   -69.2419   0.3144  0.1142 -0.3386  -0.6543
  B07   -69.1908   0.3132  0.1158 -0.2749  -0.5746
  B05   -67.7976   0.3143  0.1339 -0.2391  -0.5310
  B08   -67.5837   0.3199  0.1208 -0.3510  -0.6686
  B8A   -67.3553   0.3274  0.1073 -0.4809  -0.8292
  B04   -63.4475   0.3346  0.1433 -0.3804  -0.7053
  B12   -61.4147   0.3823  0.1042 -1.2038  -1.7221
  B11   -59.7815   0.3882  0.1023 -1.1585  -1.6664
  B03   -58.1842   0.3868  0.1226 -1.0592  -1.5496
  B01   -58.1117   0.3873  0.1124 -1.0770  -1.5704
  B02   -57.8776   0.3903  0.1146 -1.1355  -1.6447

────────  HOLD-OUT (20 %)  ─────────
Banda ganadora: B06
RMSE_test  = 0.491
R²_test    = -0.291

Conclusión sobre la selección de bandas con AIC

Al incorporar validación cruzada para evaluar la selección univariable por banda en escala logarítmica, la banda B06 se mantiene como la mejor candidata, con el AIC más bajo y el RMSE promedio más pequeño en CV (aprox. 0.31).

Sin embargo, los valores de R² ajustado en CV son negativos para todas las bandas, lo que indica que el modelo lineal simple no logra explicar bien la variabilidad fuera de la muestra, sugiriendo un ajuste pobre o posible falta de linealidad real en los datos.

La evaluación final en el conjunto hold-out (20 %) confirma esta baja capacidad predictiva, con un RMSE_test de 0.491 y un R²_test negativo (-0.291), lo que indica que el modelo predice peor que un modelo que simplemente promedie los datos.

En resumen, aunque B06 es la banda más prometedora según los criterios AIC y RMSE en entrenamiento y validación interna, la baja performance en datos nuevos apunta a la necesidad de modelos más complejos, combinaciones multivariables o transformaciones adicionales para mejorar la predicción.

2.7 Selección de bandas con Forward Selection basado en AIC

En esta etapa realizo una selección de variables (bandas espectrales) combinadas utilizando un método de Forward Selection guiado por el criterio de información de Akaike (AIC). Este método comienza con un modelo vacío e incorpora progresivamente la banda que produce la mayor mejora (reducción) en el AIC, deteniéndose cuando no se consigue una mejora significativa.

La selección se realiza sobre datos transformados en escala logarítmica para estabilizar varianzas y linealizar relaciones. Además, se valida el desempeño del modelo en cada paso mediante validación cruzada (5 folds), evaluando métricas como RMSE y R² en escala original para tener una idea más realista del ajuste.

Finalmente, se evalúa el modelo seleccionado en un conjunto independiente de test (20 % de los datos) para estimar su capacidad predictiva fuera de muestra. También se presenta la ecuación final en escala log-log y en escala original para facilitar la interpretación.

Código

# ─────────────────────────────────────────────────────────────────────
# 0 · LIBRERÍAS
# ─────────────────────────────────────────────────────────────────────
import pandas as pd
import numpy as np
import math
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split, KFold
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

# ─────────────────────────────────────────────────────────────────────
# 1 · CARGAR Y PRE-PROCESAR
# ─────────────────────────────────────────────────────────────────────
df     = pd.read_csv("datos/datos_sol_sus_acolite.csv")
bandas = [c for c in df.columns if c.startswith("B")]

mask   = (df["sol_sus"] > 0) & (df[bandas] > 0).all(axis=1)
df     = df.loc[mask].copy()

df["log_sol_sus"] = np.log(df["sol_sus"])
for b in bandas:
    df[f"log_{b}"] = np.log(df[b])

log_cols = [f"log_{b}" for b in bandas]
X_full   = df[log_cols].values
y_full   = df["log_sol_sus"].values          # todavía en log

# ─────────────────────────────────────────────────────────────────────
# 2 · HOLD-OUT 80 / 20  (orden intacto → shuffle=False)
# ─────────────────────────────────────────────────────────────────────
X_tr, X_te, y_tr, y_te = train_test_split(
    X_full, y_full, test_size=0.20, shuffle=False
)

# ─────────────────────────────────────────────────────────────────────
# 3 · FUNCIONES AUXILIARES
# ─────────────────────────────────────────────────────────────────────
def aic_gauss(model, X, y):
    n = len(y)
    k = X.shape[1] + 1
    rss = np.sum((y - model.predict(X))**2)
    return n * np.log(rss / n) + 2 * k

def r2_adj(r2, n, p):
    return 1 - (1 - r2) * (n - 1) / (n - p - 1)

def cv_metrics(Xm, ym, cv):
    """
    Devuelve métricas en ESCALA ORIGINAL para 5-fold CV:
      · rmse_mean, rmse_std, r2_mean, r2a_mean
    """
    model = LinearRegression()
    rmse_f, r2_f = [], []

    for tr_idx, vl_idx in cv.split(Xm):
        model.fit(Xm[tr_idx], ym[tr_idx])          # log-modelo
        pred_log = model.predict(Xm[vl_idx])

        # → escala original
        y_val_orig   = np.exp(ym[vl_idx])
        pred_orig    = np.exp(pred_log)

        rmse_f.append(np.sqrt(mean_squared_error(y_val_orig, pred_orig)))
        r2_f  .append(r2_score(y_val_orig, pred_orig))

    r2_mean = np.mean(r2_f)
    r2a_mean = r2_adj(r2_mean, len(ym), Xm.shape[1])
    return {
        "rmse_mean": np.mean(rmse_f),
        "rmse_std":  np.std(rmse_f, ddof=1),
        "r2_mean":   r2_mean,
        "r2a_mean":  r2a_mean
    }

kf = KFold(n_splits=5, shuffle=True, random_state=42)

# ─────────────────────────────────────────────────────────────────────
# 4 · FORWARD-SELECTION GUIADO POR AIC
# ─────────────────────────────────────────────────────────────────────
remain, select, best_aic, hist = list(range(len(bandas))), [], np.inf, []

while remain:
    best_cand = None
    for idx in remain:
        cols = select + [idx]
        mdl  = LinearRegression().fit(X_tr[:, cols], y_tr)
        aic  = aic_gauss(mdl, X_tr[:, cols], y_tr)
        if best_cand is None or aic < best_cand[0]:
            best_cand = (aic, idx)

    # mejora ≥ 2 unidades
    if best_cand[0] + 2 < best_aic:
        best_aic = best_cand[0]
        select.append(best_cand[1])
        remain.remove(best_cand[1])

        stats = cv_metrics(X_tr[:, select], y_tr, kf)
        hist.append({
            "Paso": len(select),
            "Añadido": bandas[best_cand[1]],
            "AIC": round(best_aic, 3),
            "RMSE_CV": round(stats["rmse_mean"], 3),
            "CV_std":  round(stats["rmse_std"], 3),
            "R2_CV":   round(stats["r2_mean"], 3),
            "R2aj_CV": round(stats["r2a_mean"], 3)
        })
    else:
        break

hist_df = pd.DataFrame(hist)
print("\nHistorial forward-AIC:")
print(hist_df)

sel_b = [bandas[i] for i in select]
print(f"\nBandas finales seleccionadas: {sel_b}")

# ─────────────────────────────────────────────────────────────────────
# 5 · AJUSTE FINAL Y MÉTRICAS TEST (ORIGINAL)
# ─────────────────────────────────────────────────────────────────────
Xtr_sel = X_tr[:, select]
Xte_sel = X_te[:, select]

final_model = LinearRegression().fit(Xtr_sel, y_tr)

# → escala original para test
pred_test_orig = np.exp(final_model.predict(Xte_sel))
y_test_orig    = np.exp(y_te)

rmse_test = np.sqrt(mean_squared_error(y_test_orig, pred_test_orig))
mae_test  = mean_absolute_error(y_test_orig, pred_test_orig)
r2_test   = r2_score(y_test_orig, pred_test_orig)

# ─────────────────────────────────────────────────────────────────────
# 6 · ECUACIÓN FINAL (formato fracción cuando hay exponentes negativos)
# ─────────────────────────────────────────────────────────────────────
interc  = final_model.intercept_
coefs   = final_model.coef_

pos_terms = [f"{b}^{c:.4f}"      for b, c in zip(sel_b, coefs) if c >= 0]
neg_terms = [f"{b}^{-c:.4f}"     for b, c in zip(sel_b, coefs) if c < 0]

k_factor  = math.exp(interc)
numerator = " * ".join(pos_terms) if pos_terms else "1"
denominator = " * ".join(neg_terms)

print("\nEcuación en escala log-log:")
log_eq = [f"{interc:+.4f}"] + [
    f"{coefs[i]:+0.4f}·log({sel_b[i]})" for i in range(len(sel_b))
]
print("log(sol_sus) = " + " ".join(log_eq))

print("\nEscala original:")
if denominator:
    print(f"sol_sus_orig = {k_factor:.2f} · ({numerator}) / ({denominator})")
else:
    print(f"sol_sus_orig = {k_factor:.2f} · {numerator}")

# ─────────────────────────────────────────────────────────────────────
# 7 · TABLA FINAL DE MÉTRICAS
# ─────────────────────────────────────────────────────────────────────
final_row = {
    "Bandas": "+".join(sel_b),
    "AIC_train": round(best_aic, 3),
    "RMSE_CV":   round(hist[-1]["RMSE_CV"], 3),
    "CV_std":    round(hist[-1]["CV_std"], 3),
    "R²_CV":     round(hist[-1]["R2_CV"], 3),
    "R²aj_CV":   round(hist[-1]["R2aj_CV"], 3),
    "RMSE_test": round(rmse_test, 3),
    "MAE_test":  round(mae_test, 3),
    "R²_test":   round(r2_test, 3)
}

print("\n──────── TABLA FINAL ────────")
print(pd.DataFrame([final_row]).to_string(index=False))

# ─────────────────────────────────────────────────────────────────────
# 8 · GRÁFICO REAL vs PREDICHO (train + test)
# ─────────────────────────────────────────────────────────────────────
# Train → original
pred_train_orig = np.exp(final_model.predict(Xtr_sel))
y_train_orig    = np.exp(y_tr)

plt.figure(figsize=(6, 6))
plt.scatter(y_train_orig, pred_train_orig,
            alpha=0.7, label="Train")
plt.scatter(y_test_orig, pred_test_orig,
            alpha=0.7, marker='s', edgecolors='k', label="Test")

# línea y = x
mn, mx = min(y_train_orig.min(), y_test_orig.min()), \
         max(y_train_orig.max(), y_test_orig.max())
plt.plot([mn, mx], [mn, mx], "--", color="gray")

plt.xlabel("Sólidos suspendidos observados")
plt.ylabel("Sólidos suspendidos predichos")
plt.title("Ajuste modelo (escala original)")

# ecuación en caja
eq_line = rf"$\hat{{y}} = {k_factor:.1f}"
for b, c in zip(sel_b, coefs):
    eq_line += rf"\times {b}^{{{c:.2f}}}"
eq_line += "$"
plt.text(0.05, 0.95, eq_line,
         transform=plt.gca().transAxes,
         fontsize=8, va="top",
         bbox=dict(boxstyle="round,pad=0.3", fc="white", ec="gray", lw=0.5))

plt.legend()
plt.tight_layout()
plt.show()


Historial forward-AIC:
   Paso Añadido      AIC  RMSE_CV  CV_std  R2_CV  R2aj_CV
0     1     B06  -69.242   28.935  16.235 -0.717   -0.774
1     2     B02 -114.534   13.334   3.183  0.397    0.356

Bandas finales seleccionadas: ['B06', 'B02']

Ecuación en escala log-log:
log(sol_sus) = +3.6692 +0.7513·log(B06) -0.9972·log(B02)

Escala original:
sol_sus_orig = 39.22 · (B06^0.7513) / (B02^0.9972)

──────── TABLA FINAL ────────
 Bandas  AIC_train  RMSE_CV  CV_std  R²_CV  R²aj_CV  RMSE_test  MAE_test  R²_test
B06+B02   -114.534   13.334   3.183  0.397    0.356     16.693    12.778    0.859

Conclusión del modelo seleccionado

El algoritmo de Forward Selection identificó a las bandas B06 y B02 como las más relevantes para predecir los sólidos suspendidos en escala log-log. El modelo final presenta un buen desempeño, con un RMSE_test de 16.69 y un R²_test de 0.86 en escala original, lo que indica una capacidad predictiva adecuada sobre datos no vistos. Además, la ecuación obtenida permite expresar la variable objetivo en función de las bandas seleccionadas de forma clara y directa.

Aunque el ajuste es satisfactorio, sería conveniente en futuros análisis explorar modelos no lineales o técnicas multivariadas más complejas, como regresión regularizada o árboles de decisión, que podrían captar relaciones más sutiles entre variables.

Referencias

[1]

D. C. R. Ramírez, «Método de Estimación de Sólidos Suspendidos Totales como Indicador de la Calidad del Agua Mediante Imágenes Satelitales». Universidad Nacional de Colombia, Facultad de Ciencias Agrarias, 2017.

[2]

G. D. J. L., J. Stephan, y Delance-Martinic., «Determinación del parámetro sólidos suspendidos totales (sst) mediante imágenes de sensores ópticos en un tramo de la cuenca media del río Bogotá (Colombia).», UD y la Geomática, pp. 19-27, 2014, doi: 10.14483/udistrital.jour.udgeo.2014.9.a02.

[3]

A. Cruz-Retana, C. Fonseca-Ortiz, R. Becerril-Piña, M. A. Gómez-Albores, M. Hernández-Téllez, y R. Arévalo-Mejia, «Characterization of spectral reflectance and TSS concentration in continental water bodies worldwide», vol. 1. pp. 4-18, 2023.

[4]

E. E. C. Vargas, «Modelamiento de Relaciones entre Parámetros Fisicoquímicos y Microbiológicos en Aguas de la Bahía Interior del Lago Titicaca-Puno (Perú) mediante Árboles de Predicción», Revista Tecnica De La Facultad De Ingenieria Universidad Del Zulia, vol. 44, pp. 154-168, ago. 2021, doi: 10.22209/rt.v44n3a02.

[5]

M. Moeini, A. Shojaeizadeh, y M. Geza, «Supervised machine learning for estimation of total suspended solids in urban watersheds», Water (Switzerland), vol. 13, ene. 2021, doi: 10.3390/w13020147.

[6]

L. S. Kupssinskü, T. T. Guimarães, E. M. D. Souza, y D. C. Zanotta, «A method for chlorophyll-a and suspended solids prediction through remote sensing and machine learning», Sensors (Switzerland), vol. 20, abr. 2020, doi: 10.3390/s20072125.

Notas

Aguas lénticas.↩︎
d = prueba estadística de Durbin-Watson.↩︎