Objetivo del Análisis

El objetivo es construir y comparar dos modelos distintos de aptitud climática para la caña de azúcar:

  1. Modelo de Rangos Óptimos: Un mapa booleano basado en agregados anuales de temperatura y precipitación.
  2. Modelo de Similaridad Climática: Un mapa análogo basado en la distancia climática a puntos de referencia en el Valle del Cauca, Colombia.

Configuración Inicial: Carga de Librerías

Realizaremos la carga de librerías iniciales que permita responder las preguntas del caso establecido

# Instalar

# install.packages(c("terra", "sf", "tidyverse", "rnaturalearth", "RColorBrewer", "patchwork"))

# Para manejo de datos espaciales (reemplaza a raster)

library(terra) 

# Para manejo de datos vectoriales

library(sf) 

library(tidyterra)

# Para manipulación (dplyr) y gráficos (ggplot2)

library(tidyverse) 

# Para obtener shapefiles globales

library(rnaturalearth) 

# Para paletas de colores

library(RColorBrewer)

# Para combinar gráficos de ggplot

library(patchwork)

# Configuración de R: evita la notación científica
options(scipen = 999)

Carga de Datos y Agregados

Iniciaremos con cargar la información de precipitación y temperatura promedio a 10m luego de intentar a 30s, sin exito, adicional procederemos a realizar el cargue de agregados anuales, por las capacidades de mi pc, se procede solo a cargar sudamerica, buscando tener un análisis correspondiente por territorio

path_base <- "C:/Maestria/Semestre 2/Analisis de Datos Geográficos y Espaciales/Modulo 2/Actividad"

# --- 2. Construir las rutas a las carpetas de datos ---

path_prec <- file.path(path_base, "wc2.1_10m_prec")
path_temp <- file.path(path_base, "wc2.1_10m_tavg")

# --- 3. Definir Extensión de Recorte (Sudamérica) ---

# (Longitud Min, Long Max, Lat Min, Lat Max)

extent_sa <- c(-90, -30, -60, 15)
print(paste("Cargando datos. TODO SE RECORTARÁ a la extensión:", extent_sa))
## [1] "Cargando datos. TODO SE RECORTARÁ a la extensión: -90"
## [2] "Cargando datos. TODO SE RECORTARÁ a la extensión: -30"
## [3] "Cargando datos. TODO SE RECORTARÁ a la extensión: -60"
## [4] "Cargando datos. TODO SE RECORTARÁ a la extensión: 15"
# --- 4. Carga de Precipitación (y recorte inmediato) ---

if (!dir.exists(path_prec)) stop("Error: No se encuentra el directorio: ", path_prec)
rutas_prec <- list.files(path_prec, pattern = "\\.tif$", full.names = TRUE)
if (length(rutas_prec) == 0) stop("Error: No se encontraron archivos .tif en: ", path_prec)

# Usamos 'win = extent_sa' para cargar SOLO esa ventana

prec_all_sa <- rast(rutas_prec, win = extent_sa)

# --- 5. Carga de Temperatura (y recorte inmediato) ---

if (!dir.exists(path_temp)) stop("Error: No se encuentra el directorio: ", path_temp)
rutas_temp <- list.files(path_temp, pattern = "\\.tif$", full.names = TRUE)
if (length(rutas_temp) == 0) stop("Error: No se encontraron archivos .tif en: ", path_temp)

# Usamos 'win = extent_sa'

temp_all_sa <- rast(rutas_temp, win = extent_sa)

print("Datos de WorldClim (30s) cargados y RECORTADOS a Sudamérica exitosamente.")
## [1] "Datos de WorldClim (30s) cargados y RECORTADOS a Sudamérica exitosamente."
# --- 6. Cálculo de Agregados Anuales (Sobre datos recortados) ---

# Estas operaciones ahora son ligeras

prec_anual_total_sa <- sum(prec_all_sa)
names(prec_anual_total_sa) <- "prec_anual"

tmean_anual_sa <- mean(temp_all_sa)
names(tmean_anual_sa) <- "tmean_anual"

