Hackaton

Reporte de Modelo de Predicción de Calidad del Agua

1. Introducción

1.1 Descripción del problema

El problema central consiste en predecir tres variables de calidad del agua en cuerpos de agua superficiales de Sudafrica a partir de datos de teledetección satelital y datos climáticos. Las tres variables objetivo son:

()imagenes

1.2 Fuentes de datos

Se utilizaron tres conjuntos de datos:

1. water_quality_training_dataset.csv: Contiene las mediciones de los tres parámetros de calidad. Incluye columnas de Latitude, Longitude y Sample Date como claves de cruce. Cada fila representa una muestra de agua tomada en un punto geográfico específico en una fecha determinada.

2. landsat_features_training.csv: Contiene los valores de reflectancia de superficie de las bandas Landsat 8/9, extraídas en los mismos puntos y fechas que las muestras de campo. Las bandas disponibles fueron:

   • nir: Infrarrojo cercano (~865 nm), altamente sensible a la biomasa y turbidez [^1]
   • green: Verde (~560 nm), sensible a la clorofila [^2]
   • swir16: Infrarrojo de onda corta 1 (~1610 nm), sensible a la humedad del suelo y turbidez
   • swir22: Infrarrojo de onda corta 2 (~2200 nm), partículas suspendidas en agua
   • NDMI: Índice de Diferencia Normalizada de Humedad, calculado como (NIR - SWIR16) / (NIR + SWIR16)
   • MNDWI: Índice de Agua de Diferencia Normalizada Modificado, calculado como (Green - SWIR16) / (Green + SWIR16)

3. terraclimate_features_training.csv: Contiene variables climáticas del modelo TerraClimate, incluyendo:

   • pet: Evapotranspiración potencial (mm/mes), indicador de condiciones hidrológicas y aridez local [^3]

1.3 Punto de partida: el modelo Benchmark

El notebook de referencia (Benchmark_Model_Notebook.ipynb) proporcionado por los organizadores de la competencia implementaba lo siguiente:

• Carga y cruce de los tres datasets por coordenadas y fecha, luego llena los valores faltantes con la media
• Sin escalado de características
• Modelo único por variable objetivo: HistGradientBoostingRegressor con parámetros por defecto
• Sin ingeniería de características adicional más allá de las bandas espectrales y el pet

Resultado del benchmark original: R² promedio ≈ 0.52.

# Suppress warnings
import warnings
warnings.filterwarnings('ignore')

# Visualization libraries
import matplotlib.pyplot as plt
import seaborn as sns

# Data manipulation and analysis
import numpy as np
import pandas as pd

# Multi-dimensional arrays and datasets (e.g., NetCDF, Zarr)
import xarray as xr

# Geospatial raster data handling with CRS support
import rioxarray as rxr

# Raster operations and spatial windowing
import rasterio
from rasterio.windows import Window

# Feature preprocessing and data splitting
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from scipy.spatial import cKDTree

# Machine Learning
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error

# Planetary Computer tools for STAC API access and authentication
import pystac_client
import planetary_computer as pc
from odc.stac import stac_load
from pystac.extensions.eo import EOExtension as eo

from datetime import date
from tqdm import tqdm
import os 

Water_Quality_df=pd.read_csv('water_quality_training_dataset.csv')
landsat_train_features = pd.read_csv('landsat_features_training.csv')
Terraclimate_df = pd.read_csv('terraclimate_features_training.csv')

wq_data = combine_two_datasets(Water_Quality_df, landsat_train_features, Terraclimate_df)
wq_data = wq_data.fillna(wq_data.median(numeric_only=True))

wq_data = wq_data[['swir22','NDMI','MNDWI','pet', 'Total Alkalinity', 'Electrical Conductance', 'Dissolved Reactive Phosphorus']]

def split_data(X, y, test_size=0.3, random_state=42):
    return train_test_split(X, y, test_size=test_size, random_state=random_state)

def scale_data(X_train, X_test):
    scaler = StandardScaler()
    X_train_scaled = scaler.fit_transform(X_train)
    X_test_scaled = scaler.transform(X_test)
    return X_train_scaled, X_test_scaled, scaler

def train_model(X_train_scaled, y_train):
    model = RandomForestRegressor(n_estimators=100, random_state=42)
    model.fit(X_train_scaled, y_train)
    return model

def evaluate_model(model, X_scaled, y_true, dataset_name="Test"):
    y_pred = model.predict(X_scaled)
    r2 = r2_score(y_true, y_pred)
    rmse = np.sqrt(mean_squared_error(y_true, y_pred))
    print(f"\n{dataset_name} Evaluation:")
    print(f"R²: {r2:.3f}")
    print(f"RMSE: {rmse:.3f}")
    return y_pred, r2, rmse

