Taller 5 Región Metropolitana

Magister en Teledetección. Análisis de Datos Espaciales con R. Interpolación espacial con Random Forest

Introducción:

Durante los cuatro talleres anteriores se trabajó el método geoestadístico de la manera tradicional, primero se descargó, proceso y ordeno la información, luego se generaron predictores con los cuales se realizaron regresiones lineales, que vendría siendo la parte determinística del proceso. Posteriormente, el trabajo se concentró en la autocorrelación espacial mediante los cálculos de los índices de Moran global, local y la función del variograma. Con esta función se ajustaron distintos modelos para cada uno de los meses considerados. Para finalmente, utilizar como insumos estos variogramas modelados para construir distintos tipos de Kriging que según sus características propias asignaron distintos pesos a los datos y construyeron modelos de interpolación capaces de cubrir la Región Metropolitana en su totalidad con una muy buena precisión.

El presente taller abordará el tema de la interpolación de la temperatura en la Región Metropolitana desde una perspectiva diferente. Utilizando en esta oportunidad una técnica de Machine Learning, específicamente el método denominado Random Forest el cual se utiliza principalmente para tareas de clasificación y regresión. Este es un método de Machine Learning supervisado que opera construyendo múltiples árboles de decisiones. Estos árboles entregan diferentes predicciones y el algoritmo promedia el valor de cada uno de los árboles para entregar una respuesta final. Por esta razón, el método cobra fuerza al tener un gran número de datos y varios predictores.

La Geoestadística tradicional tiene su origen en la década de los 50, pero no fue hasta los años 60 cuando el matemático Georges Matheron desarrollo nuevos métodos para modelar y predecir la variabilidad espacial de los fenómenos naturales. Entregando las bases de la disciplina que conocemos hasta el día de hoy. Los avances tecnológicos y las nuevas formas de abordar las problemáticas con las que se enfrentó Matheron en su momento abren una nueva gama de posibilidades y desafíos tanto para los profesionales de las ciencias de la tierra como para la propia disciplina de la geoestadística. En la actualidad los resultados que entregan en muchos casos las técnicas de Machine Learning las podrían convertir en un futuro próximo en el nuevo estándar en el campo del modelaje y la predicción espacial.

00. Instalación de los paquetes requeridos:

#install.packages(c('tidyverse','sf','terra','stars','tmap','remotes','geodata','rstac','earthdatalogin','viridis','gt','quarto','tinytex', 'knitr', 'rmarkdown','gstat','gridExtra','spdep','automap','ranger','tidymodels'))

01. Carga de las librerías requeridas:

library(tidymodels)      # Para manipulación de datos y modelado
library(remotes)         # Para instalar paquetes desde GitHub
library(tidyverse)       # Manipulación y visualización de datos
library(sf)              # Manejo de datos espaciales vectoriales
library(terra)           # Manejo de datos raster
library(tmap)            # Mapas temáticos
library(geodata)         # Datos geográficos
library(rstac)           # Datos STAC
library(earthdatalogin)  # Autenticación en Earthdata
library(viridis)         # Paleta de colores
library(gt)              # Tablas estilizadas
library(quarto)          # Documentos dinámicos
library(tinytex)         # Manejo de LaTeX
library(knitr)           # Generación de informes
library(rmarkdown)       # Documentos de R Markdown
library(gstat)           # Geostatística y variogramas
library(stars)           # Datos espaciales en formato stars
library(gridExtra)       # Organización de gráficos en cuadrículas
library(spdep)           # Análisis de dependencia espacial
library(automap)         # Ajuste automático de variogramas
library(ranger)          # Modelado con Random Forest

02. Preparación de los datos:

El stack solicitado en el segundo punto de los objetivos del taller fue construido con antelación durante el proceso de descarga, procesamiento y ajuste de los datos en el transcurso de las tareas del primer taller. Dicho stack corresponde a uno de los insumos que se cargan antes de comenzar con las tareas propiamente tales del presente ejercicio. Su nombre es predictores.tif y en el código corresponde al objeto “preds”.

# Cargar el objeto 'region' desde el archivo RDS
region <- readRDS("data5/region_RM_UTM.rds")

# Cargar los datos de temperatura con predictores
data_temp <- read_rds('data5/data_estas_temp_con_predictores.rds')

# Cargar los datos raster de predictores (MDE, orientación, pendiente, etc.)
preds <- rast('data5/predictores.tif')  

03. Visualización de las estaciones meteorológicas: (reales y ficticias):

# Agregar una columna que clasifica las estaciones en ficticias o reales según el station_id
data_temp <- data_temp |> 
  mutate(estacion_ficticia = ifelse(station_id >= 1000, "Ficticia", "Real"))

# Visualizar la ubicación de las estaciones de temperatura (reales y ficticias)
tmap_mode('view')  # Configurar el modo interactivo para visualización en tmap
tm_shape(region) +  # Mostrar la región
  tm_borders() +  # Dibujar los bordes de la región
  tm_shape(data_temp) +  # Añadir las estaciones de temperatura
  tm_dots(col = "estacion_ficticia", size = 0.05, 
          palette = c("blue", "red"),  # Colores para estaciones reales y ficticias
          title = "Tipo de Estación") +  # Leyenda que indica el tipo de estación
  tm_layout(legend.outside = TRUE,       # Leyenda colocada fuera del mapa
            legend.title.size = 0.8,     # Tamaño del título de la leyenda
            legend.text.size = 0.6)      # Tamaño del texto de la leyenda

04. Definición de variables predictoras y meses de análisis:

Se definió una lista de variables predictoras por cada mes de interés (Marzo, Junio, Septiembre, Diciembre).

# Definir los meses de interés
meses <- c(3, 6, 9, 12)  # Marzo, Junio, Septiembre, Diciembre
nombre_meses <- c("Marzo", "Junio", "Septiembre", "Diciembre")

# Definir las variables predictoras por mes
pred_vars_list <- list(
  c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Marzo', 'LST_Marzo'),    # Marzo
  c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Junio', 'LST_Junio'),    # Junio
  c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Septiembre', 'LST_Septiembre'),  # Septiembre
  c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Diciembre', 'LST_Diciembre')  # Diciembre
)

05. Calcular buffer de distancia (raster):

En esta sección, se aborda el cálculo de un buffer de distancia en formato raster entre la ubicación de las estaciones meteorológicas y cada uno de los píxeles del área de estudio. Este proceso es fundamental para entender cómo varían las distancias en el dominio espacial en relación con las estaciones. Se generarán tantos buffers como estaciones existan en la región, lo que permitirá analizar la influencia de cada estación en el contexto espacial del área de interés. Este enfoque es especialmente útil para modelar fenómenos geoespaciales, ya que considera la proximidad de las estaciones a los diferentes píxeles, facilitando así la interpolación de datos y la predicción de variables ambientales.

Evaluación de los resultados del modelo: La evaluación de un modelo de predicción es un paso crítico para determinar su eficacia y precisión. En esta parte del análisis, se realizará una evaluación de los resultados del modelo utilizando un enfoque de remuestreo con 10 folds, junto con comparaciones entre el conjunto de entrenamiento y el conjunto de prueba. Esta metodología permite obtener estimaciones más robustas del rendimiento del modelo, ya que se valida en diferentes subconjuntos de los datos. Se compararán métricas clave como el Error Cuadrático Medio (RMSE), el Error Absoluto Medio (MAE) y el coeficiente de determinación (R²). Estas métricas proporcionarán una visión clara de la capacidad del modelo para ajustarse a los datos y su capacidad de generalización en datos no vistos.

05.1 Preparación de datos para Marzo:

Cuando se realiza un proceso de Machine Learning lo primero que se debe tener en cuenta dentro del flujo de trabajo es la división de los datos entre un set de entrenamiento y otro de testeo. Por el lado de los datos de testeo estos una vez divididos solo se vuelven a retomar al final del proceso. Son los datos de entrenamiento los que tienen mayor participación del proceso en su totalidad. A estos datos se les puede realizar un remuestreo que divide el set de entrenamiento y con ello permitir observar la estabilidad del modelo. Entonces, el set de entrenamiento se utiliza para estimar los parámetros del modelo, mientras que el set de testeo se utiliza para encontrar una evaluación independiente del rendimiento de dicho modelo.

El remuestreo permite calcular distintas métricas, pero en distintos subconjuntos de datos y ver qué tan estables son. Hay que tener en cuenta que para trabajos de este tipo el remuestreo tiene que considerar la espacialización, ya que un remuestreo aleatorio podría dejar áreas importantes de la región sin ser consideradas. Una vez realizado el remuestreo se pueden probar diferentes modelos de Machine Learning, en este caso solo se trabajó con Random Forest.

Una vez tengo los datos divididos, seleccionado el modelo, evaluado el modelo en el set de testeo, set de entrenamiento y en el remuestreo, viene el ajuste final. Este ajuste se realiza con la función last_fist que incluye el workflow final del modelo y el data_temp_split. Estas métricas calculadas sobre el set de testeo. Con este ajuste final es posible observar las predicciones de la variable de interés y poder realizar las predicciones. Luego de hacer el ajuste final se utiliza el set de testeo para comprobar la calidad del modelo. Adicionalmente, es posible realizar un ajuste de los parámetros y volver a correr la función last_fist (con los nuevos ajustes).

# Filtrar las capas predictoras por mes
preds_mes <- subset(preds, pred_vars_list[[1]])
preds_mes_st <- st_as_stars(preds_mes) |> split('band')

# Filtrar los datos de temperatura por mes correspondiente
data_temp_filtered_st <- data_temp |> filter(month == 3)

# Obtener la extensión de preds_mes para crear la grilla
ext <- ext(preds_mes)  # Obtener los límites geográficos de preds_mes

# Crear la grilla con resolución de 500 metros a partir de la extensión calculada
grilla <- rast(ext, res = 500, crs = "EPSG:32719")  # Grilla de 500m de resolución
xy <- xyFromCell(grilla, 1:ncell(grilla))
xy_sf <- st_as_sf(data.frame(xy), coords = c('x', 'y'), crs = "EPSG:32719")

# Comprobar resultados
print(grilla)  # Imprimir detalles de la grilla

class       : SpatRaster 
dimensions  : 301, 359, 1  (nrow, ncol, nlyr)
resolution  : 500, 500  (x, y)
extent      : 249135.3, 428635.3, 6205240, 6355740  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 19S (EPSG:32719) 

# Calcular distancias entre cada estación y todos los puntos de la grilla
dists <- lapply(1:nrow(data_temp_filtered_st),\(i){
  d <- st_distance(xy_sf,data_temp_filtered_st[i,])
  values(grilla) <- d
  grilla
})

# Crear el objeto ráster con distancias
buffer_dist <- rast(dists)
names(buffer_dist) <- paste0('D', 1:59)

# Cortar el ráster a la extensión de la región
buffer_dist_crop <- crop(buffer_dist, region)

# Aplicar máscara para ajustar a la forma del polígono
buffer_dist_cortado <- mask(buffer_dist_crop, region)

# Extraer los valores de distancias de 'buffer_dist' para cada estación
buffer_dist_df <- terra::extract(buffer_dist_cortado,data_temp_filtered_st)

# Agregar los datos de distancia a la data de nuestra variable
data_temp_filtered <- cbind(data_temp_filtered_st,buffer_dist_df[,-1])

#Plotear Los 6 primeros predictores distancia como ejemplo
plot(buffer_dist_cortado[[1:6]])

# Selecionamos la variable de interes + los predictores existentes + los nuevos predictores en base a raster de distancia.
data_temp_filtered <- data_temp_filtered |> 
  select(temp_promedio, MDE, Orientación, Pendiente, DistanciaCosta, NDVI_Marzo, LST_Marzo, D1:D59)

# Eliminamos las Coordenadas para tener un data frame limpio ya que tidymodels y RF no perminten otro formato.
data_temp_filtered <- data_temp_filtered |> st_drop_geometry() # para Marzo 3

# División de los datos split

set.seed(123)
data_temp_split <- initial_split(data_temp_filtered)    #La división se realizará con una proporción de 3/4 de los datos
data_temp_split

<Training/Testing/Total>
<44/15/59>

# Extracción de los datos
data_temp_train <- training(data_temp_split) # set de datos de entrenamiento
data_temp_test <- testing(data_temp_split) # set de datos de muestreo

# Receta
tree_rec <- recipe(temp_promedio ~ . , data = data_temp_train) |> 
  step_zv(all_numeric_predictors()) |> 
  step_normalize(all_numeric_predictors())

ctrl <- control_grid(parallel_over = 'everything')

# Definiendo el Modelo
tune_spec <- rand_forest(
  trees = 1000,
  mtry = tune(),# Ajustar mtry
  min_n = tune() # Ajustar min_n
) |> 
  set_mode('regression') |> 
  set_engine('ranger', keep.inbag = TRUE)

# Creando un workflow
tune_wf <- workflow() |> 
  add_recipe(tree_rec) |> 
  add_model(tune_spec) 

# Validación cruzada

set.seed(234)
trees_folds <- vfold_cv(data_temp_train, v = 10)

set.seed(345)
tunes_res <- tune_grid(
  tune_wf,
  resample = trees_folds,
  grid = 16
)
tunes_res

# Tuning results
# 10-fold cross-validation 
# A tibble: 10 × 4
   splits         id     .metrics          .notes          
   <list>         <chr>  <list>            <list>          
 1 <split [39/5]> Fold01 <tibble [32 × 6]> <tibble [0 × 3]>
 2 <split [39/5]> Fold02 <tibble [32 × 6]> <tibble [0 × 3]>
 3 <split [39/5]> Fold03 <tibble [32 × 6]> <tibble [0 × 3]>
 4 <split [39/5]> Fold04 <tibble [32 × 6]> <tibble [0 × 3]>
 5 <split [40/4]> Fold05 <tibble [32 × 6]> <tibble [0 × 3]>
 6 <split [40/4]> Fold06 <tibble [32 × 6]> <tibble [0 × 3]>
 7 <split [40/4]> Fold07 <tibble [32 × 6]> <tibble [0 × 3]>
 8 <split [40/4]> Fold08 <tibble [32 × 6]> <tibble [0 × 3]>
 9 <split [40/4]> Fold09 <tibble [32 × 6]> <tibble [0 × 3]>
10 <split [40/4]> Fold10 <tibble [32 × 6]> <tibble [0 × 3]>

