Tidyverse 5
Extraer Datos
No olvidar cargar los paquetes con
library().
En lugar de crear variables sueltas y luego unirlas, usamos mutate() para incrustar los cálculos directamente en nuestro objeto espacial.
La función extract() de terra devuelve una tabla, por lo que seleccionamos la segunda columna (que contiene los valores) para asignarla.
Variables Terreno
Preparamos las variables para extraer. Pendiente y aspecto, además de la elevación.
Preproceso de imágen Sentinel-2.
Cálculo NDVI.
Puntos clave de la lógica Tidy aplicada aquí:
- Encadenamiento: No creamos archivos temporales. El mapa entra, se le calculan cosas y sale enriquecido.
terra::extract: Es la función puente. Toma los raster (mapa_dtm,pendiente, etc.) y la geometría vectorial (ya sea polígono o punto) y cruza la información.
extraccion <- mapa_prop |>
select(rol) |>
mutate(
elevacion = terra::extract(mapa_dtm, #variable
mapa_prop, #límites
fun = mean, #función
na.rm = TRUE)[, 2]
) |>
mutate(
pendiente = terra::extract(pendiente,
mapa_prop,
fun = mean,
na.rm = TRUE)[, 2]
) |>
mutate(
aspecto = terra::extract(aspecto,
mapa_prop,
fun = mean,
na.rm = TRUE)[, 2]
) |>
mutate(
ndvi = terra::extract(ndvi,
mapa_prop,
fun = mean,
na.rm = TRUE)[, 2]
)Predicción de Aptitud Productiva Predial
Si en lugar de decir ¿Qué cobertura tengo?, nos preguntamos: Basado en la topografía y el vigor vegetal histórico (NDVI), ¿Cuál es el potencial productivo (o valor) de cada predio?
Esto implica un flujo de trabajo avanzado:
1- Ingeniería de Características (Feature Engineering): Generar variables derivadas complejas.
2- Extracción Zonal: Pasar del mundo Raster (píxeles) al Vector (predios).
3- Modelado (Regresión): Entrenar un algoritmo para predecir una variable continua (Rendimiento, Carbono, Precio).
1- Ingeniería de Características (Más allá de lo básico).
No te quedes solo con Pendiente y NDVI. Para un nivel avanzado, genera variables de contexto:
Topografía:
- TRI (Terrain Ruggedness Index): Rugosidad del terreno (crucial para maquinaria agrícola).
- TPI (Topographic Position Index): Para saber si el predio está en un valle (acumula agua) o en una cima (erosión).
Imagen:
- Textura (GLCM): La varianza o entropía del NDVI dentro del predio. Un predio con NDVI alto pero mucha varianza puede indicar manejo desordenado o malezas.
2- La Unidad de Análisis: El Predio (Objeto).
El error de novato es modelar píxeles y luego promediar. El enfoque avanzado es:
Resumir las variables ambientales por predio (Media, Desviación Estándar).
Entrenar el modelo usando el predio como fila en la tabla de datos.
# ------------------------------------------------------------------------------
# 4. SIMULACIÓN DE LA VARIABLE RESPUESTA (GROUND TRUTH)
# ------------------------------------------------------------------------------
# Supongamos que fuiste a terreno y mediste la "Productividad" (Ton/Ha).
# Aquí la inventamos con una fórmula lógica para que el modelo funcione:
# Más NDVI y menos pendiente = Más productividad.
datos_modelo <- extraccion |>
mutate(
PRODUCTIVIDAD = (ndvi * 10) - (pendiente * 0.1) + rnorm(n(), 0, 0.5)
)
# ------------------------------------------------------------------------------
# 5. ENTRENAMIENTO DEL MODELO (MACHINE LEARNING)
# ------------------------------------------------------------------------------
message("Entrenando Random Forest...")
# A. Limpieza: Quitamos la geometría para entrenar (RF no entiende polígonos)
tabla_entrenamiento <- datos_modelo |>
st_drop_geometry() |>
select(-rol) |> # Quitar IDs
na.omit()# B. Partición Train/Test (70/30)
set.seed(99)
index <- createDataPartition(tabla_entrenamiento$PRODUCTIVIDAD, p = 0.7, list = FALSE)
train_set <- tabla_entrenamiento[index, ]
test_set <- tabla_entrenamiento[-index, ]
# C. Modelo Random Forest (Regresión)
# Formula: Productividad depende de (~) todas las otras variables (.)
rf_model <- randomForest(PRODUCTIVIDAD ~ ., data = train_set, ntree = 500, importance = TRUE)
print(rf_model)Salida
Call:
randomForest(formula = PRODUCTIVIDAD ~ ., data = train_set, ntree = 500, importance = TRUE)
Type of random forest: regression. Number of trees: 500. No. of variables tried at each split: 1.
Mean of squared residuals: 0.3078295.
% Var explained: 78.93. # ------------------------------------------------------------------------------
# 6. EVALUACIÓN Y EXPLICABILIDAD
# ------------------------------------------------------------------------------
# A. Importancia de Variables (¿Qué define el éxito de un predio?)
varImpPlot(rf_model, main = "Importancia de Variables en la Productividad")
# B. Validación con Test Set
predicciones <- predict(rf_model, test_set)
error <- postResample(pred = predicciones, obs = test_set$PRODUCTIVIDAD)
print(error) # RMSE y Rsquared
Interpretación de los errores.
El MAE me dice cuánto me equivoco en el día a día, el RMSE me avisa si tengo errores peligrosamente grandes, y el R-cuadrado me dice qué tan bien entendí el fenómeno en general.
ver detalles en código de ‘RANDOM-FOREST.R’
# ------------------------------------------------------------------------------
# 7. PREDICCIÓN ESPACIAL (MAPA DE APTITUD)
# ------------------------------------------------------------------------------
# Imaginemos que tenemos predios nuevos sin datos de productividad.
# Usamos el modelo para estimarla.
# Predecimos sobre EL MISMO objeto sf (para efectos didácticos)
datos_modelo$PREDICCION_RF <- predict(rf_model, st_drop_geometry(datos_modelo))
DETECCIÓN DE ANOMALÍAS (RESIDUAL MAPPING)
Objetivo: Identificar predios que se comportan distinto a lo esperado por el modelo.
CÁLCULO DE RESIDUOS
Predecimos sobre TODOS los datos originales para ver el comportamiento global Importante: randomForest requiere quitar la geometría para predecir.
Entendiendo los resultados
Si dices “El error es de 500 kg”, depende del cultivo. Si es trigo, es mucho. Si es caña de azúcar, es nada.
Al usar Desviaciones Estándar (Z), estandarizas el lenguaje: Este predio está a 2 sigmas por debajo de lo normal. Eso es grave estadísticamente en cualquier contexto.
datos_gestion <- datos_modelo |>
mutate(
# A. Lo que el modelo dice que debería pasar (Teórico)
PREDICHO = predict(rf_model, st_drop_geometry(datos_modelo)),
# B. La diferencia con la realidad (Residuo Bruto)
# Residuo = Realidad - Predicción
RESIDUO = PRODUCTIVIDAD - PREDICHO,
# C. Estandarización (Z-Score) - CRÍTICO PARA DETECTAR ANOMALÍAS
# Nos dice cuántas desviaciones estándar se aleja de lo normal.
# Z > 1.5 es una anomalía fuerte positiva.
# Z < -1.5 es una anomalía fuerte negativa.
RESIDUO_STD = as.numeric(scale(RESIDUO))
)# ------------------------------------------------------------------------------
# 2. CLASIFICACIÓN DE GESTIÓN (SEMÁFORO)
# ------------------------------------------------------------------------------
datos_gestion <- datos_gestion |>
mutate(
ESTADO = case_when(
RESIDUO_STD < -1.5 ~ "CRÍTICO: Sub-rendimiento", # Produce mucho menos de lo esperado
RESIDUO_STD < -0.5 ~ "Alerta Baja",
RESIDUO_STD > 1.5 ~ "EXCEPCIONAL: Sobre-rendimiento", # Produce mucho más
RESIDUO_STD > 0.5 ~ "Buen Desempeño",
TRUE ~ "Normal (Según Modelo)" # Se comporta como el modelo predice
)
)# ------------------------------------------------------------------------------
# 3. VISUALIZACIÓN DIAGNÓSTICA (MAPA DE CALOR DE RESIDUOS)
# ------------------------------------------------------------------------------
# Este mapa muestra la magnitud del error.
# ROJO = Problema (Menos de lo esperado)
# AZUL = Éxito (Más de lo esperado)
# BLANCO = El modelo acertó
g_residuos <- ggplot(datos_gestion) +
geom_sf(aes(fill = RESIDUO), color = "black", lwd = 0.2) +
# Usamos scale_fill_gradient2 para mapas divergentes
scale_fill_gradient2(
low = "#D73027", # Rojo (Negativo)
mid = "white", # Cero (Acierto)
high = "#4575B4", # Azul (Positivo)
midpoint = 0,
name = "Diferencia\n(Ton/Ha)"
) +
theme_minimal() +
labs(title = "Mapa de Residuos", subtitle = "¿Dónde falló el modelo?")# ------------------------------------------------------------------------------
# 4. MAPA DE GESTIÓN (SEMÁFORO)
# ------------------------------------------------------------------------------
# Este es el mapa que se le entrega al gerente o agrónomo.
g_gestion <- ggplot(datos_gestion) +
geom_sf(aes(fill = ESTADO), color = "white", lwd = 0.1) +
scale_fill_manual(
values = c(
"CRÍTICO: Sub-rendimiento" = "#B2182B", # Rojo oscuro
"Alerta Baja" = "#EF8A62", # Naranja
"Normal (Según Modelo)" = "grey90", # Gris
"Buen Desempeño" = "#67A9CF", # Celeste
"EXCEPCIONAL: Sobre-rendimiento" = "#2166AC" # Azul oscuro
)
) +
theme_void() +
labs(
title = "Mapa de Alertas de Gestión",
subtitle = "Identificación de predios anómalos para visita a terreno"
)
# ------------------------------------------------------------------------------
# 5. REPORTE DE LOS PEORES CASOS (Para imprimir)
# ------------------------------------------------------------------------------
peores_casos <- datos_gestion |>
st_drop_geometry() |>
filter(RESIDUO_STD < -1.5) |>
select(rol, PRODUCTIVIDAD, PREDICHO, RESIDUO) |>
arrange(RESIDUO)
print("--- PREDIOS QUE REQUIEREN VISITA URGENTE ---")
head(peores_casos)