def run_pipeline(X, y, param_name="Parameter"):
    print(f"\n{'='*60}")
    print(f"Training Model for {param_name}")
    print(f"{'='*60}")
    
    X_train, X_test, y_train, y_test = split_data(X, y)
    
    X_train_scaled, X_test_scaled, scaler = scale_data(X_train, X_test)
    
    model = train_model(X_train_scaled, y_train)
    
    y_train_pred, r2_train, rmse_train = evaluate_model(model, X_train_scaled, y_train, "Train")
    
    y_test_pred, r2_test, rmse_test = evaluate_model(model, X_test_scaled, y_test, "Test")
    
    results = {
        "Parameter": param_name,
        "R2_Train": r2_train,
        "RMSE_Train": rmse_train,
        "R2_Test": r2_test,
        "RMSE_Test": rmse_test
    }
    return model, scaler, pd.DataFrame([results])

  X = wq_data.drop(columns=['Total Alkalinity', 'Electrical Conductance', 'Dissolved Reactive Phosphorus'])

y_TA = wq_data['Total Alkalinity']
y_EC = wq_data['Electrical Conductance']
y_DRP = wq_data['Dissolved Reactive Phosphorus']

model_TA, scaler_TA, results_TA = run_pipeline(X, y_TA, "Total Alkalinity")
model_EC, scaler_EC, results_EC = run_pipeline(X, y_EC, "Electrical Conductance")
model_DRP, scaler_DRP, results_DRP = run_pipeline(X, y_DRP, "Dissolved Reactive Phosphorus")

Con los resultados:

============================================================
Training Model for Total Alkalinity
============================================================

Train Evaluation:
R²: 0.903
RMSE: 23.132

Test Evaluation:
R²: 0.546
RMSE: 50.870

============================================================
Training Model for Electrical Conductance
============================================================

Train Evaluation:
R²: 0.918
RMSE: 98.007

Test Evaluation:
R²: 0.585
RMSE: 219.999

============================================================
Training Model for Dissolved Reactive Phosphorus
============================================================

Train Evaluation:
R²: 0.882
RMSE: 17.455

Test Evaluation:
R²: 0.529
RMSE: 35.182

2. Preprocesamiento de Datos

2.1 Carga y fusión de datasets

Los tres archivos CSV se cargan y fusionan mediante un inner join usando como claves las columnas Latitude, Longitude y Sample Date. Este tipo de unión garantiza que las muestras se unan en los mismos datos

merged_df = pd.merge(wq_df, landsat_df, on=['Latitude', 'Longitude', 'Sample Date'], how='inner')
merged_df = pd.merge(merged_df, terraclimate_df, on=['Latitude', 'Longitude', 'Sample Date'], how='inner')

2.2 Eliminación de filas con targets faltantes

Como parte del preprocesamiento se eliminan las filas donde las variables objetivo/dependientes tienen valores nulos:

targets = ['Total Alkalinity', 'Electrical Conductance', 'Dissolved Reactive Phosphorus']
merged_df = merged_df.dropna(subset=targets)

No tiene sentido entrenar sobre ejemplos donde el valor a predecir es desconocido. Pues puede crear predicciones erroneas en el modelo

2.3 Imputación de features

El modelo original del benchmark usaba imputación con la media aritmética.

val_data = val_data.fillna(val_data.median(numeric_only=True))

En lugar se eso se decidio escoger la mediana, esto por que la media es muy sensible a los outliers [^4]

feature_cols = ['nir', 'green', 'swir16', 'swir22', 'NDMI', 'MNDWI', 'pet']
for col in feature_cols:
    merged_df[col] = merged_df[col].fillna(merged_df[col].median())

Esto significa que la mediana no cambia significativamente aunque se agreguen valores extremadamente grandes o pequeños

2.4 Escalado de características: StandardScaler

Se agregó un escalado en este entrenamiento.

scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_val_scaled = scaler.transform(X_test)  # transform únicamente, sin fit

El StandardScaler transforma cada feature: X_scaled = (X - μ) / σ, donde μ es la media y σ la desviación estándar del conjunto de entrenamiento. Esto es necesario por que laa latitud varía en ~10 grados, las bandas nir/green varían entre 0 y 0.5, y los ratios espectrales tienen valores superiores a 200. Los algoritmos con regularización le dan más peso a las variables de mayor escala numérica y en el StackingRegressor, el Ridge recibe como input las predicciones de 5 modelos base. Si esas predicciones tienen escalas muy diferentes (porque los targets tienen órdenes de magnitud distintos), la Ridge tiene dificultades para encontrar una combinación lineal óptima.