print("Agregados anuales calculados (sobre Sudamérica).")
## [1] "Agregados anuales calculados (sobre Sudamérica)."

Mapa de Aptitud por Rangos Optimos (Sudamérica)

Construiremos mapas de Aptitud por rangos optimos para la caña de azucar, graficando los mismos, donde los clasificaremos en 2 rangos, Apto y no Apto de acuerdo a los parametros establecidos desde el caso solicitado de una precipitación anual de 1500 y 3500 milímetros y una temperatura media de 22.5 y 28 grados centígrados

# Rangos óptimos
RANGOS_OPTIMOS <- list(
  prec_anual_min = 1500, prec_anual_max = 3500,
  tmean_anual_min = 22.5, tmean_anual_max = 28
)

# ---- Lógica de Clasificación (sobre rasters _sa) ----

aptitud_prec <- (prec_anual_total_sa >= RANGOS_OPTIMOS$prec_anual_min) & 
                (prec_anual_total_sa <= RANGOS_OPTIMOS$prec_anual_max)

aptitud_tmean <- (tmean_anual_sa >= RANGOS_OPTIMOS$tmean_anual_min) & 
                 (tmean_anual_sa <= RANGOS_OPTIMOS$tmean_anual_max)

aptitud_total_sa <- aptitud_prec * aptitud_tmean
names(aptitud_total_sa) <- "Aptitud"

# ---- Visualización del Mapa (Sudamérica) ----

aptitud_df_sa <- as.data.frame(aptitud_total_sa, xy = TRUE, na.rm = TRUE)
aptitud_df_sa$Aptitud <- factor(aptitud_df_sa$Aptitud, 
                             levels = c(0, 1), 
                             labels = c("No Apto", "Apto"))

# Cargamos el mapa de Sudamérica para el fondo

mapa_sa <- ne_countries(continent = "South America", returnclass = "sf")

# Creamos el gráfico

g_aptitud_sa <- ggplot() +
  geom_raster(data = aptitud_df_sa, aes(x = x, y = y, fill = Aptitud)) +
  geom_sf(data = mapa_sa, fill = NA, color = "gray50", size = 0.2) +
  scale_fill_manual(values = c("No Apto" = "#d7191c", "Apto" = "#1a9641")) +
  # Usamos la extensión de Sudamérica para los límites
  coord_sf(xlim = extent_sa[1:2], ylim = extent_sa[3:4], expand = FALSE) +
  labs(
    title = "Modelo 1: Aptitud Climática en Sudamérica (Agregados Anuales)",
    subtitle = "Basado en rangos booleanos (1500-3500mm & 22.5-28°C)",
    x = "Longitud", y = "Latitud", fill = "Aptitud"
  ) +
  theme_bw()

print(g_aptitud_sa)

El área verde que domina el mapa, principalmente de la cuenca del amazonas, partes de Colombia, Ecuador, Perú y el norte de brasil, siendo así la mayor parte de sudamérica tropical es lo suficientemente cálida y lo suficientemente húmeda para soporte el cultivo de caña de azúcar.

Del modo contrario se puede observar dónde falla el modelo como lo es el cono sur conocido como Argentina, chile y uruguay, sin dejar atrás la cordillera de los andes.

Sin embargo este tipo de modelo, al ser anual está sobreestimando masivamente el área de aptitud real, al ignorar la estacionalidad.

Corte de Países de Alto Potenciales

Identificaremos 2 a 3 paises , que basado en el mapa 1 generado, nuestro foco inicial será brasil, colombia y perú.

# Países de interés (dentro de Sudamérica)

paises_interes_sa <- c("Brazil", "Colombia", "Peru")

# Filtramos el shapefile de Sudamérica

paises_shape_sa <- mapa_sa %>%
  filter(name %in% paises_interes_sa)

# ---- Recorte (Crop) y Máscara (Mask) ----

aptitud_recortada_sa <- crop(aptitud_total_sa, paises_shape_sa)
aptitud_mascara_sa <- mask(aptitud_recortada_sa, paises_shape_sa)