tunes_res |> 
  collect_metrics() |> 
  filter(.metric =='rsq') |> 
  select(mean, min_n, mtry) |> 
  pivot_longer(min_n:mtry,
               values_to = 'value',
               names_to = 'parameter'
  ) |> 
  ggplot(aes(value, mean, color = parameter)) +
  geom_point(show.legend = FALSE)+
  facet_wrap(~parameter, scales = 'free_x')+
  labs(x = NULL, y = 'R2')

#Refinar el Tunnig
rf_grid <- grid_regular(
  mtry(range = c(120,150)),
  min_n(range = c(1,10)),
  levels = 5
)
rf_grid

# A tibble: 25 × 2
    mtry min_n
   <int> <int>
 1   120     1
 2   127     1
 3   135     1
 4   142     1
 5   150     1
 6   120     3
 7   127     3
 8   135     3
 9   142     3
10   150     3
# ℹ 15 more rows

set.seed(456)
regular_res <- tune_grid(
  tune_wf,
  resample = trees_folds,
  grid = rf_grid
)

regular_res 

# Tuning results
# 10-fold cross-validation 
# A tibble: 10 × 4
   splits         id     .metrics          .notes           
   <list>         <chr>  <list>            <list>           
 1 <split [39/5]> Fold01 <tibble [50 × 6]> <tibble [25 × 3]>
 2 <split [39/5]> Fold02 <tibble [50 × 6]> <tibble [25 × 3]>
 3 <split [39/5]> Fold03 <tibble [50 × 6]> <tibble [25 × 3]>
 4 <split [39/5]> Fold04 <tibble [50 × 6]> <tibble [25 × 3]>
 5 <split [40/4]> Fold05 <tibble [50 × 6]> <tibble [25 × 3]>
 6 <split [40/4]> Fold06 <tibble [50 × 6]> <tibble [25 × 3]>
 7 <split [40/4]> Fold07 <tibble [50 × 6]> <tibble [25 × 3]>
 8 <split [40/4]> Fold08 <tibble [50 × 6]> <tibble [25 × 3]>
 9 <split [40/4]> Fold09 <tibble [50 × 6]> <tibble [25 × 3]>
10 <split [40/4]> Fold10 <tibble [50 × 6]> <tibble [25 × 3]>

There were issues with some computations:

  - Warning(s) x50: 120 columns were requested but there were 65 predictors in the da...
  - Warning(s) x50: 127 columns were requested but there were 65 predictors in the da...
  - Warning(s) x50: 135 columns were requested but there were 65 predictors in the da...
  - Warning(s) x50: 142 columns were requested but there were 65 predictors in the da...
  - Warning(s) x50: 150 columns were requested but there were 65 predictors in the da...

Run `show_notes(.Last.tune.result)` for more information.

# Elegir el modelo con la mejor performance
best_rsq <- select_best(regular_res, metric = 'rsq')

# Finalizar el modelo
final_rf <- finalize_model(
  tune_spec,
  best_rsq
)

final_rf

Random Forest Model Specification (regression)

Main Arguments:
  mtry = 142
  trees = 1000
  min_n = 5

Engine-Specific Arguments:
  keep.inbag = TRUE

Computational engine: ranger 

# Creando un workflow Final
final_wf <- workflow() |> 
  add_recipe(tree_rec) |> 
  add_model(final_rf) 


final_res <- final_wf |> 
  last_fit(data_temp_split)


rf_met <- final_res |> 
  collect_metrics()

# Agregar las métricas al data.frame
(metricas_totales_marzo_RF <- rf_met)

# A tibble: 2 × 4
  .metric .estimator .estimate .config             
  <chr>   <chr>          <dbl> <chr>               
1 rmse    standard       1.40  Preprocessor1_Model1
2 rsq     standard       0.946 Preprocessor1_Model1

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_totales_marzo_RF, "data5/procesada/metricas_totales_marzo_RF.rds")

rf_fit <- final_wf |> fit(data_temp_train)

# Evaluar las predicciones

augment(rf_fit, new_data = data_temp_test)

# A tibble: 15 × 67
   .pred temp_promedio   MDE Orientación Pendiente DistanciaCosta NDVI_Marzo
   <dbl>         <dbl> <dbl>       <dbl>     <dbl>          <dbl>      <dbl>
 1  7.03          8.18 4991.       182.     36.5          155170.     0.0270
 2 10.0           8.87 3504.       145.     29.4          116020.     0.0360
 3  7.22          7.37 3833.       189.     17.6          176103.     0.0719
 4 13.1          11.6  3476.       139.     27.9          125842.     0.0694
 5  7.51          9.86 3070.        90.2    23.0          146021.     0.156 
 6 13.2          12.3  2943.       193.     10.4          125087.     0.126 
 7 11.3          12.1  2194.       195.     34.4          139848.     0.209 
 8 18.6          18.2   102.       119.      0.693         25330.     0.468 
 9 20.0          19.9   530.       193.      0.671         85680.     0.426 
10 20.1          19.5   346.       212.      1.04          64891.     0.380 
11 20.9          18.2   546.       239.      0.628         81196.     0.391 
12 18.6          18.8   155.       234.      0.962         40537.     0.445 
13 19.5          18.1   210.       192.      4.83          31664.     0.336 
14 19.5          16.8   182.       248.      1.01          50339.     0.443 
15 20.0          20.0   538.       223.      0.641         85495.     0.413 
# ℹ 60 more variables: LST_Marzo <dbl>, D1 <dbl>, D2 <dbl>, D3 <dbl>, D4 <dbl>,
#   D5 <dbl>, D6 <dbl>, D7 <dbl>, D8 <dbl>, D9 <dbl>, D10 <dbl>, D11 <dbl>,
#   D12 <dbl>, D13 <dbl>, D14 <dbl>, D15 <dbl>, D16 <dbl>, D17 <dbl>,
#   D18 <dbl>, D19 <dbl>, D20 <dbl>, D21 <dbl>, D22 <dbl>, D23 <dbl>,
#   D24 <dbl>, D25 <dbl>, D26 <dbl>, D27 <dbl>, D28 <dbl>, D29 <dbl>,
#   D30 <dbl>, D31 <dbl>, D32 <dbl>, D33 <dbl>, D34 <dbl>, D35 <dbl>,
#   D36 <dbl>, D37 <dbl>, D38 <dbl>, D39 <dbl>, D40 <dbl>, D41 <dbl>, …

data_ev <- augment(rf_fit, new_data = data_temp_test)
ggplot(data_ev,aes(.pred,temp_promedio)) + 
  geom_point() + 
  geom_abline(slope =1,lty='dashed', color = 'red') + 
  theme_bw()

# Evaluar las predicciones en el set de entrenamiento

metrics <- augment(rf_fit, new_data = data_temp_train) |> 
  metrics(truth = .pred, 
          estimate = temp_promedio)

# Agregar las métricas Evaluacion Train al data.frame
(metricas_train_marzo_RF <- metrics)

# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       0.563
2 rsq     standard       0.991
3 mae     standard       0.443

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_train_marzo_RF, "data5/procesada/metricas_train_marzo_RF.rds")


# Evaluar las predicciones en el set de testeo

metrics <- augment(rf_fit, new_data = data_temp_test) |> 
  metrics(truth = .pred, 
          estimate = temp_promedio)

# Agregar las métricas Evaluacion Muestreo al data.frame
(metricas_test_marzo_RF <- metrics)

# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       1.40 
2 rsq     standard       0.944
3 mae     standard       1.08 

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_test_marzo_RF, "data5/procesada/metricas_test_marzo_RF.rds")


# Otra forma de evaluar las predicciones es a traves de un remuestreo     # Sirve para corroborar la estabilidad de las metricas de las predicciones

set.seed(567)
data_temp_folds <- vfold_cv(data_temp_train, v = 10)
data_temp_folds

#  10-fold cross-validation 
# A tibble: 10 × 2
   splits         id    
   <list>         <chr> 
 1 <split [39/5]> Fold01
 2 <split [39/5]> Fold02
 3 <split [39/5]> Fold03
 4 <split [39/5]> Fold04
 5 <split [40/4]> Fold05
 6 <split [40/4]> Fold06
 7 <split [40/4]> Fold07
 8 <split [40/4]> Fold08
 9 <split [40/4]> Fold09
10 <split [40/4]> Fold10

# Usando un workflow Final
final_wf <- workflow() |> 
  add_recipe(tree_rec) |> 
  add_model(final_rf)

data_temp_res <- fit_resamples(final_wf , data_temp_folds)
data_temp_res

# Resampling results
# 10-fold cross-validation 
# A tibble: 10 × 4
   splits         id     .metrics         .notes          
   <list>         <chr>  <list>           <list>          
 1 <split [39/5]> Fold01 <tibble [2 × 4]> <tibble [1 × 3]>
 2 <split [39/5]> Fold02 <tibble [2 × 4]> <tibble [1 × 3]>
 3 <split [39/5]> Fold03 <tibble [2 × 4]> <tibble [1 × 3]>
 4 <split [39/5]> Fold04 <tibble [2 × 4]> <tibble [1 × 3]>
 5 <split [40/4]> Fold05 <tibble [2 × 4]> <tibble [1 × 3]>
 6 <split [40/4]> Fold06 <tibble [2 × 4]> <tibble [1 × 3]>
 7 <split [40/4]> Fold07 <tibble [2 × 4]> <tibble [1 × 3]>
 8 <split [40/4]> Fold08 <tibble [2 × 4]> <tibble [1 × 3]>
 9 <split [40/4]> Fold09 <tibble [2 × 4]> <tibble [1 × 3]>
10 <split [40/4]> Fold10 <tibble [2 × 4]> <tibble [1 × 3]>

There were issues with some computations:

  - Warning(s) x10: 142 columns were requested but there were 65 predictors in the da...

Run `show_notes(.Last.tune.result)` for more information.

data_temp_res |>     
  unnest(.metrics)    # unnest me entrega el detalle de todos los folds

# A tibble: 20 × 7
   splits         id     .metric .estimator .estimate .config           .notes  
   <list>         <chr>  <chr>   <chr>          <dbl> <chr>             <list>  
 1 <split [39/5]> Fold01 rmse    standard       1.68  Preprocessor1_Mo… <tibble>
 2 <split [39/5]> Fold01 rsq     standard       0.938 Preprocessor1_Mo… <tibble>
 3 <split [39/5]> Fold02 rmse    standard       1.83  Preprocessor1_Mo… <tibble>
 4 <split [39/5]> Fold02 rsq     standard       0.833 Preprocessor1_Mo… <tibble>
 5 <split [39/5]> Fold03 rmse    standard       1.51  Preprocessor1_Mo… <tibble>
 6 <split [39/5]> Fold03 rsq     standard       0.985 Preprocessor1_Mo… <tibble>
 7 <split [39/5]> Fold04 rmse    standard       1.21  Preprocessor1_Mo… <tibble>
 8 <split [39/5]> Fold04 rsq     standard       0.961 Preprocessor1_Mo… <tibble>
 9 <split [40/4]> Fold05 rmse    standard       1.09  Preprocessor1_Mo… <tibble>
10 <split [40/4]> Fold05 rsq     standard       0.980 Preprocessor1_Mo… <tibble>
11 <split [40/4]> Fold06 rmse    standard       1.59  Preprocessor1_Mo… <tibble>
12 <split [40/4]> Fold06 rsq     standard       0.818 Preprocessor1_Mo… <tibble>
13 <split [40/4]> Fold07 rmse    standard       1.12  Preprocessor1_Mo… <tibble>
14 <split [40/4]> Fold07 rsq     standard       0.943 Preprocessor1_Mo… <tibble>
15 <split [40/4]> Fold08 rmse    standard       1.51  Preprocessor1_Mo… <tibble>
16 <split [40/4]> Fold08 rsq     standard       0.924 Preprocessor1_Mo… <tibble>
17 <split [40/4]> Fold09 rmse    standard       1.45  Preprocessor1_Mo… <tibble>
18 <split [40/4]> Fold09 rsq     standard       1.00  Preprocessor1_Mo… <tibble>
19 <split [40/4]> Fold10 rmse    standard       2.24  Preprocessor1_Mo… <tibble>
20 <split [40/4]> Fold10 rsq     standard       0.987 Preprocessor1_Mo… <tibble>

metrics <-data_temp_res |> 
  collect_metrics()   # Collect_metrics me muestra el promedio de las métricas de todo los folds

# Agregar las métricas Evaluacion Muestreo al data.frame
(metricas_vfold_cv_marzo_RF <- metrics)

# A tibble: 2 × 6
  .metric .estimator  mean     n std_err .config             
  <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
1 rmse    standard   1.52     10  0.110  Preprocessor1_Model1
2 rsq     standard   0.937    10  0.0201 Preprocessor1_Model1

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_vfold_cv_marzo_RF, "data5/procesada/metricas_vfold_cv_marzo_RF.rds")


# Combinar predictores y distancias
preds_marzo <- c(preds_mes, buffer_dist)

names(preds_marzo)

 [1] "MDE"            "Orientación"    "Pendiente"      "DistanciaCosta"
 [5] "NDVI_Marzo"     "LST_Marzo"      "D1"             "D2"            
 [9] "D3"             "D4"             "D5"             "D6"            
[13] "D7"             "D8"             "D9"             "D10"           
[17] "D11"            "D12"            "D13"            "D14"           
[21] "D15"            "D16"            "D17"            "D18"           
[25] "D19"            "D20"            "D21"            "D22"           
[29] "D23"            "D24"            "D25"            "D26"           
[33] "D27"            "D28"            "D29"            "D30"           
[37] "D31"            "D32"            "D33"            "D34"           
[41] "D35"            "D36"            "D37"            "D38"           
[45] "D39"            "D40"            "D41"            "D42"           
[49] "D43"            "D44"            "D45"            "D46"           
[53] "D47"            "D48"            "D49"            "D50"           
[57] "D51"            "D52"            "D53"            "D54"           
[61] "D55"            "D56"            "D57"            "D58"           
[65] "D59"           

# Convertir a tibble
preds_val <- values(preds_marzo) |> as_tibble() |> 
  mutate(
    across(1:6, \(x) replace_na(x, 1)),  # Reemplaza los NA con 1 para las primeras variables
  )

# Renombrar columnas
names(preds_val)[1:6] <- c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Marzo', 'LST_Marzo')

# Realizar la predicción con el modelo de bosque aleatorio
pred_df_RF <- predict(rf_fit, preds_val)

# Asignar los valores de predicción a la grilla
values(grilla) <- pred_df_RF$.pred