El scaler se ajusta únicamente sobre el conjunto de entrenamiento (fit_transform) y se aplica (sin reajustar) sobre validación y test (transform). Si el scaler viera los datos de validación durante fit, la media y σ que calcularía incorporarían información “del futuro”, haciendo que los resultados de validación sean artificialmente optimistas y no representativos del rendimiento real en producción. El scaler se guarda como scaler.pkl para garantizar la misma transformación en inferencia.

3. Ingeniería de Características (Feature Engineering)

La ingeniería de características fue el componente con mayor impacto individual en el rendimiento del modelo. Partiendo solo de las 7 features del benchmark (4 bandas Landsat + NDMI + MNDWI + pet), se llegó a un total de 30 features.

3.1 Características temporales

3.1.1 Mes y día del año

merged_df['Month'] = merged_df['Sample Date'].dt.month
merged_df['DayOfYear'] = merged_df['Sample Date'].dt.dayofyear

3.1.2 Codificación cíclica: sin_doy y cos_doy

merged_df['sin_doy'] = np.sin(2 * np.pi * merged_df['DayOfYear'] / 365)
merged_df['cos_doy'] = np.cos(2 * np.pi * merged_df['DayOfYear'] / 365)

¿Por qué codificación cíclica? Si se usa DayOfYear directamente como feature numérica, el modelo aprende que el día 1 (1 de enero) y el día 365 (31 de diciembre) están a una “distancia” de 364 días, cuando en realidad están a un solo día de diferencia (el año es cíclico). Los árboles de decisión no tienen este problema por naturaleza, pero la codificación cíclica facilita que el modelo encuetre patrones periódicos más fácilmente. [^5]

Las funciones sin y cos del día del año crean un espacio bidimensional donde puntos cercanos en el tiempo circular están también cercanos en el espacio de features: el 31 de diciembre y el 1 de enero quedan a distancia pequeña en (sin_doy, cos_doy).

3.2 Características espaciales

3.2.1 Clustering espacial con K-Means

from sklearn.cluster import KMeans
coords = merged_df[['Latitude', 'Longitude']]
kmeans = KMeans(n_clusters=10, random_state=42, n_init=10)
merged_df['Cluster'] = kmeans.fit_predict(coords)

El K-Means con 10 clusters divide el territorio geográfico en 10 regiones. La etiqueta del cluster (Cluster = 0 a 9) permite al modelo aprender efectos regionales sin tener que descubrir las fronteras por sí solo a partir de la latitud y longitud continuas.

Se escogieron 10 clusters por que considero que tiene la granularidad suficiente para capturar variación regional sin introducir overfitting

Cuando se introdujo la caracteristica de Cluster el modelo mejoró el DRP de 0.6964 a 0.7077.

3.2.3 Interacciones espaciales-espectrales

merged_df['nir_lat'] = merged_df['nir'] * merged_df['Latitude']
merged_df['green_lat'] = merged_df['green'] * merged_df['Latitude']
merged_df['swir16_lat'] = merged_df['swir16'] * merged_df['Latitude']

la interacción banda × latitud permite al modelo capturar variaciones espaciales en la relación entre reflectancia y la variable objetivo.

Estos términos permiten modelar efectos donde la influencia de una variable depende de otra, lo cual no puede representarse mediante modelos puramente aditivos. [^6]

3.3 Transformaciones logarítmicas

for col in ['green', 'nir', 'swir16']:
    merged_df[f'log_{col}'] = np.log1p(merged_df[col])

Las distribuciones de las bandas espectrales son asimétricas. La transformación log1p(x) = log(1 + x) comprime los valores extremos y expande los valores pequeños, haciendo la distribución más simétrica.

Se usa log1p en lugar de log para manejar valores igual a cero sin producir -inf.

3.4 Ratios espectrales

merged_df['nir_green_ratio']       = merged_df['nir']    / (merged_df['green']  + 1e-6)
merged_df['swir16_swir22_ratio']   = merged_df['swir16'] / (merged_df['swir22'] + 1e-6)
merged_df['green_swir16_ratio']    = merged_df['green']  / (merged_df['swir16'] + 1e-6)
merged_df['green_nir_ratio']       = merged_df['green']  / (merged_df['nir']    + 1e-6)
merged_df['swir16_nir_ratio']      = merged_df['swir16'] / (merged_df['nir']    + 1e-6)