# ---- Visualización de Países ----

aptitud_paises_df_sa <- as.data.frame(aptitud_mascara_sa, xy = TRUE, na.rm = TRUE)

if(nrow(aptitud_paises_df_sa) > 0) {
  names(aptitud_paises_df_sa)[3] <- "Aptitud"
  aptitud_paises_df_sa <- aptitud_paises_df_sa %>%
    mutate(Aptitud = factor(Aptitud, levels = c(0, 1), labels = c("No Apto", "Apto")))
} else {
  print("ADVERTENCIA: No se encontraron píxeles de aptitud (0 o 1) dentro de los países seleccionados.")
}

# ---- GRÁFICO----

g_aptitud_paises_sa <- ggplot(data = paises_shape_sa) +
  
  geom_raster(data = aptitud_paises_df_sa, aes(x = x, y = y, fill = Aptitud)) +
  geom_sf(fill = NA, color = "black", size = 0.5) +
  
  scale_fill_manual(values = c("No Apto" = "#d7191c", "Apto" = "#1a9641")) +
  coord_sf(expand = FALSE) + 
  
  labs(
    
    title = "Aptitud (Modelo 1) en Paises Seleccionados",
    subtitle = paste(paises_interes_sa, collapse = ", "),
    x = "Longitud", y = "Latitud"
  ) +
  theme_bw() +
  facet_wrap(~name)

print(g_aptitud_paises_sa)

De esta manera delimitamos los limites geograficos de los 3 paises evaluados, dependiendo del objetivo del análisis podriamos separar las gráficas de cada pais individual para mejorar su visualización

Puntos de Muesta (Valle del Cauca)

Capturamos la evaluación de 3 puntos estráteticos del valle del cauca, como lo es palmira, cali y tuluá, con el objetivo de verificar el comportamiento de cada uno de ellos,donde evaluaremos la temperatura mensual y la precipitación del mismo

# Puntos (Longitud, Latitud)

puntos_valle <- data.frame(
  Nombre = c("Palmira", "Cali", "Tuluá"),
  lon = c(-76.30, -76.52, -76.19),
  lat = c(3.53, 3.45, 4.08)
)

# ---- LÍNEA CRÍTICA ----

puntos_vect <- vect(puntos_valle, geom=c("lon", "lat"), crs="EPSG:4326")

# ---- Extracción de Datos (del stack _sa) ----
# Combinamos los 12 meses de temp y 12 de prec en un stack de 24 capas
# Los objetos 'temp_all_sa' y 'prec_all_sa' deben existir del chunk 'load_data'

climate_stack_full_sa <- c(temp_all_sa, prec_all_sa)
names(climate_stack_full_sa) <- c(
  paste0("Tmean_", 1:12), 
  paste0("Prec_", 1:12)
)

# Extraemos los perfiles climáticos

extracted_data <- terra::extract(climate_stack_full_sa, puntos_vect)
extracted_df <- cbind(puntos_valle, extracted_data[, -1])

# ---- Preparación de Datos para Gráficos ----

data_long <- extracted_df %>%
  pivot_longer(
    cols = starts_with("Tmean_") | starts_with("Prec_"),
    names_to = "Variable", values_to = "Valor"
  ) %>%
  separate(Variable, into = c("Tipo", "Mes"), sep = "_") %>%
  mutate(Mes = as.numeric(Mes))

data_temp <- data_long %>% filter(Tipo == "Tmean")
data_prec <- data_long %>% filter(Tipo == "Prec")

# ---- Gráfico de Series de Tiempo (Temperatura) ----

g_temp_series <- ggplot(data_temp, aes(x = Mes, y = Valor, color = Nombre, group = Nombre)) +
  geom_line(size = 1) + geom_point(size = 2) +
  scale_x_continuous(breaks = 1:12, labels = month.abb) +
  labs(
    title = "Perfil de Temperatura Media Mensual",
    subtitle = "Puntos de Referencia en Valle del Cauca",
    x = "Mes", y = "Temperatura Media (°C)", color = "Sitio"
  ) + theme_minimal()