# Cortar el ráster a la extensión de la región
temp_pred_rast <- crop(grilla, st_as_sf(region))

# Aplicar máscara para ajustar a la forma del polígono
temp_pred_rast_marzo_RF <- mask(grilla, st_as_sf(region))

# Guardar el raster como archivo TIFF
writeRaster(temp_pred_rast_marzo_RF, "data5/procesada/RF_marzo_predicciones_temp_promedio.tif", overwrite = TRUE)

names(temp_pred_rast_marzo_RF) <- 'Predicción Marzo'

# Visualizar el resultado usando tmap
tm_shape(temp_pred_rast_marzo_RF) +
  tm_raster(style = "quantile", palette = viridis::viridis(20)) +
  tm_layout(main.title = "Predicción de Temperatura con RF")

05.2 Preparación de datos para Junio:

# Filtrar las capas predictoras por mes
preds_mes <- subset(preds, pred_vars_list[[2]])
preds_mes_st <- st_as_stars(preds_mes) |> split('band')

# Filtrar los datos de temperatura por mes correspondiente
data_temp_filtered_st <- data_temp |> filter(month == 6)

# Obtener la extensión de preds_mes para crear la grilla
ext <- ext(preds_mes)  # Obtener los límites geográficos de preds_mes

# Crear la grilla con resolución de 500 metros a partir de la extensión calculada
grilla <- rast(ext, res = 500, crs = "EPSG:32719")  # Grilla de 500m de resolución
xy <- xyFromCell(grilla, 1:ncell(grilla))
xy_sf <- st_as_sf(data.frame(xy), coords = c('x', 'y'), crs = "EPSG:32719")

# Comprobar resultados
print(grilla)  # Imprimir detalles de la grilla

class       : SpatRaster 
dimensions  : 301, 359, 1  (nrow, ncol, nlyr)
resolution  : 500, 500  (x, y)
extent      : 249135.3, 428635.3, 6205240, 6355740  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 19S (EPSG:32719) 

# Calcular distancias entre cada estación y todos los puntos de la grilla
dists <- lapply(1:nrow(data_temp_filtered_st),\(i){
  d <- st_distance(xy_sf,data_temp_filtered_st[i,])
  values(grilla) <- d
  grilla
})

# Crear el objeto ráster con distancias
buffer_dist <- rast(dists)
names(buffer_dist) <- paste0('D', 1:58)

# Cortar el ráster a la extensión de la región
buffer_dist_crop <- crop(buffer_dist, region)

# Aplicar máscara para ajustar a la forma del polígono
buffer_dist_cortado <- mask(buffer_dist_crop, region)

# Extraer los valores de distancias de 'buffer_dist' para cada estación
buffer_dist_df <- terra::extract(buffer_dist_cortado,data_temp_filtered_st)

# Agregar los datos de distancia a la data de nuestra variable
data_temp_filtered <- cbind(data_temp_filtered_st,buffer_dist_df[,-1])

#Plotear Los 6 primeros predictores distancia como ejemplo
plot(buffer_dist_cortado[[1:6]])

# Selecionamos la variable de interes + los predictores existentes + los nuevos predictores en base a raster de distancia.
data_temp_filtered <- data_temp_filtered |> 
  select(temp_promedio, MDE, Orientación, Pendiente, DistanciaCosta, NDVI_Junio, LST_Junio, D1:D58)

# Eliminamos las Coordenadas para tener un data frame limpio ya que tidymodels y RF no perminten otro formato.
data_temp_filtered <- data_temp_filtered |> st_drop_geometry() # para Junio 6

# División de los datos split

set.seed(123)
data_temp_split <- initial_split(data_temp_filtered)    #La división se realizará con una proporción de 3/4 de los datos
data_temp_split

<Training/Testing/Total>
<43/15/58>

# Extracción de los datos
data_temp_train <- training(data_temp_split) # set de datos de entrenamiento
data_temp_test <- testing(data_temp_split) # set de datos de muestreo

# Receta
tree_rec <- recipe(temp_promedio ~ . , data = data_temp_train) |> 
  step_zv(all_numeric_predictors()) |> 
  step_normalize(all_numeric_predictors())

ctrl <- control_grid(parallel_over = 'everything')

# Definiendo el Modelo
tune_spec <- rand_forest(
  trees = 1000,
  mtry = tune(),# Ajustar mtry
  min_n = tune() # Ajustar min_n
) |> 
  set_mode('regression') |> 
  set_engine('ranger', keep.inbag = TRUE)

# Creando un workflow
tune_wf <- workflow() |> 
  add_recipe(tree_rec) |> 
  add_model(tune_spec) 

# Validación cruzada

set.seed(234)
trees_folds <- vfold_cv(data_temp_train, v = 10)

set.seed(345)
tunes_res <- tune_grid(
  tune_wf,
  resample = trees_folds,
  grid = 16
)
tunes_res

# Tuning results
# 10-fold cross-validation 
# A tibble: 10 × 4
   splits         id     .metrics          .notes          
   <list>         <chr>  <list>            <list>          
 1 <split [38/5]> Fold01 <tibble [32 × 6]> <tibble [1 × 3]>
 2 <split [38/5]> Fold02 <tibble [32 × 6]> <tibble [1 × 3]>
 3 <split [38/5]> Fold03 <tibble [32 × 6]> <tibble [1 × 3]>
 4 <split [39/4]> Fold04 <tibble [32 × 6]> <tibble [0 × 3]>
 5 <split [39/4]> Fold05 <tibble [32 × 6]> <tibble [0 × 3]>
 6 <split [39/4]> Fold06 <tibble [32 × 6]> <tibble [0 × 3]>
 7 <split [39/4]> Fold07 <tibble [32 × 6]> <tibble [0 × 3]>
 8 <split [39/4]> Fold08 <tibble [32 × 6]> <tibble [0 × 3]>
 9 <split [39/4]> Fold09 <tibble [32 × 6]> <tibble [0 × 3]>
10 <split [39/4]> Fold10 <tibble [32 × 6]> <tibble [0 × 3]>

There were issues with some computations:

  - Warning(s) x3: A correlation computation is required, but `estimate` is constant...

Run `show_notes(.Last.tune.result)` for more information.

tunes_res |> 
  collect_metrics() |> 
  filter(.metric =='rsq') |> 
  select(mean, min_n, mtry) |> 
  pivot_longer(min_n:mtry,
               values_to = 'value',
               names_to = 'parameter'
  ) |> 
  ggplot(aes(value, mean, color = parameter)) +
  geom_point(show.legend = FALSE)+
  facet_wrap(~parameter, scales = 'free_x')+
  labs(x = NULL, y = 'R2')

#Refinar el Tunnig
rf_grid <- grid_regular(
  mtry(range = c(120,150)),
  min_n(range = c(1,10)),
  levels = 5
)
rf_grid

# A tibble: 25 × 2
    mtry min_n
   <int> <int>
 1   120     1
 2   127     1
 3   135     1
 4   142     1
 5   150     1
 6   120     3
 7   127     3
 8   135     3
 9   142     3
10   150     3
# ℹ 15 more rows

set.seed(456)
regular_res <- tune_grid(
  tune_wf,
  resample = trees_folds,
  grid = rf_grid
)

regular_res 

# Tuning results
# 10-fold cross-validation 
# A tibble: 10 × 4
   splits         id     .metrics          .notes           
   <list>         <chr>  <list>            <list>           
 1 <split [38/5]> Fold01 <tibble [50 × 6]> <tibble [25 × 3]>
 2 <split [38/5]> Fold02 <tibble [50 × 6]> <tibble [25 × 3]>
 3 <split [38/5]> Fold03 <tibble [50 × 6]> <tibble [25 × 3]>
 4 <split [39/4]> Fold04 <tibble [50 × 6]> <tibble [25 × 3]>
 5 <split [39/4]> Fold05 <tibble [50 × 6]> <tibble [25 × 3]>
 6 <split [39/4]> Fold06 <tibble [50 × 6]> <tibble [25 × 3]>
 7 <split [39/4]> Fold07 <tibble [50 × 6]> <tibble [25 × 3]>
 8 <split [39/4]> Fold08 <tibble [50 × 6]> <tibble [25 × 3]>
 9 <split [39/4]> Fold09 <tibble [50 × 6]> <tibble [25 × 3]>
10 <split [39/4]> Fold10 <tibble [50 × 6]> <tibble [25 × 3]>

There were issues with some computations:

  - Warning(s) x50: 120 columns were requested but there were 64 predictors in the da...
  - Warning(s) x50: 127 columns were requested but there were 64 predictors in the da...
  - Warning(s) x50: 135 columns were requested but there were 64 predictors in the da...
  - Warning(s) x50: 142 columns were requested but there were 64 predictors in the da...
  - Warning(s) x50: 150 columns were requested but there were 64 predictors in the da...

Run `show_notes(.Last.tune.result)` for more information.

# Elegir el modelo con la mejor performance
best_rsq <- select_best(regular_res, metric = 'rsq')

# Finalizar el modelo
final_rf <- finalize_model(
  tune_spec,
  best_rsq
)

final_rf

Random Forest Model Specification (regression)

Main Arguments:
  mtry = 127
  trees = 1000
  min_n = 1

Engine-Specific Arguments:
  keep.inbag = TRUE

Computational engine: ranger 

# Creando un workflow Final
final_wf <- workflow() |> 
  add_recipe(tree_rec) |> 
  add_model(final_rf) 


final_res <- final_wf |> 
  last_fit(data_temp_split)


rf_met <- final_res |> 
  collect_metrics()

# Agregar las métricas al data.frame
(metricas_totales_junio_RF <- rf_met)

# A tibble: 2 × 4
  .metric .estimator .estimate .config             
  <chr>   <chr>          <dbl> <chr>               
1 rmse    standard       1.48  Preprocessor1_Model1
2 rsq     standard       0.961 Preprocessor1_Model1

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_totales_junio_RF, "data5/procesada/metricas_totales_junio_RF.rds")

rf_fit <- final_wf |> fit(data_temp_train)

# Evaluar las predicciones

augment(rf_fit, new_data = data_temp_test)

# A tibble: 15 × 66
   .pred temp_promedio   MDE Orientación Pendiente DistanciaCosta NDVI_Junio
   <dbl>         <dbl> <dbl>       <dbl>     <dbl>          <dbl>      <dbl>
 1 -7.67         -7.58 4991.       182.     36.5          155170.   -0.0438 
 2 -2.51         -3.17 3504.       145.     29.4          116020.    0.0210 
 3  1.52         -1.72 3476.       139.     27.9          125842.   -0.0178 
 4 -7.11         -5.31 3070.        90.2    23.0          146021.    0.00170
 5  1.57         -1.08 2943.       193.     10.4          125087.   -0.0323 
 6 -3.06         -2.75 2194.       195.     34.4          139848.    0.326  
 7 -2.84         -3.87 2206.       297.     28.1          144998.    0.280  
 8  9.07          9.92  553.       196.      0.785         66952.    0.260  
 9  8.22          9.45  667.       215.      1.04          95408.    0.462  
10 10.4          10.8   530.       193.      0.671         85680.    0.358  
11  9.91         10.9   439.       281.      1.48          95222.    0.507  
12 10.7          10.5   470.       228.      0.834         75014.    0.431  
13 11.3          11.6   155.       234.      0.962         40537.    0.464  
14 10.4          10.5   538.       223.      0.641         85495.    0.339  
15  8.57         11.6   721.       269.      1.51          83577.    0.284  
# ℹ 59 more variables: LST_Junio <dbl>, D1 <dbl>, D2 <dbl>, D3 <dbl>, D4 <dbl>,
#   D5 <dbl>, D6 <dbl>, D7 <dbl>, D8 <dbl>, D9 <dbl>, D10 <dbl>, D11 <dbl>,
#   D12 <dbl>, D13 <dbl>, D14 <dbl>, D15 <dbl>, D16 <dbl>, D17 <dbl>,
#   D18 <dbl>, D19 <dbl>, D20 <dbl>, D21 <dbl>, D22 <dbl>, D23 <dbl>,
#   D24 <dbl>, D25 <dbl>, D26 <dbl>, D27 <dbl>, D28 <dbl>, D29 <dbl>,
#   D30 <dbl>, D31 <dbl>, D32 <dbl>, D33 <dbl>, D34 <dbl>, D35 <dbl>,
#   D36 <dbl>, D37 <dbl>, D38 <dbl>, D39 <dbl>, D40 <dbl>, D41 <dbl>, …

data_ev <- augment(rf_fit, new_data = data_temp_test)
ggplot(data_ev,aes(.pred,temp_promedio)) + 
  geom_point() + 
  geom_abline(slope =1,lty='dashed', color = 'red') + 
  theme_bw()

# Evaluar las predicciones en el set de entrenamiento

metrics <- augment(rf_fit, new_data = data_temp_train) |> 
  metrics(truth = .pred, 
          estimate = temp_promedio)

# Agregar las métricas Evaluacion Train al data.frame
(metricas_train_junio_RF <- metrics)

# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       0.535
2 rsq     standard       0.996
3 mae     standard       0.393

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_train_junio_RF, "data5/procesada/metricas_train_junio_RF.rds")


# Evaluar las predicciones en el set de testeo

metrics <- augment(rf_fit, new_data = data_temp_test) |> 
  metrics(truth = .pred, 
          estimate = temp_promedio)

# Agregar las métricas Evaluacion Muestreo al data.frame
(metricas_test_junio_RF <- metrics)

# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       1.53 
2 rsq     standard       0.958
3 mae     standard       1.12 

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_test_junio_RF, "data5/procesada/metricas_test_junio_RF.rds")


# Otra forma de evaluar las predicciones es a traves de un remuestreo     # Sirve para corroborar la estabilidad de las metricas de las predicciones

set.seed(567)
data_temp_folds <- vfold_cv(data_temp_train, v = 10)
data_temp_folds

#  10-fold cross-validation 
# A tibble: 10 × 2
   splits         id    
   <list>         <chr> 
 1 <split [38/5]> Fold01
 2 <split [38/5]> Fold02
 3 <split [38/5]> Fold03
 4 <split [39/4]> Fold04
 5 <split [39/4]> Fold05
 6 <split [39/4]> Fold06
 7 <split [39/4]> Fold07
 8 <split [39/4]> Fold08
 9 <split [39/4]> Fold09
10 <split [39/4]> Fold10

# Usando un workflow Final
final_wf <- workflow() |> 
  add_recipe(tree_rec) |> 
  add_model(final_rf)