Los ratios de bandas son indicadores espectrales ampliamente usados en teledetección de calidad del agua. Son invariantes ante cambios de iluminación (ángulo solar, condiciones atmosféricas) porque afectan a numerador y denominador por igual, cancelándose. Se usa +1e-6 en el denominador para evitar divisiones por cero sin alterar significativamente el resultado.

3.5 Features polinómicas y de interacción

# Cuadráticas (captura relaciones no lineales)
merged_df['nir_sq']    = merged_df['nir']    ** 2
merged_df['green_sq']  = merged_df['green']  ** 2
merged_df['swir16_sq'] = merged_df['swir16'] ** 2

# Interacciones (captura efectos conjuntos entre bandas)
merged_df['nir_green_inter']  = merged_df['nir']  * merged_df['green']
merged_df['nir_swir16_inter'] = merged_df['nir']  * merged_df['swir16']

Los modelos de árbol pueden aproximar términos polinómicos a través de splits sucesivos, pero agregar explícitamente las features cuadráticas les permite descubrir esas relaciones con menos splits y de forma más eficiente (menor profundidad requerida, menos riesgo de overfitting).

Las interacciones capturan efectos sinérgicos: la relación entre nir y green juntos puede predecir la clorofila mejor que cada banda por separado, porque la clorofila absorbe fuertemente en rojo pero no en verde ni en NIR. [^11]

3.6 Resumen del espacio de features final

4. Evolución de la Arquitectura del Modelo

4.1 Estado inicial

El benchmark usaba un único HistGradientBoostingRegressor (HGBR) por variable objetivo, con los parámetros por defecto de scikit-learn. Este es un modelo de gradient boosting basado en histogramas (similar a LightGBM), que es eficiente y robusto, pero con hiperparámetros sin optimizar y sin features adicionales su rendimiento estaba limitado.

Resultado inicial:
  Total Alkalinity:             R² ≈ 0.52
  Electrical Conductance:       R² ≈ 0.55
  Dissolved Reactive Phosphorus: R² ≈ 0.49
  Promedio:                     R² ≈ 0.52

4.2 Primera iteración: VotingRegressor (HGBR + RF) → R² ≈ 0.79

El primer paso fue crear un ensemble con VotingRegressor combinando HGBR y RandomForestRegressor:

VotingRegressor([
    ('hgb', HistGradientBoostingRegressor(...)),
    ('rf',  RandomForestRegressor(...))
])

El VotingRegressor promedia las predicciones de múltiples modelos. Cuando dos modelos cometen errores en diferentes observaciones (sus errores no están perfectamente correlacionados), el promedio de sus predicciones tiene un error total menor que cualquiera de los dos por separado. Este efecto se llama reducción de varianza por ensemble.

Esto resulto en un R² ≈ 0.79. Un salto de 27 puntos porcentuales en referencia al benchmark, principalmente gracias a la combinación de los patrones que cada modelo captura de forma diferente.

4.3 Segunda iteración: Agregado de ExtraTreesRegressor → R² ≈ 0.80

Se agrega un tercer estimador al VotingRegressor: ExtraTreesRegressor (Extremely Randomized Trees):

VotingRegressor([
    ('hgb', HistGradientBoostingRegressor(...)),
    ('rf',  RandomForestRegressor(...)),
    ('et',  ExtraTreesRegressor(...))
])

¿Qué aporta ExtraTrees? A diferencia del RandomForest (que selecciona el mejor split entre un subconjunto aleatorio de features), ExtraTrees selecciona tanto las features como los umbrales de split de forma completamente aleatoria. Esto produce modelos individuales con más varianza, pero el promedio de muchos de ellos converge a una función de menor sesgo. Para datasets ruidosos como el DRP, la aleatoriedad extra puede ser beneficiosa.

Impacto: R² ≈ 0.80 (mejora marginal de +0.01), pero establece la base para el Stacking posterior.

4.4 Tercer iteración: StackingRegressor → R² ≈ 0.8123

Se reemplaza el VotingRegressor por un StackingRegressor con Ridge como meta-estimador:

StackingRegressor(
    estimators=[('hgb', best_hgb), ('rf', best_rf), ('et', best_et)],
    final_estimator=Ridge(alpha=1.0),
    cv=5,
    n_jobs=-1
)

¿Por qué Stacking supera a Voting?

El VotingRegressor asigna pesos iguales a todos los modelos. Si un modelo es significativamente mejor que otros para un subconjunto de observaciones (por ejemplo, RF es mejor para datos espacialmente concentrados, mientras HGBR es mejor para datos temporales), el promedio simple no aprovecha esta especialización.[^12]