# ---- Gráfico de Series de Tiempo (Precipitación) ----

g_prec_series <- ggplot(data_prec, aes(x = Mes, y = Valor, fill = Nombre)) +
  geom_col(position = "dodge") + 
  scale_x_continuous(breaks = 1:12, labels = month.abb) +
  labs(
    # Quitamos tildes para evitar errores de codificación
    
    title = "Perfil de Precipitacion Mensual", 
    subtitle = "Puntos de Referencia en Valle del Cauca",
    x = "Mes", y = "Precipitacion Total (mm)", fill = "Sitio"
  ) + theme_minimal()

# Mostramos los gráficos

print(g_temp_series / g_prec_series)

Para el caso de la temperatura entendemos mejor la estacionalidad observando el pico más bajo en los meses de octubre, noviembre y diciembre, siendo su caso contrario para túlua el mes de Julio siendo el más alto, seguido de agosto para cali, palmira.

Mapa de Similaridad Climática

Buscaremos lugares que se comportan climáticamente como el valle del cauca en este caso de palmira, mes a mes, para no cesgarnos por un promedio anual que puede llegar a ocultar, temas releventas para el análisis requerido y nuestro principal objetivo, para ello lo haremos con la distancia euclidiana y resolver la pregunta, ¿ Dónde más en sudamérica puedo encontrar un clima casi idéntico para replicar ese éxito? esto asumiendo que en palmira crece perfectamente la caña

# Perfil de referencia (Palmira) - 24 variables

perfil_referencia <- extracted_df %>%
  filter(Nombre == "Palmira") %>%
  select(starts_with("Tmean_") | starts_with("Prec_")) %>%
  as.numeric()

# ---- Cálculo de Similaridad  ----

fn_dist_euclidiana <- function(v) {
  distancia <- sqrt(sum((v - perfil_referencia)^2, na.rm = TRUE))
  return(distancia)
}
print(paste("Iniciando cálculo de similaridad... (sobre Sudamérica)"))
## [1] "Iniciando cálculo de similaridad... (sobre Sudamérica)"
mapa_distancia_sa <- app(climate_stack_full_sa, fun = fn_dist_euclidiana)
print("Cálculo de similaridad terminado.")
## [1] "Cálculo de similaridad terminado."
mapa_similaridad_sa <- 1 / (1 + mapa_distancia_sa)
names(mapa_similaridad_sa) <- "Similaridad"

# ---- Visualización del Mapa ----

# Creamos el gráfico de ggplot

g_similaridad_sa <- ggplot() +
  
  # 1. Usamos geom_spatraster para dibujar el raster directamente
  
  geom_spatraster(data = mapa_similaridad_sa, aes(fill = Similaridad)) +
  
  # 2. Añadimos las fronteras (mapa_sa viene del Paso 1)
  
  geom_sf(data = mapa_sa, fill = NA, color = "gray50", size = 0.2) +
  
  # 3. Definimos la escala de color
  
  scale_fill_distiller(palette = "RdYlGn", direction = 1, na.value = NA) +
  
  coord_sf(xlim = extent_sa[1:2], ylim = extent_sa[3:4], expand = FALSE) +
  labs(
    title = "Modelo 2: Similaridad Climatica con Palmira (Sudamerica)", # Sin tilde
    subtitle = "Basado en Distancia Euclidiana (Enfoque Analogo)", # Sin tilde
    x = "Longitud", y = "Latitud",
    fill = "Similaridad\n(Más = Mejor)"
  ) +
  theme_bw()

print(g_similaridad_sa)

De acuerdo a la visualización podemos observar que el clima de palmira es extremadamente especifico y poco común

# Usamos 'patchwork' para poner los dos mapas (g_aptitud_sa y g_similaridad_sa)

g_comparacion_sa <- g_aptitud_sa + g_similaridad_sa + 
  plot_layout(guides = 'collect') # Alinea las leyendas

print(g_comparacion_sa)