data_temp_res <- fit_resamples(final_wf , data_temp_folds)
data_temp_res

# Resampling results
# 10-fold cross-validation 
# A tibble: 10 × 4
   splits         id     .metrics         .notes          
   <list>         <chr>  <list>           <list>          
 1 <split [38/5]> Fold01 <tibble [2 × 4]> <tibble [1 × 3]>
 2 <split [38/5]> Fold02 <tibble [2 × 4]> <tibble [1 × 3]>
 3 <split [38/5]> Fold03 <tibble [2 × 4]> <tibble [1 × 3]>
 4 <split [39/4]> Fold04 <tibble [2 × 4]> <tibble [1 × 3]>
 5 <split [39/4]> Fold05 <tibble [2 × 4]> <tibble [1 × 3]>
 6 <split [39/4]> Fold06 <tibble [2 × 4]> <tibble [1 × 3]>
 7 <split [39/4]> Fold07 <tibble [2 × 4]> <tibble [1 × 3]>
 8 <split [39/4]> Fold08 <tibble [2 × 4]> <tibble [1 × 3]>
 9 <split [39/4]> Fold09 <tibble [2 × 4]> <tibble [1 × 3]>
10 <split [39/4]> Fold10 <tibble [2 × 4]> <tibble [1 × 3]>

There were issues with some computations:

  - Warning(s) x10: 127 columns were requested but there were 64 predictors in the da...

Run `show_notes(.Last.tune.result)` for more information.

data_temp_res |>     
  unnest(.metrics)    # unnest me entrega el detalle de todos los folds

# A tibble: 20 × 7
   splits         id     .metric .estimator .estimate .config           .notes  
   <list>         <chr>  <chr>   <chr>          <dbl> <chr>             <list>  
 1 <split [38/5]> Fold01 rmse    standard     1.37    Preprocessor1_Mo… <tibble>
 2 <split [38/5]> Fold01 rsq     standard     0.980   Preprocessor1_Mo… <tibble>
 3 <split [38/5]> Fold02 rmse    standard     0.807   Preprocessor1_Mo… <tibble>
 4 <split [38/5]> Fold02 rsq     standard     0.991   Preprocessor1_Mo… <tibble>
 5 <split [38/5]> Fold03 rmse    standard     2.37    Preprocessor1_Mo… <tibble>
 6 <split [38/5]> Fold03 rsq     standard     0.923   Preprocessor1_Mo… <tibble>
 7 <split [39/4]> Fold04 rmse    standard     1.70    Preprocessor1_Mo… <tibble>
 8 <split [39/4]> Fold04 rsq     standard     0.00247 Preprocessor1_Mo… <tibble>
 9 <split [39/4]> Fold05 rmse    standard     0.865   Preprocessor1_Mo… <tibble>
10 <split [39/4]> Fold05 rsq     standard     0.983   Preprocessor1_Mo… <tibble>
11 <split [39/4]> Fold06 rmse    standard     1.71    Preprocessor1_Mo… <tibble>
12 <split [39/4]> Fold06 rsq     standard     0.968   Preprocessor1_Mo… <tibble>
13 <split [39/4]> Fold07 rmse    standard     1.52    Preprocessor1_Mo… <tibble>
14 <split [39/4]> Fold07 rsq     standard     0.977   Preprocessor1_Mo… <tibble>
15 <split [39/4]> Fold08 rmse    standard     1.78    Preprocessor1_Mo… <tibble>
16 <split [39/4]> Fold08 rsq     standard     0.956   Preprocessor1_Mo… <tibble>
17 <split [39/4]> Fold09 rmse    standard     1.11    Preprocessor1_Mo… <tibble>
18 <split [39/4]> Fold09 rsq     standard     0.982   Preprocessor1_Mo… <tibble>
19 <split [39/4]> Fold10 rmse    standard     0.453   Preprocessor1_Mo… <tibble>
20 <split [39/4]> Fold10 rsq     standard     1.00    Preprocessor1_Mo… <tibble>

metrics <-data_temp_res |> 
  collect_metrics()   # Collect_metrics me muestra el promedio de las métricas de todo los folds

# Agregar las métricas Evaluacion Muestreo al data.frame
(metricas_vfold_cv_junio_RF <- metrics)

# A tibble: 2 × 6
  .metric .estimator  mean     n std_err .config             
  <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
1 rmse    standard   1.37     10  0.180  Preprocessor1_Model1
2 rsq     standard   0.876    10  0.0973 Preprocessor1_Model1

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_vfold_cv_junio_RF, "data5/procesada/metricas_vfold_cv_junio_RF.rds")


# Combinar predictores y distancias
preds_junio <- c(preds_mes, buffer_dist)

names(preds)

 [1] "MDE"             "Orientación"     "Pendiente"       "DistanciaCosta" 
 [5] "NDVI_Marzo"      "NDVI_Junio"      "NDVI_Septiembre" "NDVI_Diciembre" 
 [9] "LST_Marzo"       "LST_Junio"       "LST_Septiembre"  "LST_Diciembre"  

# Convertir a tibble
preds_val <- values(preds_junio) |> as_tibble() |> 
  mutate(
    across(1:6, \(x) replace_na(x, 1)),  # Reemplaza los NA con 1 para las primeras variables
  )

# Renombrar columnas
names(preds_val)[1:6] <- c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Junio', 'LST_Junio')

# Realizar la predicción con el modelo de bosque aleatorio
pred_df_RF <- predict(rf_fit, preds_val)

# Asignar los valores de predicción a la grilla
values(grilla) <- pred_df_RF$.pred

# Cortar el ráster a la extensión de la región
temp_pred_rast <- crop(grilla, st_as_sf(region))

# Aplicar máscara para ajustar a la forma del polígono
temp_pred_rast_junio_RF <- mask(grilla, st_as_sf(region))

# Guardar el raster como archivo TIFF
writeRaster(temp_pred_rast_junio_RF, "data5/procesada/RF_junio_predicciones_temp_promedio.tif", overwrite = TRUE)

names(temp_pred_rast_junio_RF) <- 'Predicción Junio'

# Visualizar el resultado usando tmap
tm_shape(temp_pred_rast_junio_RF) +
  tm_raster(style = "quantile", palette = viridis::viridis(20)) +
  tm_layout(main.title = "Predicción de Temperatura con RF")

05.3 Preparación de datos para Septiembre:

# Filtrar las capas predictoras por mes

preds_mes <- subset(preds, pred_vars_list[[3]])
preds_mes_st <- st_as_stars(preds_mes) |> split('band')

# Filtrar los datos de temperatura por mes correspondiente
data_temp_filtered_st <- data_temp |> filter(month == 9)

# Obtener la extensión de preds_mes para crear la grilla
ext <- ext(preds_mes)  # Obtener los límites geográficos de preds_mes

# Crear la grilla con resolución de 500 metros a partir de la extensión calculada
grilla <- rast(ext, res = 500, crs = "EPSG:32719")  # Grilla de 500m de resolución
xy <- xyFromCell(grilla, 1:ncell(grilla))
xy_sf <- st_as_sf(data.frame(xy), coords = c('x', 'y'), crs = "EPSG:32719")

# Comprobar resultados
print(grilla)  # Imprimir detalles de la grilla

class       : SpatRaster 
dimensions  : 301, 359, 1  (nrow, ncol, nlyr)
resolution  : 500, 500  (x, y)
extent      : 249135.3, 428635.3, 6205240, 6355740  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 19S (EPSG:32719) 

# Calcular distancias entre cada estación y todos los puntos de la grilla
dists <- lapply(1:nrow(data_temp_filtered_st),\(i){
  d <- st_distance(xy_sf,data_temp_filtered_st[i,])
  values(grilla) <- d
  grilla
})

# Crear el objeto ráster con distancias
buffer_dist <- rast(dists)
names(buffer_dist) <- paste0('D', 1:57)

# Cortar el ráster a la extensión de la región
buffer_dist_crop <- crop(buffer_dist, region)

# Aplicar máscara para ajustar a la forma del polígono
buffer_dist_cortado <- mask(buffer_dist_crop, region)

# Extraer los valores de distancias de 'buffer_dist' para cada estación
buffer_dist_df <- terra::extract(buffer_dist_cortado,data_temp_filtered_st)

# Agregar los datos de distancia a la data de nuestra variable
data_temp_filtered <- cbind(data_temp_filtered_st,buffer_dist_df[,-1])

#Plotear Los 6 primeros predictores distancia como ejemplo
plot(buffer_dist_cortado[[1:6]])

# Selecionamos la variable de interes + los predictores existentes + los nuevos predictores en base a raster de distancia.
data_temp_filtered <- data_temp_filtered |> 
  select(temp_promedio, MDE, Orientación, Pendiente, DistanciaCosta, NDVI_Septiembre, LST_Septiembre, D1:D57)

# Eliminamos las Coordenadas para tener un data frame limpio ya que tidymodels y RF no perminten otro formato.
data_temp_filtered <- data_temp_filtered |> st_drop_geometry() # para Septiembre 9

# División de los datos split

set.seed(123)
data_temp_split <- initial_split(data_temp_filtered)    #La división se realizará con una proporción de 3/4 de los datos
data_temp_split

<Training/Testing/Total>
<42/15/57>

# Extracción de los datos
data_temp_train <- training(data_temp_split) # set de datos de entrenamiento
data_temp_test <- testing(data_temp_split) # set de datos de muestreo

# Receta
tree_rec <- recipe(temp_promedio ~ . , data = data_temp_train) |> 
  step_zv(all_numeric_predictors()) |> 
  step_normalize(all_numeric_predictors())

ctrl <- control_grid(parallel_over = 'everything')

# Definiendo el Modelo
tune_spec <- rand_forest(
  trees = 1000,
  mtry = tune(),# Ajustar mtry
  min_n = tune() # Ajustar min_n
) |> 
  set_mode('regression') |> 
  set_engine('ranger', keep.inbag = TRUE)

# Creando un workflow
tune_wf <- workflow() |> 
  add_recipe(tree_rec) |> 
  add_model(tune_spec) 

# Validación cruzada

set.seed(234)
trees_folds <- vfold_cv(data_temp_train, v = 10)

set.seed(345)
tunes_res <- tune_grid(
  tune_wf,
  resample = trees_folds,
  grid = 16
)
tunes_res

# Tuning results
# 10-fold cross-validation 
# A tibble: 10 × 4
   splits         id     .metrics          .notes          
   <list>         <chr>  <list>            <list>          
 1 <split [37/5]> Fold01 <tibble [32 × 6]> <tibble [2 × 3]>
 2 <split [37/5]> Fold02 <tibble [32 × 6]> <tibble [2 × 3]>
 3 <split [38/4]> Fold03 <tibble [32 × 6]> <tibble [1 × 3]>
 4 <split [38/4]> Fold04 <tibble [32 × 6]> <tibble [1 × 3]>
 5 <split [38/4]> Fold05 <tibble [32 × 6]> <tibble [1 × 3]>
 6 <split [38/4]> Fold06 <tibble [32 × 6]> <tibble [1 × 3]>
 7 <split [38/4]> Fold07 <tibble [32 × 6]> <tibble [1 × 3]>
 8 <split [38/4]> Fold08 <tibble [32 × 6]> <tibble [1 × 3]>
 9 <split [38/4]> Fold09 <tibble [32 × 6]> <tibble [1 × 3]>
10 <split [38/4]> Fold10 <tibble [32 × 6]> <tibble [1 × 3]>

There were issues with some computations:

  - Warning(s) x2: 38 samples were requested but there were 37 rows in the data. 37 ...
  - Warning(s) x10: A correlation computation is required, but `estimate` is constant...

Run `show_notes(.Last.tune.result)` for more information.

tunes_res |> 
  collect_metrics() |> 
  filter(.metric =='rsq') |> 
  select(mean, min_n, mtry) |> 
  pivot_longer(min_n:mtry,
               values_to = 'value',
               names_to = 'parameter'
  ) |> 
  ggplot(aes(value, mean, color = parameter)) +
  geom_point(show.legend = FALSE)+
  facet_wrap(~parameter, scales = 'free_x')+
  labs(x = NULL, y = 'R2')

#Refinar el Tunnig
rf_grid <- grid_regular(
  mtry(range = c(120,150)),
  min_n(range = c(1,10)),
  levels = 5
)
rf_grid

# A tibble: 25 × 2
    mtry min_n
   <int> <int>
 1   120     1
 2   127     1
 3   135     1
 4   142     1
 5   150     1
 6   120     3
 7   127     3
 8   135     3
 9   142     3
10   150     3
# ℹ 15 more rows

set.seed(456)
regular_res <- tune_grid(
  tune_wf,
  resample = trees_folds,
  grid = rf_grid
)

regular_res 

# Tuning results
# 10-fold cross-validation 
# A tibble: 10 × 4
   splits         id     .metrics          .notes           
   <list>         <chr>  <list>            <list>           
 1 <split [37/5]> Fold01 <tibble [50 × 6]> <tibble [25 × 3]>
 2 <split [37/5]> Fold02 <tibble [50 × 6]> <tibble [25 × 3]>
 3 <split [38/4]> Fold03 <tibble [50 × 6]> <tibble [25 × 3]>
 4 <split [38/4]> Fold04 <tibble [50 × 6]> <tibble [25 × 3]>
 5 <split [38/4]> Fold05 <tibble [50 × 6]> <tibble [25 × 3]>
 6 <split [38/4]> Fold06 <tibble [50 × 6]> <tibble [25 × 3]>
 7 <split [38/4]> Fold07 <tibble [50 × 6]> <tibble [25 × 3]>
 8 <split [38/4]> Fold08 <tibble [50 × 6]> <tibble [25 × 3]>
 9 <split [38/4]> Fold09 <tibble [50 × 6]> <tibble [25 × 3]>
10 <split [38/4]> Fold10 <tibble [50 × 6]> <tibble [25 × 3]>

There were issues with some computations:

  - Warning(s) x50: 120 columns were requested but there were 63 predictors in the da...
  - Warning(s) x50: 127 columns were requested but there were 63 predictors in the da...
  - Warning(s) x50: 135 columns were requested but there were 63 predictors in the da...
  - Warning(s) x50: 142 columns were requested but there were 63 predictors in the da...
  - Warning(s) x50: 150 columns were requested but there were 63 predictors in the da...

Run `show_notes(.Last.tune.result)` for more information.

# Elegir el modelo con la mejor performance
best_rsq <- select_best(regular_res, metric = 'rsq')

# Finalizar el modelo
final_rf <- finalize_model(
  tune_spec,
  best_rsq
)