El StackingRegressor entrena un meta-estimador que aprende a combinar las predicciones de los modelos. [^13]

El StackingRegresor funciona primero entrenando los 3 modelos base que se entrenan con validación cruzada de 5 pliegues (cv=5), que consiste en que para cada fold, se generan predicciones sobre datos en los que el modelo no se entreno, estas predicciones se usan como features para entrenar la Ridge, este modelo aprende los pesos óptimos para combinar los 3 modelos, con regularización para evitar overfitting

Resultados con Stacking (3 modelos):
  Total Alkalinity:             R² = 0.8490
  Electrical Conductance:       R² = 0.8702
  Dissolved Reactive Phosphorus: R² = 0.7178
  Promedio:                     R² = 0.8123

Se uso Ridge como meta-estimador pues este cuenta con regularización L2, que evita que el meta-estimador asigne pesos extremos a un solo modelo base

4.5 Cuarta iteración: XGBoost + LightGBM + Stacking de 5 modelos (estado actual)

Se integran XGBoost y LightGBM como modelos base adicionales:

StackingRegressor(
    estimators=[
        ('xgb',  best_xgb),
        ('lgbm', best_lgbm),
        ('hgb',  best_hgb),
        ('rf',   best_rf),
        ('et',   best_et),
    ],
    final_estimator=Ridge(alpha=1.0),
    cv=5
)

Se agregaron los modelos de XGBoost y LightGBM para sustituir al HGBR

  Total Alkalinity:             R² = 0.8540
  Electrical Conductance:       R² = 0.8702
  Dissolved Reactive Phosphorus: R² = 0.7278
  Promedio:                     R² = 0.8173

4.6 Quinta iteración: Nuevas Features

Despues de la iteracion anterior, note que existe un limite en los resultados, específicamente en DRP, asi que descidí aumentar las features del sistema a las siguientes:

2.1.1 Features Espectrales

NDWI = (green - nir) / (green + nir + ε)
Turbidity = green / (swir16 + ε)
NIR_GREEN_ratio = nir / (green + ε)
SWIR16_NIR_ratio = swir16 / (nir + ε)
SWIR22_SWIR16_ratio = swir22 / (swir16 + ε)
NDTI = (green - swir16) / (green + swir16 + ε)
WRI = (green + nir) / (swir16 + swir22 + ε)
AWEI = 4(green - swir16) - (0.25·nir + 2.75·swir22)
Brightness = √(nir² + green² + swir16² + swir22²)
Spectral_Variance = var(nir, green, swir16, swir22)

• NDWI: Discrimina agua de vegetación/suelo
• Turbidity: Proxy de sedimentos suspendidos (correlaciona con alcalinidad)
• Ratios espectrales: Capturan relaciones no lineales
• AWEI: Índice de agua mejorado para aguas turbias
• Brightness & Variance: Capturan propiedades ópticas generales

2.1.2 Features Temporales

Month, Year, DayOfYear
sin_doy = sin(2π · DayOfYear / 365)
cos_doy = cos(2π · DayOfYear / 365)
Season = ((Month % 12 + 3) // 3)
Is_Summer, Is_Winter

• Ciclicidad: sin/cos capturan estacionalidad circular
• Season: Algunos procesos químicos varían por estación
• Month/Year: Tendencias temporales

2.1.3 Features Geoespaciales

Cluster = KMeans(n_clusters=10)
Dist_to_Center = √((lat - lat_mean)² + (lon - lon_mean)²)
Dist_to_Cluster_Center

2.1.4 Features de Interacción

NDMI_MNDWI_product
NDWI_squared, NDWI_cubed
Turbidity_log
PET_Month, PET_sin_doy
Month_squared, Year_scaled

• Productos: Capturan interacciones multiplicativas
• Polinomios: Relaciones no lineales
• Log transforms: Manejan distribuciones asimétricas

2.1.5 Por Clusters

Cluster_mean_nir, nir_deviation
Cluster_mean_green, green_deviation
Cluster_mean_swir16, swir16_deviation
Cluster_mean_NDWI, NDWI_deviation
Cluster_mean_Turbidity, Turbidity_deviation

• Cluster means: Contexto local
• Deviations: Anomalías respecto al cluster

Iteración 6 — Bagging Agresivo y Simplificación de Arquitectura → R² = 0.8331

Después de experimentar con los hiperparametros de Stacking y observar que solo empeoraba el R2 y que varios modelos (RF, ET, HGBR) daban predicciones muy similares entre sí, se simplifico la arquitectura enfocándose solo en los dos algoritmos de mejor precision.

Analizando como se comportaba el modelo, note que los modelos de RandomForest, ExtraTrees y HistGradientBoostingRegressor no tenian gran efecto en el modelo final, por lo que se eliminaron

En este punto se decidio escoger una tecnica de bagging[^14], la cual consisten en entrenar diferentes modelos en diferentes subsets de los datos y usarlo para hacer una prediccion.

B) Bagging de XGBoost y LightGBM:

# Ahora: 10 XGBoost + 10 LightGBM con diferentes seeds
xgb_models = [
    XGBRegressor(**best_params, random_state=42+i) 
    for i in range(10)
]
lgb_models = [
    LGBMRegressor(**best_params, random_state=42+i) 
    for i in range(10)
]

# Predicción final = promedio de las 20 predicciones
y_pred = np.mean([model.predict(X) for model in xgb_models + lgb_models], axis=0)

C) Expansión masiva del espacio de features:

• De 30 features → 55 features
• Se agregaron features espectrales avanzadas (NDWI, Turbidity, NDTI, WRI, AWEI, Brightness, Spectral_Variance)
• Features temporales estacionales (Season, Is_Summer, Is_Winter)
• Features por cluster (cluster means y deviations para múltiples variables)
• Clustering espacial aumentado de k=10 → k=100 para capturar micro-regiones
• Features de distancia a hotspots/coldspots

D) Hiperparámetros específicos por target: Para DRP (el target más difícil), se usaron configuraciones mucho más agresivas:

# DRP: Alta complejidad + regularización fuerte
n_estimators=2800, learning_rate=0.007, max_depth=15
reg_lambda=3.0  # regularización L2 muy fuerte

# TA/EC: Menos complejidad
n_estimators=2200, learning_rate=0.009, max_depth=12-13
reg_lambda=1.8-2.2

Resultados:

Con esto se mostro que el bagging superó al stacking con modelos diversos,

La estrategia de entrenar 10 y 10 modelos produjo mejores resultados que todos juntos. y aque cada modelo con su entrenamiento diferente, genera predicciones independientes cuyo promedio tiene menor varianza.

El análisis de correlaciones reveló que la feature con mayor correlación con DRP era solo r=0.365 (Cluster_Density), mientras que para Conductance era r=0.687 (Turbidity). Esta diferencia daba a entender que el problema no era del de modelo, sino de la información disponible,

Se decidió explorar otro datasets que pudieran aportar features.

Se muestra que DRP reducia el promedio del R2, asi que se analizo mejor estos datos:

Mean:    51.32 μg/L
Std:     51.31 μg/L  (casi igual a la media = alta varianza)
Min:     2.10 μg/L
Q1:      17.50 μg/L
Median:  34.00 μg/L
Q3:      67.00 μg/L
Max:     485.00 μg/L  (outlier extremo, 9.5 x mediana)
Skewness: 2.84  (muy sesgada a la derecha)

[Mostrar grafica]

6.2 El problema del DRP

El Dissolved Reactive Phosphorus se mantuvo consistentemente como el objetivo con la menor prescicion, el más difícil de predecir, con una R² entre 0.70 y 0.72 en los mejores modelos. El DRP tiene una distribucion en la que la mayoría de los valores están concentrados en el rango 0–50 µg/L, pero existen outliers de hasta 500+ µg/L. El R² penaliza fuertemente los errores cuadráticos en estos outliers, haciendo que sea difícil mejorar la métrica sin capturar bien los extremos.

6.2.2 Relación débil con las bandas espectrales

El DRP es un ion disuelto en el agua. A diferencia de la turbidez o la clorofila (que tienen señales espectrales directas y bien documentadas), el ortofosfato disuelto no absorbe ni refleja en las longitudes de onda del satélite Landsat[^15] de forma significativa. La relación entre los valores de reflectancia y el DRP es indirecta [^15]

El análisis de importancias (ExtraTrees sobre conjunto de entrenamiento) reveló:

La dominancia de las features espaciales (Latitude, Longitude, Cluster con importancia combinada = 0.56) sobre las espectrales sugiere que el modelo está principalmente aprendiendo dónde hay alta concentración de DRP en el territorio, más que el valor espectral que corresponde a ese nivel de fósforo.

8. Conclusiones y Próximos Pasos

8.1 Resumen de avances logrados

Se logró un incremento de +29 puntos porcentuales en R² promedio respecto al benchmark.

Evaluando los resultados

Con estos modelos, empece a enviar los resultados en los valores de validacion que proporciona el concurso, con esto me encontre que estos modelos tuvieron un ~0.1 en la evaluacion del concurso a pesar de tener ~0.8 de R2 en la validacion local.