final_rf

Random Forest Model Specification (regression)

Main Arguments:
  mtry = 150
  trees = 1000
  min_n = 10

Engine-Specific Arguments:
  keep.inbag = TRUE

Computational engine: ranger 

# Creando un workflow Final
final_wf <- workflow() |> 
  add_recipe(tree_rec) |> 
  add_model(final_rf) 


final_res <- final_wf |> 
  last_fit(data_temp_split)


rf_met <- final_res |> 
  collect_metrics()

# Agregar las métricas al data.frame
(metricas_totales_septiembre_RF <- rf_met)

# A tibble: 2 × 4
  .metric .estimator .estimate .config             
  <chr>   <chr>          <dbl> <chr>               
1 rmse    standard       1.27  Preprocessor1_Model1
2 rsq     standard       0.972 Preprocessor1_Model1

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_totales_septiembre_RF, "data5/procesada/metricas_totales_septiembre_RF.rds")

rf_fit <- final_wf |> fit(data_temp_train)

# Evaluar las predicciones

augment(rf_fit, new_data = data_temp_test)

# A tibble: 15 × 65
    .pred temp_promedio   MDE Orientación Pendiente DistanciaCosta
    <dbl>         <dbl> <dbl>       <dbl>     <dbl>          <dbl>
 1 -7.54         -7.03  4991.       182.     36.5          155170.
 2 -3.54         -3.66  3504.       145.     29.4          116020.
 3  0.438        -1.62  3476.       139.     27.9          125842.
 4 -7.00         -5.23  3070.        90.2    23.0          146021.
 5  0.646        -0.969 2943.       193.     10.4          125087.
 6 -3.36         -2.46  2194.       195.     34.4          139848.
 7 10.8          11.8    553.       196.      0.785         66952.
 8 10.0          11.1    667.       215.      1.04          95408.
 9 12.0          12.0    530.       193.      0.671         85680.
10 12.2          12.0    470.       228.      0.834         75014.
11 10.9          10.8    627.       222.      0.813         90594.
12 14.1          11.8    143.       179.      0.664         28408.
13 12.1          11.7    538.       223.      0.641         85495.
14  9.71         11.6    721.       269.      1.51          83577.
15 11.9          12.1    492.       139.      0.769         73461.
# ℹ 59 more variables: NDVI_Septiembre <dbl>, LST_Septiembre <dbl>, D1 <dbl>,
#   D2 <dbl>, D3 <dbl>, D4 <dbl>, D5 <dbl>, D6 <dbl>, D7 <dbl>, D8 <dbl>,
#   D9 <dbl>, D10 <dbl>, D11 <dbl>, D12 <dbl>, D13 <dbl>, D14 <dbl>, D15 <dbl>,
#   D16 <dbl>, D17 <dbl>, D18 <dbl>, D19 <dbl>, D20 <dbl>, D21 <dbl>,
#   D22 <dbl>, D23 <dbl>, D24 <dbl>, D25 <dbl>, D26 <dbl>, D27 <dbl>,
#   D28 <dbl>, D29 <dbl>, D30 <dbl>, D31 <dbl>, D32 <dbl>, D33 <dbl>,
#   D34 <dbl>, D35 <dbl>, D36 <dbl>, D37 <dbl>, D38 <dbl>, D39 <dbl>, …

data_ev <- augment(rf_fit, new_data = data_temp_test)
ggplot(data_ev,aes(.pred,temp_promedio)) + 
  geom_point() + 
  geom_abline(slope =1,lty='dashed', color = 'red') + 
  theme_bw()

# Evaluar las predicciones en el set de entrenamiento

metrics <- augment(rf_fit, new_data = data_temp_train) |> 
  metrics(truth = .pred, 
          estimate = temp_promedio)

# Agregar las métricas Evaluacion Train al data.frame
(metricas_train_septiembre_RF <- metrics)

# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       0.674
2 rsq     standard       0.995
3 mae     standard       0.499

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_train_septiembre_RF, "data5/procesada/metricas_train_septiembre_RF.rds")


# Evaluar las predicciones en el set de testeo

metrics <- augment(rf_fit, new_data = data_temp_test) |> 
  metrics(truth = .pred, 
          estimate = temp_promedio)

# Agregar las métricas Evaluacion Muestreo al data.frame
(metricas_test_septiembre_RF <- metrics)

# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       1.22 
2 rsq     standard       0.975
3 mae     standard       0.936

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_test_septiembre_RF, "data5/procesada/metricas_test_septiembre_RF.rds")


# Otra forma de evaluar las predicciones es a traves de un remuestreo     # Sirve para corroborar la estabilidad de las metricas de las predicciones

set.seed(567)
data_temp_folds <- vfold_cv(data_temp_train, v = 10)
data_temp_folds

#  10-fold cross-validation 
# A tibble: 10 × 2
   splits         id    
   <list>         <chr> 
 1 <split [37/5]> Fold01
 2 <split [37/5]> Fold02
 3 <split [38/4]> Fold03
 4 <split [38/4]> Fold04
 5 <split [38/4]> Fold05
 6 <split [38/4]> Fold06
 7 <split [38/4]> Fold07
 8 <split [38/4]> Fold08
 9 <split [38/4]> Fold09
10 <split [38/4]> Fold10

# Usando un workflow Final
final_wf <- workflow() |> 
  add_recipe(tree_rec) |> 
  add_model(final_rf)

data_temp_res <- fit_resamples(final_wf , data_temp_folds)
data_temp_res

# Resampling results
# 10-fold cross-validation 
# A tibble: 10 × 4
   splits         id     .metrics         .notes          
   <list>         <chr>  <list>           <list>          
 1 <split [37/5]> Fold01 <tibble [2 × 4]> <tibble [1 × 3]>
 2 <split [37/5]> Fold02 <tibble [2 × 4]> <tibble [1 × 3]>
 3 <split [38/4]> Fold03 <tibble [2 × 4]> <tibble [1 × 3]>
 4 <split [38/4]> Fold04 <tibble [2 × 4]> <tibble [1 × 3]>
 5 <split [38/4]> Fold05 <tibble [2 × 4]> <tibble [1 × 3]>
 6 <split [38/4]> Fold06 <tibble [2 × 4]> <tibble [1 × 3]>
 7 <split [38/4]> Fold07 <tibble [2 × 4]> <tibble [1 × 3]>
 8 <split [38/4]> Fold08 <tibble [2 × 4]> <tibble [1 × 3]>
 9 <split [38/4]> Fold09 <tibble [2 × 4]> <tibble [1 × 3]>
10 <split [38/4]> Fold10 <tibble [2 × 4]> <tibble [1 × 3]>

There were issues with some computations:

  - Warning(s) x10: 150 columns were requested but there were 63 predictors in the da...

Run `show_notes(.Last.tune.result)` for more information.

data_temp_res |>     
  unnest(.metrics)    # unnest me entrega el detalle de todos los folds

# A tibble: 20 × 7
   splits         id     .metric .estimator .estimate .config           .notes  
   <list>         <chr>  <chr>   <chr>          <dbl> <chr>             <list>  
 1 <split [37/5]> Fold01 rmse    standard       1.73  Preprocessor1_Mo… <tibble>
 2 <split [37/5]> Fold01 rsq     standard       0.988 Preprocessor1_Mo… <tibble>
 3 <split [37/5]> Fold02 rmse    standard       1.76  Preprocessor1_Mo… <tibble>
 4 <split [37/5]> Fold02 rsq     standard       0.890 Preprocessor1_Mo… <tibble>
 5 <split [38/4]> Fold03 rmse    standard       3.17  Preprocessor1_Mo… <tibble>
 6 <split [38/4]> Fold03 rsq     standard       0.998 Preprocessor1_Mo… <tibble>
 7 <split [38/4]> Fold04 rmse    standard       1.24  Preprocessor1_Mo… <tibble>
 8 <split [38/4]> Fold04 rsq     standard       0.994 Preprocessor1_Mo… <tibble>
 9 <split [38/4]> Fold05 rmse    standard       0.715 Preprocessor1_Mo… <tibble>
10 <split [38/4]> Fold05 rsq     standard       0.997 Preprocessor1_Mo… <tibble>
11 <split [38/4]> Fold06 rmse    standard       1.87  Preprocessor1_Mo… <tibble>
12 <split [38/4]> Fold06 rsq     standard       0.961 Preprocessor1_Mo… <tibble>
13 <split [38/4]> Fold07 rmse    standard       3.07  Preprocessor1_Mo… <tibble>
14 <split [38/4]> Fold07 rsq     standard       0.989 Preprocessor1_Mo… <tibble>
15 <split [38/4]> Fold08 rmse    standard       1.54  Preprocessor1_Mo… <tibble>
16 <split [38/4]> Fold08 rsq     standard       0.900 Preprocessor1_Mo… <tibble>
17 <split [38/4]> Fold09 rmse    standard       0.636 Preprocessor1_Mo… <tibble>
18 <split [38/4]> Fold09 rsq     standard       0.224 Preprocessor1_Mo… <tibble>
19 <split [38/4]> Fold10 rmse    standard       0.613 Preprocessor1_Mo… <tibble>
20 <split [38/4]> Fold10 rsq     standard       0.994 Preprocessor1_Mo… <tibble>

metrics <-data_temp_res |> 
  collect_metrics()   # Collect_metrics me muestra el promedio de las métricas de todo los folds

# Agregar las métricas Evaluacion Muestreo al data.frame
(metricas_vfold_cv_septiembre_RF <- metrics)

# A tibble: 2 × 6
  .metric .estimator  mean     n std_err .config             
  <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
1 rmse    standard   1.63     10  0.290  Preprocessor1_Model1
2 rsq     standard   0.893    10  0.0755 Preprocessor1_Model1

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_vfold_cv_septiembre_RF, "data5/procesada/metricas_vfold_cv_septiembre_RF.rds")


# Combinar predictores y distancias
preds_septiembre <- c(preds_mes, buffer_dist)

names(preds)

 [1] "MDE"             "Orientación"     "Pendiente"       "DistanciaCosta" 
 [5] "NDVI_Marzo"      "NDVI_Junio"      "NDVI_Septiembre" "NDVI_Diciembre" 
 [9] "LST_Marzo"       "LST_Junio"       "LST_Septiembre"  "LST_Diciembre"  

# Convertir a tibble
preds_val <- values(preds_septiembre) |> as_tibble() |> 
  mutate(
    across(1:6, \(x) replace_na(x, 1)),  # Reemplaza los NA con 1 para las primeras variables
  )

# Renombrar columnas
names(preds_val)[1:6] <- c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Septiembre', 'LST_Septiembre')

# Realizar la predicción con el modelo de bosque aleatorio
pred_df_RF <- predict(rf_fit, preds_val)

# Asignar los valores de predicción a la grilla
values(grilla) <- pred_df_RF$.pred

# Cortar el ráster a la extensión de la región
temp_pred_rast <- crop(grilla, st_as_sf(region))

# Aplicar máscara para ajustar a la forma del polígono
temp_pred_rast_septiembre_RF <- mask(grilla, st_as_sf(region))

# Guardar el raster como archivo TIFF
writeRaster(temp_pred_rast_septiembre_RF, "data5/procesada/RF_septiembre_predicciones_temp_promedio.tif", overwrite = TRUE)

names(temp_pred_rast_septiembre_RF) <- 'Predicción Septiembre'

# Visualizar el resultado usando tmap
tm_shape(temp_pred_rast_septiembre_RF) +
  tm_raster(style = "quantile", palette = viridis::viridis(20)) +
  tm_layout(main.title = "Predicción de Temperatura con RF")

05.4 Preparación de datos para Diciembre:

# Filtrar las capas predictoras por mes
preds_mes <- subset(preds, pred_vars_list[[4]])
preds_mes_st <- st_as_stars(preds_mes) |> split('band')

# Filtrar los datos de temperatura por mes correspondiente
data_temp_filtered_st <- data_temp |> filter(month == 12)

# Obtener la extensión de preds_mes para crear la grilla
ext <- ext(preds_mes)  # Obtener los límites geográficos de preds_mes

# Crear la grilla con resolución de 500 metros a partir de la extensión calculada
grilla <- rast(ext, res = 500, crs = "EPSG:32719")  # Grilla de 500m de resolución
xy <- xyFromCell(grilla, 1:ncell(grilla))
xy_sf <- st_as_sf(data.frame(xy), coords = c('x', 'y'), crs = "EPSG:32719")

# Comprobar resultados
print(grilla)  # Imprimir detalles de la grilla

class       : SpatRaster 
dimensions  : 301, 359, 1  (nrow, ncol, nlyr)
resolution  : 500, 500  (x, y)
extent      : 249135.3, 428635.3, 6205240, 6355740  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 19S (EPSG:32719) 

# Calcular distancias entre cada estación y todos los puntos de la grilla
dists <- lapply(1:nrow(data_temp_filtered_st),\(i){
  d <- st_distance(xy_sf,data_temp_filtered_st[i,])
  values(grilla) <- d
  grilla
})

# Crear el objeto ráster con distancias
buffer_dist <- rast(dists)
names(buffer_dist) <- paste0('D', 1:58)

# Cortar el ráster a la extensión de la región
buffer_dist_crop <- crop(buffer_dist, region)

# Aplicar máscara para ajustar a la forma del polígono
buffer_dist_cortado <- mask(buffer_dist_crop, region)

# Extraer los valores de distancias de 'buffer_dist' para cada estación
buffer_dist_df <- terra::extract(buffer_dist_cortado,data_temp_filtered_st)

# Agregar los datos de distancia a la data de nuestra variable
data_temp_filtered <- cbind(data_temp_filtered_st,buffer_dist_df[,-1])

#Plotear Los 6 primeros predictores distancia como ejemplo
plot(buffer_dist_cortado[[1:6]])

# Selecionamos la variable de interes + los predictores existentes + los nuevos predictores en base a raster de distancia.
data_temp_filtered <- data_temp_filtered |> 
  select(temp_promedio, MDE, Orientación, Pendiente, DistanciaCosta, NDVI_Diciembre, LST_Diciembre, D1:D58)

# Eliminamos las Coordenadas para tener un data frame limpio ya que tidymodels y RF no perminten otro formato.
data_temp_filtered <- data_temp_filtered |> st_drop_geometry() # para Diciembre 12

# División de los datos split

set.seed(123)
data_temp_split <- initial_split(data_temp_filtered)    #La división se realizará con una proporción de 3/4 de los datos
data_temp_split

<Training/Testing/Total>
<43/15/58>

# Extracción de los datos
data_temp_train <- training(data_temp_split) # set de datos de entrenamiento
data_temp_test <- testing(data_temp_split) # set de datos de muestreo

# Receta
tree_rec <- recipe(temp_promedio ~ . , data = data_temp_train) |> 
  step_zv(all_numeric_predictors()) |> 
  step_normalize(all_numeric_predictors())

ctrl <- control_grid(parallel_over = 'everything')

# Definiendo el Modelo
tune_spec <- rand_forest(
  trees = 1000,
  mtry = tune(),# Ajustar mtry
  min_n = tune() # Ajustar min_n
) |> 
  set_mode('regression') |> 
  set_engine('ranger', keep.inbag = TRUE)

# Creando un workflow
tune_wf <- workflow() |> 
  add_recipe(tree_rec) |> 
  add_model(tune_spec) 

# Validación cruzada

set.seed(234)
trees_folds <- vfold_cv(data_temp_train, v = 10)

set.seed(345)
tunes_res <- tune_grid(
  tune_wf,
  resample = trees_folds,
  grid = 15
)
tunes_res

# Tuning results
# 10-fold cross-validation 
# A tibble: 10 × 4
   splits         id     .metrics          .notes          
   <list>         <chr>  <list>            <list>          
 1 <split [38/5]> Fold01 <tibble [30 × 6]> <tibble [2 × 3]>
 2 <split [38/5]> Fold02 <tibble [30 × 6]> <tibble [2 × 3]>
 3 <split [38/5]> Fold03 <tibble [30 × 6]> <tibble [2 × 3]>
 4 <split [39/4]> Fold04 <tibble [30 × 6]> <tibble [2 × 3]>
 5 <split [39/4]> Fold05 <tibble [30 × 6]> <tibble [2 × 3]>
 6 <split [39/4]> Fold06 <tibble [30 × 6]> <tibble [2 × 3]>
 7 <split [39/4]> Fold07 <tibble [30 × 6]> <tibble [2 × 3]>
 8 <split [39/4]> Fold08 <tibble [30 × 6]> <tibble [2 × 3]>
 9 <split [39/4]> Fold09 <tibble [30 × 6]> <tibble [2 × 3]>
10 <split [39/4]> Fold10 <tibble [30 × 6]> <tibble [2 × 3]>

There were issues with some computations:

  - Warning(s) x3: 40 samples were requested but there were 38 rows in the data. 38 ...
  - Warning(s) x7: 40 samples were requested but there were 39 rows in the data. 39 ...
  - Warning(s) x10: A correlation computation is required, but `estimate` is constant...

Run `show_notes(.Last.tune.result)` for more information.

tunes_res |> 
  collect_metrics() |> 
  filter(.metric =='rsq') |> 
  select(mean, min_n, mtry) |> 
  pivot_longer(min_n:mtry,
               values_to = 'value',
               names_to = 'parameter'
  ) |> 
  ggplot(aes(value, mean, color = parameter)) +
  geom_point(show.legend = FALSE)+
  facet_wrap(~parameter, scales = 'free_x')+
  labs(x = NULL, y = 'R2')

#Refinar el Tunnig
rf_grid <- grid_regular(
  mtry(range = c(120,150)),
  min_n(range = c(1,10)),
  levels = 5
)
rf_grid

# A tibble: 25 × 2
    mtry min_n
   <int> <int>
 1   120     1
 2   127     1
 3   135     1
 4   142     1
 5   150     1
 6   120     3
 7   127     3
 8   135     3
 9   142     3
10   150     3
# ℹ 15 more rows

set.seed(456)
regular_res <- tune_grid(
  tune_wf,
  resample = trees_folds,
  grid = rf_grid
)

regular_res 

# Tuning results
# 10-fold cross-validation 
# A tibble: 10 × 4
   splits         id     .metrics          .notes           
   <list>         <chr>  <list>            <list>           
 1 <split [38/5]> Fold01 <tibble [50 × 6]> <tibble [25 × 3]>
 2 <split [38/5]> Fold02 <tibble [50 × 6]> <tibble [25 × 3]>
 3 <split [38/5]> Fold03 <tibble [50 × 6]> <tibble [25 × 3]>
 4 <split [39/4]> Fold04 <tibble [50 × 6]> <tibble [25 × 3]>
 5 <split [39/4]> Fold05 <tibble [50 × 6]> <tibble [25 × 3]>
 6 <split [39/4]> Fold06 <tibble [50 × 6]> <tibble [25 × 3]>
 7 <split [39/4]> Fold07 <tibble [50 × 6]> <tibble [25 × 3]>
 8 <split [39/4]> Fold08 <tibble [50 × 6]> <tibble [25 × 3]>
 9 <split [39/4]> Fold09 <tibble [50 × 6]> <tibble [25 × 3]>
10 <split [39/4]> Fold10 <tibble [50 × 6]> <tibble [25 × 3]>

There were issues with some computations:

  - Warning(s) x50: 120 columns were requested but there were 64 predictors in the da...
  - Warning(s) x50: 127 columns were requested but there were 64 predictors in the da...
  - Warning(s) x50: 135 columns were requested but there were 64 predictors in the da...
  - Warning(s) x50: 142 columns were requested but there were 64 predictors in the da...
  - Warning(s) x50: 150 columns were requested but there were 64 predictors in the da...

Run `show_notes(.Last.tune.result)` for more information.

# Elegir el modelo con la mejor performance
best_rsq <- select_best(regular_res, metric = 'rsq')

# Finalizar el modelo
final_rf <- finalize_model(
  tune_spec,
  best_rsq
)

final_rf

Random Forest Model Specification (regression)

Main Arguments:
  mtry = 142
  trees = 1000
  min_n = 3

Engine-Specific Arguments:
  keep.inbag = TRUE

Computational engine: ranger 

# Creando un workflow Final
final_wf <- workflow() |> 
  add_recipe(tree_rec) |> 
  add_model(final_rf) 


final_res <- final_wf |> 
  last_fit(data_temp_split)


rf_met <- final_res |> 
  collect_metrics()

# Agregar las métricas al data.frame
(metricas_totales_diciembre_RF <- rf_met)

# A tibble: 2 × 4
  .metric .estimator .estimate .config             
  <chr>   <chr>          <dbl> <chr>               
1 rmse    standard       1.61  Preprocessor1_Model1
2 rsq     standard       0.950 Preprocessor1_Model1

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_totales_diciembre_RF, "data5/procesada/metricas_totales_diciembre_RF.rds")

rf_fit <- final_wf |> fit(data_temp_train)

# Evaluar las predicciones

augment(rf_fit, new_data = data_temp_test)

# A tibble: 15 × 66
   .pred temp_promedio   MDE Orientación Pendiente DistanciaCosta NDVI_Diciembre
   <dbl>         <dbl> <dbl>       <dbl>     <dbl>          <dbl>          <dbl>
 1  2.09          2.90 4991.       182.     36.5          155170.        -0.0651
 2  6.24          5.52 3504.       145.     29.4          116020.         0.0313
 3 10.8           8.36 3476.       139.     27.9          125842.        -0.0766
 4  2.71          4.52 3070.        90.2    23.0          146021.         0.132 
 5 10.9           9.14 2943.       193.     10.4          125087.        -0.0535
 6  7.33          7.57 2194.       195.     34.4          139848.         0.264 
 7  7.43          5.89 2206.       297.     28.1          144998.         0.339 
 8 19.4          20.7   553.       196.      0.785         66952.         0.277 
 9 15.3          18.7   667.       215.      1.04          95408.         0.522 
10 18.5          18.8   530.       193.      0.671         85680.         0.486 
11 17.9          19.8   439.       281.      1.48          95222.         0.563 
12 18.7          18.6   470.       228.      0.834         75014.         0.511 
13 18.3          18.1   155.       234.      0.962         40537.         0.472 
14 18.5          18.7   538.       223.      0.641         85495.         0.512 
15 18.0          20.7   721.       269.      1.51          83577.         0.307 
# ℹ 59 more variables: LST_Diciembre <dbl>, D1 <dbl>, D2 <dbl>, D3 <dbl>,
#   D4 <dbl>, D5 <dbl>, D6 <dbl>, D7 <dbl>, D8 <dbl>, D9 <dbl>, D10 <dbl>,
#   D11 <dbl>, D12 <dbl>, D13 <dbl>, D14 <dbl>, D15 <dbl>, D16 <dbl>,
#   D17 <dbl>, D18 <dbl>, D19 <dbl>, D20 <dbl>, D21 <dbl>, D22 <dbl>,
#   D23 <dbl>, D24 <dbl>, D25 <dbl>, D26 <dbl>, D27 <dbl>, D28 <dbl>,
#   D29 <dbl>, D30 <dbl>, D31 <dbl>, D32 <dbl>, D33 <dbl>, D34 <dbl>,
#   D35 <dbl>, D36 <dbl>, D37 <dbl>, D38 <dbl>, D39 <dbl>, D40 <dbl>, …

data_ev <- augment(rf_fit, new_data = data_temp_test)
ggplot(data_ev,aes(.pred,temp_promedio)) + 
  geom_point() + 
  geom_abline(slope =1,lty='dashed', color = 'red') + 
  theme_bw()

# Evaluar las predicciones en el set de entrenamiento

metrics <- augment(rf_fit, new_data = data_temp_train) |> 
  metrics(truth = .pred, 
          estimate = temp_promedio)

# Agregar las métricas Evaluacion Train al data.frame
(metricas_train_diciembre_RF <- metrics)

# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       0.552
2 rsq     standard       0.995
3 mae     standard       0.383

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_train_diciembre_RF, "data5/procesada/metricas_train_diciembre_RF.rds")


# Evaluar las predicciones en el set de testeo

metrics <- augment(rf_fit, new_data = data_temp_test) |> 
  metrics(truth = .pred, 
          estimate = temp_promedio)

# Agregar las métricas Evaluacion Muestreo al data.frame
(metricas_test_diciembre_RF <- metrics)

# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       1.64 
2 rsq     standard       0.946
3 mae     standard       1.30 

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_test_diciembre_RF, "data5/procesada/metricas_test_diciembre_RF.rds")


# Otra forma de evaluar las predicciones es a traves de un remuestreo     # Sirve para corroborar la estabilidad de las metricas de las predicciones

set.seed(567)
data_temp_folds <- vfold_cv(data_temp_train, v = 10)
data_temp_folds

#  10-fold cross-validation 
# A tibble: 10 × 2
   splits         id    
   <list>         <chr> 
 1 <split [38/5]> Fold01
 2 <split [38/5]> Fold02
 3 <split [38/5]> Fold03
 4 <split [39/4]> Fold04
 5 <split [39/4]> Fold05
 6 <split [39/4]> Fold06
 7 <split [39/4]> Fold07
 8 <split [39/4]> Fold08
 9 <split [39/4]> Fold09
10 <split [39/4]> Fold10

# Usando un workflow Final
final_wf <- workflow() |> 
  add_recipe(tree_rec) |> 
  add_model(final_rf)

data_temp_res <- fit_resamples(final_wf , data_temp_folds)
data_temp_res

# Resampling results
# 10-fold cross-validation 
# A tibble: 10 × 4
   splits         id     .metrics         .notes          
   <list>         <chr>  <list>           <list>          
 1 <split [38/5]> Fold01 <tibble [2 × 4]> <tibble [1 × 3]>
 2 <split [38/5]> Fold02 <tibble [2 × 4]> <tibble [1 × 3]>
 3 <split [38/5]> Fold03 <tibble [2 × 4]> <tibble [1 × 3]>
 4 <split [39/4]> Fold04 <tibble [2 × 4]> <tibble [1 × 3]>
 5 <split [39/4]> Fold05 <tibble [2 × 4]> <tibble [1 × 3]>
 6 <split [39/4]> Fold06 <tibble [2 × 4]> <tibble [1 × 3]>
 7 <split [39/4]> Fold07 <tibble [2 × 4]> <tibble [1 × 3]>
 8 <split [39/4]> Fold08 <tibble [2 × 4]> <tibble [1 × 3]>
 9 <split [39/4]> Fold09 <tibble [2 × 4]> <tibble [1 × 3]>
10 <split [39/4]> Fold10 <tibble [2 × 4]> <tibble [1 × 3]>

There were issues with some computations:

  - Warning(s) x10: 142 columns were requested but there were 64 predictors in the da...

Run `show_notes(.Last.tune.result)` for more information.

data_temp_res |>     
  unnest(.metrics)    # unnest me entrega el detalle de todos los folds

# A tibble: 20 × 7
   splits         id     .metric .estimator .estimate .config           .notes  
   <list>         <chr>  <chr>   <chr>          <dbl> <chr>             <list>  
 1 <split [38/5]> Fold01 rmse    standard       1.32  Preprocessor1_Mo… <tibble>
 2 <split [38/5]> Fold01 rsq     standard       0.976 Preprocessor1_Mo… <tibble>
 3 <split [38/5]> Fold02 rmse    standard       0.634 Preprocessor1_Mo… <tibble>
 4 <split [38/5]> Fold02 rsq     standard       0.993 Preprocessor1_Mo… <tibble>
 5 <split [38/5]> Fold03 rmse    standard       2.92  Preprocessor1_Mo… <tibble>
 6 <split [38/5]> Fold03 rsq     standard       0.818 Preprocessor1_Mo… <tibble>
 7 <split [39/4]> Fold04 rmse    standard       1.43  Preprocessor1_Mo… <tibble>
 8 <split [39/4]> Fold04 rsq     standard       0.550 Preprocessor1_Mo… <tibble>
 9 <split [39/4]> Fold05 rmse    standard       0.888 Preprocessor1_Mo… <tibble>
10 <split [39/4]> Fold05 rsq     standard       0.976 Preprocessor1_Mo… <tibble>
11 <split [39/4]> Fold06 rmse    standard       1.16  Preprocessor1_Mo… <tibble>
12 <split [39/4]> Fold06 rsq     standard       0.978 Preprocessor1_Mo… <tibble>
13 <split [39/4]> Fold07 rmse    standard       1.37  Preprocessor1_Mo… <tibble>
14 <split [39/4]> Fold07 rsq     standard       0.976 Preprocessor1_Mo… <tibble>
15 <split [39/4]> Fold08 rmse    standard       2.08  Preprocessor1_Mo… <tibble>
16 <split [39/4]> Fold08 rsq     standard       0.987 Preprocessor1_Mo… <tibble>
17 <split [39/4]> Fold09 rmse    standard       1.24  Preprocessor1_Mo… <tibble>
18 <split [39/4]> Fold09 rsq     standard       0.975 Preprocessor1_Mo… <tibble>
19 <split [39/4]> Fold10 rmse    standard       0.941 Preprocessor1_Mo… <tibble>
20 <split [39/4]> Fold10 rsq     standard       0.996 Preprocessor1_Mo… <tibble>