6. Fase de Robustez y Transformación de Datos

Para esto se usaron técnicas para manejar valores atípicos (outliers). En estos cambios se uso RobustScaler en lugar de StandardScaler, el cual escala usando la mediana y el rango intercuartílico (IQR) en lugar de la media y desviación estándar, haciéndolo inmune a outliers extremos y el clipping de valores de reflectancia negativos o físicamente imposibles (valores fuera del rango [0, 1] para bandas normalizadas).

from sklearn.preprocessing import RobustScaler

# Corregir reflectancias físicamente imposibles
for col in ['nir', 'green', 'swir16', 'swir22']:
    merged_df[col] = merged_df[col].clip(lower=0.0, upper=1.0)

scaler = RobustScaler()
X_train_scaled = scaler.fit_transform(X_train)

Los datos de calidad del agua suelen tener distribuciones fuertemente sesgadas. Predecir el valor crudo hace que el modelo concentre su capacidad en los valores extremos y pierda precisión en el rango más frecuente.

make_log_target_sub.py

Para este paso se aplicó una transformación logarítmica directamente sobre el objetivo (np.log1p). El modelo se entrena para predecir log(1 + y) en lugar de y directamente, y al momento de generar las predicciones finales se aplica la transformación inversa con np.expm1. Adicionalmente se combinaron XGBoost, LightGBM y GradientBoosting en un ensamble, y se añadieron características de Planetary Computer (lulc, dem_slope) junto con la codificación cíclica temporal (sin_doy).

# Transformar target
y_train_log = np.log1p(y_train)
model.fit(X_train, y_train_log)

# Predicción con transformación inversa
y_pred = np.expm1(model.predict(X_test))

La transformación logarítmica convierte la distribución sesgada a la derecha en una aproximadamente normal, lo que permite que el modelo optimice el MSE de forma simétrica en lugar de concentrar gradientes en los outliers extremos. Aunque esta transformación estabilizó la varianza, el modelo aún podía mejorar mediante una validación más rigurosa.

7. Fase de Optimización y Validación Cruzada

make_kfold_sub.py

En lugar de una sola partición train/test, se adoptó una estrategia de K-Fold Cross Validation con 10 pliegues. El modelo se entrena 10 veces con diferentes subconjuntos de datos y las predicciones sobre el set de prueba se promedian entre todos los pliegues.

from sklearn.model_selection import KFold

kf = KFold(n_splits=10, shuffle=True, random_state=42)
test_preds = np.zeros(len(X_test))

for fold, (train_idx, val_idx) in enumerate(kf.split(X_train)):
    X_tr, X_val = X_train[train_idx], X_train[val_idx]
    y_tr, y_val = y_train[train_idx], y_train[val_idx]

    model.fit(X_tr, y_tr)
    test_preds += model.predict(X_test) / kf.n_splits

El promedio de las 10 predicciones sobre test reduce drásticamente la varianza del estimador final, ya que cada fold expone al modelo a una distribución diferente de los datos, reduciendo el riesgo de que el resultado dependa de una partición particular.

make_huber_sub.py

Con la validación cruzada establecida, el siguiente paso fue explorar funciones de pérdida alternativas para reducir la sensibilidad a errores extremos. Se adoptó la Huber Loss, un híbrido entre MSE y MAE: para errores pequeños se comporta como MSE (derivada suave, convergencia rápida) y para errores grandes se comporta como MAE (gradiente acotado, penaliza menos los outliers). El parámetro delta controla la frontera entre ambos regímenes.

# XGBoost con Huber loss
xgb_model = XGBRegressor(
    objective='reg:pseudohubererror',
    huber_slope=1.0,
)

# LightGBM con Huber loss
lgb_model = LGBMRegressor(
    objective='huber',
    alpha=0.9,
)

8. Fase de Ensembles Avanzados

make_catboost_ensemble.py

Para aumentar la diversidad del ensamble se incorporaron CatBoost y HistGradientBoosting, llegando a un total de 5 modelos base: XGBoost, LightGBM, GradientBoosting, CatBoost e HistGradientBoosting. Paralelamente se refinó el espacio de features, seleccionando las 12 más robustas mediante la importancia promediada sobre los 5 modelos y eliminando las que añadían ruido sin mejorar la generalización.

models = [
    ('xgb',  XGBRegressor(**xgb_params)),
    ('lgbm', LGBMRegressor(**lgbm_params)),
    ('gbm',  GradientBoostingRegressor(**gbm_params)),
    ('cat',  CatBoostRegressor(**cat_params, verbose=0)),
    ('hgb',  HistGradientBoostingRegressor(**hgb_params)),
]