metrics <-data_temp_res |> 
  collect_metrics()   # Collect_metrics me muestra el promedio de las métricas de todo los folds

# Agregar las métricas Evaluacion Muestreo al data.frame
(metricas_vfold_cv_diciembre_RF <- metrics)

# A tibble: 2 × 6
  .metric .estimator  mean     n std_err .config             
  <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
1 rmse    standard   1.40     10  0.209  Preprocessor1_Model1
2 rsq     standard   0.923    10  0.0445 Preprocessor1_Model1

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_vfold_cv_diciembre_RF, "data5/procesada/metricas_vfold_cv_diciembre_RF.rds")


# Combinar predictores y distancias
preds_diciembre <- c(preds_mes, buffer_dist)

names(preds)

 [1] "MDE"             "Orientación"     "Pendiente"       "DistanciaCosta" 
 [5] "NDVI_Marzo"      "NDVI_Junio"      "NDVI_Septiembre" "NDVI_Diciembre" 
 [9] "LST_Marzo"       "LST_Junio"       "LST_Septiembre"  "LST_Diciembre"  

# Convertir a tibble
preds_val <- values(preds_diciembre) |> as_tibble() |> 
  mutate(
    across(1:6, \(x) replace_na(x, 1)),  # Reemplaza los NA con 1 para las primeras variables
  )

# Renombrar columnas
names(preds_val)[1:6] <- c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Diciembre', 'LST_Diciembre')

# Realizar la predicción con el modelo de bosque aleatorio
pred_df_RF <- predict(rf_fit, preds_val)

# Asignar los valores de predicción a la grilla
values(grilla) <- pred_df_RF$.pred

# Cortar el ráster a la extensión de la región
temp_pred_rast <- crop(grilla, st_as_sf(region))

# Aplicar máscara para ajustar a la forma del polígono
temp_pred_rast_diciembre_RF <- mask(grilla, st_as_sf(region))

# Guardar el raster como archivo TIFF
writeRaster(temp_pred_rast_diciembre_RF, "data5/procesada/RF_diciembre_predicciones_temp_promedio.tif", overwrite = TRUE)

names(temp_pred_rast_diciembre_RF) <- 'Predicción Diciembre'

# Visualizar el resultado usando tmap
tm_shape(temp_pred_rast_diciembre_RF) +
  tm_raster(style = "quantile", palette = viridis::viridis(20)) +
  tm_layout(main.title = "Predicción de Temperatura con RF")

Los datos de entrenamiento del mes de marzo muestran que el modelo se ajusta bien a los datos de entrenamiento con un bajo error, por su parte los datos de prueba aunque tienen un error mayor, su R2 sigue siendo alto. Para el caso de junio el modelo se ajusta muy bien a los datos de entrenamiento con muy buenos valores en todas sus métricas. Al igual que en marzo el R2 para los datos de prueba indica un muy buen ajuste.

Al observar septiembre es posible identificar que los datos siguen con la misma tendencia de los meses anteriores y se destaca de manera positiva en términos de rendimiento en el conjunto de prueba, con el RMSE más bajo y un alto valor de R2. Finalmente, el caso de diciembre muestra el peor desempeño en términos de error en el conjunto de prueba, con un RMSE elevado y el MAE más elevado de todo el conjunto de datos.

06. Evaluacion modelo de Random Forest (RF) para la temperatura de los meses asignados:

En esta sección, se evaluó el modelo de Random Forest (RF) para predecir la temperatura en marzo, junio, septiembre y diciembre. Se utilizaron métricas de rendimiento como el Error Cuadrático Medio (RMSE) y el coeficiente de determinación (R²) para analizar la precisión del modelo. Las métricas de cada mes se combinaron en un data.frame y se almacenaron en un archivo RDS para un análisis posterior. Además, se generaron gráficos de barras para visualizar el RMSE y el R², facilitando la identificación de tendencias en el desempeño del modelo a lo largo de los meses. Esta evaluación ofreció una comprensión más clara del comportamiento del modelo RF en diferentes condiciones temporales y ayuda a identificar el mejor ajuste para la predicción de temperaturas en la región.

# Crear el data.frame combinado con las métricas de cada mes para el Modelo predictivo RF, aplicado
(metricas_totales_RF <- bind_rows(
  Marzo = metricas_totales_marzo_RF,
  Junio = metricas_totales_junio_RF,
  Septiembre = metricas_totales_septiembre_RF,
  Diciembre = metricas_totales_diciembre_RF,
  .id = "Mes"
))

# A tibble: 8 × 5
  Mes        .metric .estimator .estimate .config             
  <chr>      <chr>   <chr>          <dbl> <chr>               
1 Marzo      rmse    standard       1.40  Preprocessor1_Model1
2 Marzo      rsq     standard       0.946 Preprocessor1_Model1
3 Junio      rmse    standard       1.48  Preprocessor1_Model1
4 Junio      rsq     standard       0.961 Preprocessor1_Model1
5 Septiembre rmse    standard       1.27  Preprocessor1_Model1
6 Septiembre rsq     standard       0.972 Preprocessor1_Model1
7 Diciembre  rmse    standard       1.61  Preprocessor1_Model1
8 Diciembre  rsq     standard       0.950 Preprocessor1_Model1

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_totales_RF, "data5/procesada/metricas_totales_RF.rds")


# Filtrar para RMSE
metricas_rmse <- metricas_totales_RF |> 
  filter(.metric == "rmse")

# Gráfico comparativo de RMSE entre meses
ggplot(metricas_rmse, aes(x = Mes, y = .estimate, fill = Mes)) +
  geom_bar(stat = "identity", position = "dodge", na.rm = TRUE) +
  geom_text(aes(label = round(.estimate, 2)), position = position_dodge(width = 0.9), vjust = -0.3, na.rm = TRUE) +
  ggtitle(" Modelo Predictivo RF -- Comparación de RMSE entre Modelos por Mes") +
  xlab("MESES") +
  ylab("RMSE") +
  scale_fill_manual(name = "Mes", values = c("Marzo" = "steelblue", "Junio" = "orange", "Septiembre" = "green", "Diciembre" = "red")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Filtrar para R²
metricas_rsq <- metricas_totales_RF |> 
  filter(.metric == "rsq")

# Gráfico comparativo de R² entre meses
ggplot(metricas_rsq, aes(x = Mes, y = .estimate, fill = Mes)) +
  geom_bar(stat = "identity", position = "dodge", na.rm = TRUE) +
  geom_text(aes(label = round(.estimate, 2)), position = position_dodge(width = 0.9), vjust = -0.3, na.rm = TRUE) +
  ggtitle("Modelo Predictivo RF -- Comparación de R² entre Modelos por Mes") +
  xlab("MESES") +
  ylab("R²") +
  scale_fill_manual(name = "Mes", values = c("Marzo" = "steelblue", "Junio" = "orange", "Septiembre" = "green", "Diciembre" = "red")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Graficar manteniendo el orden de los modelos por mes
data_FR_long <- metricas_totales_RF |> 
  pivot_longer(cols = .estimate, names_to = "Métrica", values_to = "Valor") 


# Modelo RF (RMSE y R² entre meses)
plot_RF <- ggplot(metricas_totales_RF, aes(x = Mes, y = .estimate, fill = .metric)) +
  geom_col(alpha = 0.7, na.rm = TRUE, position = "dodge") +  # Usar position = "dodge" para separar las barras
  ggtitle("Modelo Predictivo RF -- (RMSE y R² entre meses)") +
  xlab("MESES") +
  ylab("Valor de la Métrica") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_fill_manual(name = "Métrica", values = c("rmse" = "steelblue", "rsq" = "orange"))  # Colores por métrica

# Mostrar el gráfico
plot_RF

07. Evaluación de los resultados del modelo considerando: Set de entrenamiento:

En esta sección, se evaluó el rendimiento del modelo de Random Forest (RF) en el set de entrenamiento. Se combinaron las métricas de rendimiento, incluyendo RMSE, R² y MAE, para los meses de marzo, junio, septiembre y diciembre. Este análisis es esencial para entender la capacidad predictiva del modelo en los datos utilizados para su entrenamiento.Se presentará un gráfico de barras que comparará las métricas entre los diferentes meses, facilitando la visualización de cómo varía el rendimiento del modelo. Este enfoque no solo valida la calidad del modelo, sino que también sienta las bases para futuras comparaciones con el set de testeo y otras configuraciones.

# Crear el data.frame combinado con las métricas Train del Modelo predictivo RF.
(metricas_train_RF <- bind_rows(
  marzo = metricas_train_marzo_RF,
  junio = metricas_train_junio_RF,
  septiembre = metricas_train_septiembre_RF,
  diciembre = metricas_train_diciembre_RF,
  .id = "Mes"
))

# A tibble: 12 × 4
   Mes        .metric .estimator .estimate
   <chr>      <chr>   <chr>          <dbl>
 1 marzo      rmse    standard       0.563
 2 marzo      rsq     standard       0.991
 3 marzo      mae     standard       0.443
 4 junio      rmse    standard       0.535
 5 junio      rsq     standard       0.996
 6 junio      mae     standard       0.393
 7 septiembre rmse    standard       0.674
 8 septiembre rsq     standard       0.995
 9 septiembre mae     standard       0.499
10 diciembre  rmse    standard       0.552
11 diciembre  rsq     standard       0.995
12 diciembre  mae     standard       0.383

# Graficar manteniendo el orden de los modelos por mes
data_train_FR_long <- metricas_train_RF |> 
  pivot_longer(cols = .estimate, names_to = "Métrica", values_to = "Valor") 

# Train RF (RMSE y R² y MAE entre meses)
ggplot(metricas_train_RF, aes(x = Mes, y = .estimate, fill = .metric)) +
  geom_col(alpha = 0.7, na.rm = TRUE, position = "dodge") +  # Usar position = "dodge" para separar las barras
  ggtitle("Evaluacion Set de Entrenamiento RF -- (RMSE y R² y MAE entre meses)") +
  xlab("MESES") +
  ylab("Valor de la Métrica") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_fill_manual(name = "Métrica", values = c("rmse" = "steelblue", "rsq" = "orange", 'mae'= 'grey'))  # Colores por métrica

08. Evaluación de los resultados del modelo considerando: Set de muestreo:

En esta sección, se evaluó el rendimiento del modelo de Random Forest (RF) utilizando el set de muestreo. Se combinaron las métricas de rendimiento, incluyendo RMSE, R² y MAE, para los meses de marzo, junio, septiembre y diciembre. Esta evaluación es crucial para determinar la capacidad del modelo para generalizar a datos no vistos durante el entrenamiento. Se presentará un gráfico comparativo que ilustrará cómo varían estas métricas entre los diferentes meses, proporcionando una visión clara del rendimiento del modelo en el set de testeo. Este análisis permitirá identificar posibles áreas de mejora y ajustar el modelo según sea necesario.

# Crear el data.frame combinado con las métricas Test del Modelo predictivo RF.
(metricas_test_RF <- bind_rows(
  marzo = metricas_test_marzo_RF,
  junio = metricas_test_junio_RF,
  septiembre = metricas_test_septiembre_RF,
  diciembre = metricas_test_diciembre_RF,
  .id = "Mes"
))

# A tibble: 12 × 4
   Mes        .metric .estimator .estimate
   <chr>      <chr>   <chr>          <dbl>
 1 marzo      rmse    standard       1.40 
 2 marzo      rsq     standard       0.944
 3 marzo      mae     standard       1.08 
 4 junio      rmse    standard       1.53 
 5 junio      rsq     standard       0.958
 6 junio      mae     standard       1.12 
 7 septiembre rmse    standard       1.22 
 8 septiembre rsq     standard       0.975
 9 septiembre mae     standard       0.936
10 diciembre  rmse    standard       1.64 
11 diciembre  rsq     standard       0.946
12 diciembre  mae     standard       1.30 

# Graficar manteniendo el orden de los modelos por mes
data_test_FR_long <- metricas_test_RF |> 
  pivot_longer(cols = .estimate, names_to = "Métrica", values_to = "Valor") 

# Test RF (RMSE y R² y MAE entre meses)
ggplot(metricas_test_RF, aes(x = Mes, y = .estimate, fill = .metric)) +
  geom_col(alpha = 0.7, na.rm = TRUE, position = "dodge") +  # Usar position = "dodge" para separar las barras
  ggtitle("Evaluacion Set de Muestreo RF -- (RMSE y R² y MAE entre meses)") +
  xlab("MESES") +
  ylab("Valor de la Métrica") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_fill_manual(name = "Métrica", values = c("rmse" = "steelblue", "rsq" = "orange", 'mae'= 'grey'))  # Colores por métrica

09. Evaluación de los resultados del modelo considerando: Remuestreo con 10 fold:

En esta sección, se llevó a cabo una evaluación exhaustiva del modelo de Random Forest (RF) mediante un remuestreo con 10 folds. Este enfoque permitió validar la robustez del modelo al dividir los datos en 10 subconjuntos, lo que facilitó la estimación del rendimiento del modelo en diversas configuraciones. Se combinarán las métricas de rendimiento, como RMSE y R², Promedios de los fold, para los meses de marzo, junio, septiembre y diciembre, y se presentará un gráfico que compara estos valores entre meses. Esta evaluación es esencial para determinar la estabilidad y la capacidad de generalización del modelo a través de diferentes particiones de datos, proporcionando una base sólida para la interpretación de los resultados.

# Crear el data.frame combinado con las métricas Remuestreo v= 10 del Modelo predictivo RF.
(metricas_vfold_cv_RF <- bind_rows(
  marzo = metricas_vfold_cv_marzo_RF,
  junio = metricas_vfold_cv_junio_RF,
  septiembre = metricas_vfold_cv_septiembre_RF,
  diciembre = metricas_vfold_cv_diciembre_RF,
  .id = "Mes"
))

# A tibble: 8 × 7
  Mes        .metric .estimator  mean     n std_err .config             
  <chr>      <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
1 marzo      rmse    standard   1.52     10  0.110  Preprocessor1_Model1
2 marzo      rsq     standard   0.937    10  0.0201 Preprocessor1_Model1
3 junio      rmse    standard   1.37     10  0.180  Preprocessor1_Model1
4 junio      rsq     standard   0.876    10  0.0973 Preprocessor1_Model1
5 septiembre rmse    standard   1.63     10  0.290  Preprocessor1_Model1
6 septiembre rsq     standard   0.893    10  0.0755 Preprocessor1_Model1
7 diciembre  rmse    standard   1.40     10  0.209  Preprocessor1_Model1
8 diciembre  rsq     standard   0.923    10  0.0445 Preprocessor1_Model1

# Graficar manteniendo el orden de los modelos por mes
data_vfold_cv_FR_long <- metricas_vfold_cv_RF |> 
  pivot_longer(cols = mean, names_to = "Métrica", values_to = "Valor") 

# Test RF (RMSE y R² entre meses)
ggplot(metricas_vfold_cv_RF, aes(x = Mes, y = mean, fill = .metric)) +
  geom_col(alpha = 0.7, na.rm = TRUE, position = "dodge") +  # Usar position = "dodge" para separar las barras
  ggtitle("Evaluacion Set de Muestreo RF, con 10 folds -- (RMSE y R² entre meses)") +
  xlab("MESES") +
  ylab("Valor de la Métrica") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_fill_manual(name = "Métrica", values = c("rmse" = "steelblue", "rsq" = "orange"))  # Colores por métrica

10. Comparación de los resultados RMSE y R² (Kriging Ordinario – Regresión-Kriging – Random Forest):

En esta sección, se realizó una comparación de métricas de rendimiento entre tres métodos de modelado: Kriging Ordinario (OK), Regresión-Kriging (RK) y Random Forest (RF). Se analizarón las métricas RMSE (Error Cuadrático Medio) y R² (Coeficiente de Determinación) para cada mes, lo que permitió evaluar la precisión y eficacia de cada enfoque. Los datos de rendimiento se combinarán y se visualizarán en un gráfico que permite observar cómo varían estas métricas entre los métodos y los meses de análisis. Esta comparación es crucial para identificar cuál de los métodos proporciona un mejor ajuste a los datos y una mayor capacidad predictiva en el contexto del estudio de la temperatura.

# Cargar las métricas totales por método
metricas_totales_OK <- read_rds('data5/procesada/metricas_totales_OK.rds')
metricas_totales_RK <- read_rds('data5/procesada/metricas_totales_RK.rds')
metricas_totales_RF <- read_rds('data5/procesada/metricas_totales_RF.rds')

# Convertir las métricas de RF a formato ancho
metricas_RF_wide <- metricas_totales_RF %>% 
  select(Mes, .metric, .estimate) %>% 
  pivot_wider(names_from = .metric, values_from = .estimate, 
              values_fill = NA) %>%
  rename(RMSE = rmse, R2 = rsq)

# Unir las métricas OK y RK en un solo dataframe
(metricas_combinadas <- bind_rows(
  mutate(metricas_totales_OK, Metodo = "OK"),
  mutate(metricas_totales_RK, Metodo = "RK"),
  mutate(metricas_RF_wide, Metodo = "RF") %>% 
    select(Mes, RMSE, R2, Metodo)  # Asegúrate de tener las mismas columnas
))

          Mes     RMSE        R2 Metodo
1       Marzo 1.503599 0.9172096     OK
2       Junio 1.789455 0.9488037     OK
3  Septiembre 1.742595 0.9587974     OK
4   Diciembre 1.397218 0.9596065     OK
5       Marzo 1.667656 0.8983197     RK
6       Junio 1.361883 0.9693096     RK
7  Septiembre 1.576472 0.9656508     RK
8   Diciembre 1.632773 0.9448681     RK
9       Marzo 1.402282 0.9456539     RF
10      Junio 1.477279 0.9609757     RF
11 Septiembre 1.272083 0.9719306     RF
12  Diciembre 1.608564 0.9499205     RF

# Convertir a formato largo para facilitar la visualización
metricas_long <- metricas_combinadas %>% 
  pivot_longer(cols = c(RMSE, R2), names_to = "Métrica", values_to = "Valor")

# Graficar las métricas (RMSE y R²) entre meses para cada método
ggplot(metricas_long, aes(x = Mes, y = Valor, fill = Métrica)) +
  geom_col(alpha = 0.7, na.rm = TRUE, position = "dodge") +  
  ggtitle("Comparación de Métricas (Kriging Ordinario -- Random Forest -- Regresión-Kriging)") +
  xlab("MESES") +
  ylab("Valor de la Métrica") +
  facet_wrap(~ Metodo) +  # Crear un panel por cada método (OK, RK, RF)
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_fill_manual(name = "Métrica", values = c("RMSE" = "steelblue", "R2" = "orange"))

11. Comparacion de los modelos predictivos (Kriging Ordinario – Random Forest – Regresión-Kriging) para los meses trabajados:

En esta sección, se presentarán los mapas de predicción de temperatura para la Región Metropolitana, generados a partir de tres métodos de modelado: Kriging Ordinario (OK), Random Forest (RF) y Regresión-Kriging (RK), para los meses de marzo, junio, septiembre y diciembre. Se han creado mapas para cada método, utilizando rásters de predicción correspondientes y aplicando una paleta de colores viridis en seis clases, que permite visualizar las variaciones de temperatura de manera clara. Los mapas incluyen elementos como bordes de la región, una rosa de los vientos y títulos que especifican el método y el mes analizado. Finalmente, se combinan los mapas en una visualización horizontal para cada mes, facilitando la comparación visual del desempeño de los diferentes modelos en la predicción de temperatura.

# Mostrar Mapas de Predicción de temperatura en la RM.

# Crear una lista para almacenar los mapas de predicción
mapas <- list()
# Definir los nombres de los meses
nombre_meses <- c("Marzo", "Junio", "Septiembre", "Diciembre")

# Cargar los mapas desde los archivos guardados
for (i in 1:4) {
  mes_nombre <- nombre_meses[i]
  
  # Cargar el ráster de predicción para cada método
  pred_raster_OK <- rast(paste0("data5/procesada/OK_", tolower(mes_nombre), "_predicciones_temp_promedio.tif"))
  pred_raster_RF <- rast(paste0("data5/procesada/RF_", tolower(mes_nombre), "_predicciones_temp_promedio.tif"))
  pred_raster_RK <- rast(paste0("data5/procesada/RK_", tolower(mes_nombre), "_predicciones_temp_promedio.tif"))
  
  # Crear mapas de predicción para cada método
  mapa_pred_OK <- tm_shape(pred_raster_OK[[1]]) +  # Usar la primera capa
    tm_raster(style = 'quantile', 
              palette = viridis(6),  # Cambiar a 6 clases
              title = "Pred. OK (°C)",  # Cambiar título de la leyenda
              colorNA = "white",  # Cambiar a blanco para los NA
              alpha = 0.9, 
              showNA = FALSE) +  # Eliminar "missing" de la leyenda
    tm_shape(region) +
    tm_borders(lwd = 2) +
    tm_compass(position = c("left", "top"), type = "4star", size = 1) +
    tm_layout(
      main.title = paste("Predicción T° OK -", mes_nombre),  
      main.title.size = 1,  
      main.title.position = "center",  
      main.title.fontface = "bold",  
      legend.outside = FALSE,  # Leyenda dentro del mapa
      legend.title.size = 0.9,  # Aumentar tamaño del título de la leyenda
      legend.text.size = 0.7,   # Aumentar tamaño del texto de la leyenda
      legend.position = c("right", "bottom"),  # Posición de la leyenda dentro del mapa
      frame = FALSE  
    )
  
  mapa_pred_RF <- tm_shape(pred_raster_RF) +  # Usar la única capa
    tm_raster(style = 'quantile', 
              palette = viridis(6),  # Cambiar a 6 clases
              title = "Pred. RF (°C)",  # Cambiar título de la leyenda
              colorNA = "white",  # Cambiar a blanco para los NA
              alpha = 0.9, 
              showNA = FALSE) +  # Eliminar "missing" de la leyenda
    tm_shape(region) +
    tm_borders(lwd = 2) +
    tm_compass(position = c("left", "top"), type = "4star", size = 1) +
    tm_layout(
      main.title = paste("Predicción T° RF -", mes_nombre),  
      main.title.size = 1,  
      main.title.position = "center",  
      main.title.fontface = "bold",  
      legend.outside = FALSE,  # Leyenda dentro del mapa
      legend.title.size = 0.9,  # Aumentar tamaño del título de la leyenda
      legend.text.size = 0.7,   # Aumentar tamaño del texto de la leyenda
      legend.position = c("right", "bottom"),  # Posición de la leyenda dentro del mapa
      frame = FALSE  
    )
  
  mapa_pred_RK <- tm_shape(pred_raster_RK[[1]]) +  # Usar la primera capa
    tm_raster(style = 'quantile', 
              palette = viridis(6),  # Cambiar a 6 clases
              title = "Pred. RK (°C)",  # Cambiar título de la leyenda
              colorNA = "white",  # Cambiar a blanco para los NA
              alpha = 0.9, 
              showNA = FALSE) +  # Eliminar "missing" de la leyenda
    tm_shape(region) +
    tm_borders(lwd = 2) +
    tm_compass(position = c("left", "top"), type = "4star", size = 1) +
    tm_layout(
      main.title = paste("Predicción T° RK -", mes_nombre),  
      main.title.size = 1,  
      main.title.position = "center",  
      main.title.fontface = "bold",  
      legend.outside = FALSE,  # Leyenda dentro del mapa
      legend.title.size = 0.9,  # Aumentar tamaño del título de la leyenda
      legend.text.size = 0.7,   # Aumentar tamaño del texto de la leyenda
      legend.position = c("right", "bottom"),  # Posición de la leyenda dentro del mapa
      frame = FALSE  
    )
  
  # Almacenar mapas en la lista
  mapas[[paste("OK -", mes_nombre)]] <- mapa_pred_OK
  mapas[[paste("RF -", mes_nombre)]] <- mapa_pred_RF
  mapas[[paste("RK -", mes_nombre)]] <- mapa_pred_RK
}

# Cambiar a modo de visualización
tmap_mode("plot")  

# Mostrar mapas en una sola ventana horizontal para cada mes
for (i in 1:4) {
  mes_nombre <- nombre_meses[i]
  
  # Combinar mapas en una sola ventana
  mapa_comb <- tmap_arrange(
    mapas[[paste("OK -", mes_nombre)]], 
    mapas[[paste("RF -", mes_nombre)]], 
    mapas[[paste("RK -", mes_nombre)]],
    ncol = 3  # Establecer el número de columnas
  )
  
  # Imprimir el mapa combinado
  print(mapa_comb)
}

Discución:

A lo largo del curso se utilizaron distintos métodos de análisis de datos espaciales con características y alcances diferentes. En un principio se efectuaron modelamientos con métodos no geoestadísticos de carácter determinístico, como los son el vecino más cercano y la interpolación de distancia inversa ponderada. Posteriormente, se aumentó el nivel de complejidad en el análisis incorporando el método geoestadístico del Kriging con sus variantes de Kriging Ordinario y Regresión-Kriging. Donde este último necesariamente requiere de un análisis variográfico para poder ser utilizado.

Finalmente, se alcanzó a utilizar la tecnología de los procesos de Machine Learning, utilizando específicamente el método de Random Forest que bajo su algoritmo de aprendizaje automático logra combinar varios modelos individuales y manejar muchas variables predictoras a la vez. Una de las grandes ventajas de este método en comparación con el análisis de datos espaciales geoestadísticos clásicos es que es posible evitarse el proceso completo del análisis variográfico.

Considerando los resultados entregados por las métricas de RMSE y R2 (parámetros que sirven como evaluadores de los distintos modelos) para cada uno de estas técnicas en los distintos meses trabajados se tiene que las técnicas no geoestadísticas mostraron los peores resultados en general, siendo el IDW la técnica con peor desempeño con un RMSE que se mueve entre 1,94 y 2,39. Mostrando el error más alto de todos los métodos considerados. Los valores de R2 fueron constantemente inferiores a los otros métodos, alcanzando un mínimo de 0,895 en marzo.

Para el caso del vecino más próximo, el segundo peor estimador, se observó un desempeño uniforme a lo largo de los meses trabajados con un RMSE entre 1,44 y 1,55. En cuanto al R2, a pesar de no ser la mejor opción para los meses de junio y septiembre, obtuvo valores bastante altos del orden de 0,966 y 0,968.

Entrando en los métodos geoestadísticos, tenemos al Kriging Ordinario que fue el método que mostró los mejores números para el mes de diciembre, con el RMSE más bajo (1,39) y el R2 mas alto (0,959). En general se posicionó como un buen método con valores de RMSE moderados (1,39 a 1,78) y altos valores de R2. Por su parte, el Regresión-Kriging fue la mejor alternativa para el mes de junio con un RMSE de 1,36 y un R2 de 0,969. Aunque para los otros meses sus errores tiende a ser mayores, destacando marzo con 1,66.

Para el caso de la técnica de Random Forest se puede afirmar con los datos obtenidos que es el método que más destacó obteniendo los mejores resultados para dos meses, marzo y septiembre. En el caso de marzo mostró un buen desempeño con un RMSE de 1,40 y un R2 de 0,945. Ambos valores muy superiores a todos los demás métodos. En el caso de Septiembre su RMSE resultó de 1,27 y el R2 de 0,971.

En resumen los métodos geoestadísticos continúan siendo una buena alternativa a tener en consideración a la hora de realizar análisis de datos espaciales. Con el método de Regresión-Kriging destacando especialmente en junio. Por su parte, los métodos no geoestadísticos no lograron destacar en ninguno de los modelos de los cuatro meses trabajados y el método de IDW obtuvo el peor rendimiento de todos, que según este taller no sería muy recomendable para la estimación de variables como la temperatura en áreas tan extensas como la Región Metropolitana.

En definitiva, el método de Random Forest se muestra como la alternativa más robusta para el análisis espacial y la estimación de datos mediante interpolación. Esto debido a que dentro de su análisis incorpora una gran cantidad de predictores, como fue el caso de este ejercicio donde se trabajó con cerca de 60 predictores diferentes para cada mes que fortalecieron sus estimaciones y redujeron la varianza en toda el área de estudio. Además posee una serie de procesos que van analizando y validando la información en pequeños paquetes de datos que proporciona información sobre la importancia de las diferentes características en la predicción de la variable objetivo.

l’esercizio è finito