test_preds = np.mean([m.predict(X_test) for _, m in models], axis=0)

CatBoost aporta diversidad real al ensamble porque su implementación interna del boosting (ordered boosting) difiere algorítmicamente de XGBoost y LightGBM, produciendo errores poco correlacionados con los de los otros modelos, lo cual es justamente lo que hace útil el promedio.

9. Fase Final: Técnicas Post-Entrenamiento

make_pseudolabel_sub.py

Con un ensamble sólido y estable, la única vía de mejora sustancial era incorporar información del propio set de prueba. Se implementó pseudo-etiquetado (semi-supervised learning): las predicciones más confiables de la mejor sumisión anterior se tratan como si fueran datos reales y se usan para re-entrenar el modelo junto con los datos originales.

# Filtrar predicciones de alta confianza (menor desacuerdo entre modelos)
pseudo_mask = ensemble_std < confidence_threshold
X_pseudo = X_test[pseudo_mask]
y_pseudo = best_submission_preds[pseudo_mask]

# Re-entrenar con datos originales + pseudo-etiquetas
X_augmented = np.vstack([X_train, X_pseudo])
y_augmented = np.concatenate([y_train, y_pseudo])
model.fit(X_augmented, y_augmented)

Adicionalmente se aplicaron restricciones físicas de dominio sobre las predicciones finales para descartar valores hidroquímicamente imposibles, como que la conductividad eléctrica no puede exceder 1.5 veces la alcalinidad total, o que las concentraciones no pueden superar límites físicos observados en la región.

predictions['EC'] = np.where(
    predictions['EC'] > 1.5 * predictions['TA'],
    1.5 * predictions['TA'],
    predictions['EC']
)
predictions['TA']  = predictions['TA'].clip(upper=1200)
predictions['EC']  = predictions['EC'].clip(upper=5000)
predictions['DRP'] = predictions['DRP'].clip(upper=600)

10. Análisis de Progresión de Envíos

Para evaluar objetivamente si las iteraciones estaban convergiendo o simplemente oscilando, se desarrolló un programa que carga el historial completo de envíos a la competencia y grafica la evolución del R² público a lo largo del tiempo. Esto permitió identificar visualmente qué fases generaron saltos reales en el score público, cuáles producían mejoras en validación local que no se transferían al set de prueba, y el plateau final donde los envíos comenzaron a oscilar alrededor de un valor máximo sin mejora sostenida, señal de que se había alcanzado el límite de la información disponible en las features espectrales para predecir DRP.

fig, ax = plt.subplots(figsize=(14, 5))

ax.plot(submissions['timestamp'], submissions['r2_public'],
        marker='o', linewidth=1.5, color='steelblue', zorder=2)

for phase, color in phase_colors.items():
    mask = submissions['phase'] == phase
    ax.scatter(submissions.loc[mask, 'timestamp'],
               submissions.loc[mask, 'r2_public'],
               color=color, s=60, label=phase, zorder=3)

ax.axhline(y=0.52, linestyle='--', color='red', alpha=0.6, label='Benchmark (0.52)')
ax.set_xlabel('Fecha del envío')
ax.set_ylabel('R² público')
ax.set_title('Evolución del score por envío')
ax.legend()
plt.tight_layout()

[^1] https://en.wikipedia.org/wiki/Near-infrared_spectroscopy [^2] https://www.gisandbeers.com/caracteristicas-las-imagenes-satelite-landsat-9/ [^3] https://en.wikipedia.org/wiki/Potential_evapotranspiration [^4] https://en.wikipedia.org/wiki/Robust_statistics [^5] https://deepwiki.com/WillKoehrsen/Data-Analysis/3.2-cyclical-feature-encoding [^6] https://feat.engineering/01-introduction.html [^7] https://science.nasa.gov/earth/earth-observatory/measuring-vegetation-ndvi-evi/ [^9] https://swiftgeospatial.solutions/2025/03/12/the-basics-of-spectral-bands-nir-swir-and-rgb/#toc_Shortwave_Infrared_SWIR [^10] https://www.sciencedirect.com/science/article/abs/pii/S0034425796000673 [^11] https://www.sciencedirect.com/science/article/abs/pii/S0176161704704034 [^12] https://www.geeksforgeeks.org/machine-learning/voting-regressor/ [^13] https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.StackingRegressor.html [^14] https://en.wikipedia.org/wiki/Bootstrap_aggregating [^15] https://www.mdpi.com/2073-4441/13/12/1704