Introducción

El presente proyecto desarrolla un análisis agroclimático espacial diseñado para evaluar el potencial global y regional del cultivo de la caña de azúcar mediante la contrastación de dos enfoques metodológicos complementarios. Inicialmente, se construye un modelo deductivo fundamentado en los rangos fisiológicos óptimos de temperatura y precipitación para mapear la aptitud climática mundial y delimitar países de alto rendimiento. Posteriormente, el estudio transita hacia un enfoque empírico de analogía climática, tomando como referente ecosistémico el Valle del Cauca; a partir de la extracción de series temporales en puntos específicos de esta zona —reconocida por su alta productividad—, se aplican métricas de similaridad (como la distancia euclidiana) para rastrear y proyectar regiones alrededor del mundo que repliquen estas condiciones exactas. La comparación analítica de ambos mapas —el teórico basado en umbrales y el análogo basado en similitud matemática— permite establecer conclusiones robustas sobre la viabilidad del cultivo, validando áreas estratégicas mediante el cruce de la teoría biológica con el comportamiento climático de un nicho exitoso comprobado.

# Instalar los paquetes necesarios si no están disponibles en el entorno
# install.packages(c("terra", "geodata", "ggplot2", "tidyterra"))

library(terra)
## Warning: package 'terra' was built under R version 4.5.2
## terra 1.8.93
library(geodata)
## Warning: package 'geodata' was built under R version 4.5.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.2
library(tidyterra)
## Warning: package 'tidyterra' was built under R version 4.5.3
## 
## Adjuntando el paquete: 'tidyterra'
## The following object is masked from 'package:stats':
## 
##     filter
# ---------------------------------------------------------
# 1. Descarga de datos de clima base (WorldClim) a nivel global
# ---------------------------------------------------------
# Usamos una resolución de 10 minutos de arco para un procesamiento global eficiente
cat("Descargando datos climáticos...\n")
## Descargando datos climáticos...
tavg <- worldclim_global(var = "tavg", res = 10, path = tempdir())
## Cached as: C:\Users\ACER\AppData\Local\Temp\Rtmp8WaS6M/climate/wc2.1_10m/wc2.1_10m_tavg.zip
prec <- worldclim_global(var = "prec", res = 10, path = tempdir())
## Cached as: C:\Users\ACER\AppData\Local\Temp\Rtmp8WaS6M/climate/wc2.1_10m/wc2.1_10m_prec.zip
# ---------------------------------------------------------
# 2. Procesamiento de variables bioclimáticas anuales
# ---------------------------------------------------------
# Calcular Temperatura Media Anual
tmean_annual <- mean(tavg)

# Calcular Precipitación Total Anual
prec_annual <- sum(prec)


# ---------------------------------------------------------
# 3. Definición de reglas de reclasificación (Modelo tipo EcoCrop)
# ---------------------------------------------------------
# Clases de aptitud: 0 = No Apto, 1 = Marginal/Moderado, 2 = Óptimo

# Matriz para Temperatura (°C)
m_temp <- c(-50, 15, 0,  # Menor a 15: No apto
            15,  24, 1,  # Entre 15 y 24: Marginal
            24,  30, 2,  # Entre 24 y 30: Óptimo
            30,  38, 1,  # Entre 30 y 38: Marginal
            38,  60, 0)  # Mayor a 38: No apto
rcl_temp <- matrix(m_temp, ncol=3, byrow=TRUE)
aptitud_temp <- classify(tmean_annual, rcl_temp)

# Matriz para Precipitación (mm)
m_prec <- c(-1,   1000,  0,  # Menor a 1000: No apto
            1000, 1500,  1,  # Entre 1000 y 1500: Marginal
            1500, 2500,  2,  # Entre 1500 y 2500: Óptimo
            2500, 3000,  1,  # Entre 2500 y 3000: Marginal
            3000, 10000, 0)  # Mayor a 3000: No apto
rcl_prec <- matrix(m_prec, ncol=3, byrow=TRUE)
aptitud_prec <- classify(prec_annual, rcl_prec)

# ---------------------------------------------------------
# 4. Cálculo de Aptitud Final
# ---------------------------------------------------------
# Se aplica el Principio de Factores Limitantes (Ley de Liebig):
# El factor más restrictivo determina la aptitud general.
aptitud_final <- min(aptitud_temp, aptitud_prec)

# Convertir a factor para facilitar la asignación de colores en ggplot
aptitud_final <- as.factor(aptitud_final)
levels(aptitud_final) <- data.frame(
  ID = c(0, 1, 2), 
  Aptitud = c("No Apto", "Aptitud Marginal", "Aptitud Óptima")
)

# ---------------------------------------------------------
# 5. Visualización del Mapa a Nivel Mundial
# ---------------------------------------------------------
ggplot() +
  geom_spatraster(data = aptitud_final) +
  scale_fill_manual(
    values = c(
      "No Apto" = "#E0E0E0",           # Gris claro
      "Aptitud Marginal" = "#FFC107",  # Amarillo/Ámbar
      "Aptitud Óptima" = "#4CAF50"     # Verde
    ),
    na.translate = FALSE # Evitar graficar zonas sin datos (ej. océanos)
  ) +
  labs(
    title = "Aptitud Climática Global para el Cultivo de Caña de Azúcar",
    subtitle = "Basado en rangos de Temperatura y Precipitación (FAO EcoCrop)",
    fill = "Nivel de Aptitud",
    x = "Longitud",
    y = "Latitud"
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold", size = 14),
    panel.grid.major = element_line(color = "gray90")
  )
## <SpatRaster> resampled to 5e+05 cells.

Aptitud Climática Global para el Cultivo de Caña de Azúcar

El modelo geoespacial evidencia que la aptitud climática para el cultivo de la caña de azúcar (Saccharum officinarum) está fuertemente delimitada por el cinturón tropical y subtropical del planeta, respondiendo a la alta demanda térmica e hídrica de esta especie (metabolismo C4). El mapa segmenta claramente el globo en tres macro-regiones:

Zonas de Aptitud Óptima (Verde): Se concentran de manera casi exclusiva en las latitudes ecuatoriales (aproximadamente entre los 15°N y 15°S). Destacan la cuenca del Amazonas y el Orinoco en Sudamérica, la cuenca del Congo en África Central y gran parte del Sudeste Asiático Insular (Indonesia, Malasia, Filipinas). Estas áreas garantizan de forma natural el régimen pluviométrico (1500 - 2500 mm anuales) y las temperaturas medias sostenidas (24 - 30°C) necesarias para maximizar el rendimiento de biomasa y acumulación de sacarosa sin necesidad de riego suplementario crítico.

Zonas de Aptitud Marginal o Moderada (Amarillo): Actúan como franjas de transición hacia los subtrópicos (hasta los 30°-35° de latitud N y S) e incluyen regiones históricamente productoras como el centro-sur de Brasil, la India, el sureste de China, el este de Australia, Centroamérica y el sureste de Estados Unidos. En estas zonas, el cultivo es comercialmente viable, pero se ve limitado estacionalmente por uno de los dos factores: inviernos más frescos que ralentizan el crecimiento, o precipitaciones insuficientes que obligan a la implementación de sistemas de riego para alcanzar niveles óptimos de producción.

Zonas No Aptas (Gris): Comprenden las latitudes medias y altas, así como los grandes biomas áridos (Sahara, Medio Oriente, centro de Australia). Las limitantes aquí son severas y excluyentes: temperaturas gélidas o heladas (la caña es altamente susceptible al daño por congelación) o un déficit hídrico extremo que hace agronómica y económicamente inviable el cultivo a gran escala bajo condiciones de secano.

Implicación Agronómica: Este mapa de línea base corrobora que, bajo condiciones estrictamente de secano (sin riego) y analizando promedios anuales, la expansión natural y más rentable de la caña de azúcar compite directamente con los ecosistemas de bosque tropical húmedo. Por lo tanto, en las zonas verdes, el desafío principal no es climático, sino de sostenibilidad ambiental y uso del suelo; mientras que en las zonas amarillas, el reto es tecnológico (riego y manejo agronómico).

Top 3 paises con mejor rendimiento climático para la producción de caña de azucar a nivel mundial

cat("Descargando límites administrativos (shapefiles) de los países...\n")
## Descargando límites administrativos (shapefiles) de los países...
# Usamos gadm() de la librería geodata para obtener las fronteras (nivel 0 = país completo)
bra <- gadm(country = "BRA", level = 0, path = tempdir())
## Cached as: C:\Users\ACER\AppData\Local\Temp\Rtmp8WaS6M/gadm/gadm41_BRA_0_pk.rds
ind <- gadm(country = "IND", level = 0, path = tempdir())
## Cached as: C:\Users\ACER\AppData\Local\Temp\Rtmp8WaS6M/gadm/gadm41_IND_0_pk.rds
col <- gadm(country = "COL", level = 0, path = tempdir())
## Cached as: C:\Users\ACER\AppData\Local\Temp\Rtmp8WaS6M/gadm/gadm41_COL_0_pk.rds
# Unir los tres polígonos en un solo vector espacial
paises_top3 <- rbind(bra, ind, col)

# ---------------------------------------------------------
# 1. Corte (Crop) y Enmascarado (Mask)
# ---------------------------------------------------------
cat("Realizando el corte espacial...\n")
## Realizando el corte espacial...
# Crop reduce el marco delimitador (bounding box) al área de los países
aptitud_crop <- crop(aptitud_final, paises_top3)

# Mask elimina los valores que caen fuera de las fronteras exactas (ej. océano o países vecinos)
aptitud_top3 <- mask(aptitud_crop, paises_top3)

# ---------------------------------------------------------
# 2. Visualización Específica
# ---------------------------------------------------------
ggplot() +
  # Capa raster recortada
  geom_spatraster(data = aptitud_top3) +
  # Superponer las fronteras de los países en negro
  geom_spatvector(data = paises_top3, fill = NA, color = "black", linewidth = 0.4) +
  scale_fill_manual(
    values = c(
      "No Apto" = "#E0E0E0",
      "Aptitud Marginal" = "#FFC107",
      "Aptitud Óptima" = "#4CAF50"
    ),
    na.translate = FALSE 
  ) +
  labs(
    title = "Aptitud Climática para Caña de Azúcar en Países Clave",
    subtitle = "Brasil, India y Colombia (Corte de Línea Base Global)",
    fill = "Nivel de Aptitud",
    x = "Longitud",
    y = "Latitud"
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold", size = 14),
    panel.grid.major = element_line(color = "gray90")
  )

Para este análisis, se seleccionaron Brasil, India y Colombia como casos de estudio estratégicos. Brasil e India fueron elegidos por ser los mayores productores mundiales de biomasa de caña de azúcar, concentrando las extensiones territoriales de aptitud climática óptima más grandes del planeta. Por su parte, Colombia se incluyó debido a que, a pesar de contar con una menor extensión total cultivada, posee nichos ecológicos específicos —como el Valle del Cauca— que registran algunos de los rendimientos agroindustriales por hectárea más altos a nivel global, lo cual resulta ideal para establecer un modelo empírico de alta eficiencia en los posteriores análisis de similaridad climática.

# Instalar patchwork si no lo tienes: install.packages("patchwork")
library(terra)
library(geodata)
library(ggplot2)
library(tidyterra)
library(patchwork) # Excelente para unir gráficos de ggplot
## 
## Adjuntando el paquete: 'patchwork'
## The following object is masked from 'package:terra':
## 
##     area
# Asumimos que los objetos 'tmean_annual', 'prec_annual', 'bra', 'ind' y 'col' 
# ya están cargados en tu entorno a partir de los pasos anteriores.

# 1. Crear una lista con los polígonos de los tres países
lista_paises <- list(
  "Brasil" = bra, 
  "India" = ind, 
  "Colombia" = col
)

# 2. Bucle para generar y graficar los mapas por país
for(nombre_pais in names(lista_paises)) {
  
  # Extraer el polígono del país actual
  poligono <- lista_paises[[nombre_pais]]
  
  cat("Procesando mapas para:", nombre_pais, "\n")
  
  # --- CORTE Y ENMASCARADO ---
  # Temperatura
  temp_pais <- mask(crop(tmean_annual, poligono), poligono)
  # Precipitación
  prec_pais <- mask(crop(prec_annual, poligono), poligono)
  
  # --- GRÁFICO DE TEMPERATURA ---
  # Usamos la paleta "inferno" (negro-morado-naranja-amarillo) ideal para calor
  mapa_temp <- ggplot() +
    geom_spatraster(data = temp_pais) +
    geom_spatvector(data = poligono, fill = NA, color = "black", linewidth = 0.3) +
    scale_fill_viridis_c(
      option = "inferno", 
      name = "Temp (°C)",
      na.value = "transparent"
    ) +
    labs(
      title = "Temperatura Media Anual",
      x = "Longitud", y = "Latitud"
    ) +
    theme_minimal() +
    theme(legend.position = "right", plot.title = element_text(size = 11, face = "bold"))
  
  # --- GRÁFICO DE PRECIPITACIÓN ---
  # Usamos la paleta "mako" invertida (azules oscuros a claros) ideal para agua
  mapa_prec <- ggplot() +
    geom_spatraster(data = prec_pais) +
    geom_spatvector(data = poligono, fill = NA, color = "black", linewidth = 0.3) +
    scale_fill_viridis_c(
      option = "mako", 
      direction = -1, # Invertir para que más lluvia sea más oscuro/azul
      name = "Prec (mm)",
      na.value = "transparent"
    ) +
    labs(
      title = "Precipitación Total Anual",
      x = "Longitud", y = "Latitud"
    ) +
    theme_minimal() +
    theme(legend.position = "right", plot.title = element_text(size = 11, face = "bold"))
  
  # --- UNIÓN Y VISUALIZACIÓN ---
  # Usamos patchwork para ponerlos lado a lado
  grafico_combinado <- mapa_temp + mapa_prec + 
    plot_annotation(
      title = paste("Climatología de Línea Base (1970-2000) -", nombre_pais),
      theme = theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.4))
    )
  
  # Imprimir el gráfico resultante en la ventana de Plots
  print(grafico_combinado)
}
## Procesando mapas para: Brasil

## Procesando mapas para: India

## Procesando mapas para: Colombia

## 1. Brasil: El gigante de las planicies y la transición climática La climatología de Brasil evidencia por qué es el líder mundial en producción, pero también revela sus limitantes geográficas.

Análisis Térmico: Gran parte del territorio (norte, centro y costa nororiental) mantiene temperaturas óptimas (superiores a 24°C). Sin embargo, el mapa muestra un claro enfriamiento hacia el extremo sur (región Sur y partes de São Paulo/Paraná), donde las temperaturas medias caen por debajo de los 20°C, marcando la frontera térmica del cultivo debido al riesgo de heladas invernales que dañan los tejidos de la planta.

Análisis Hídrico: La cuenca amazónica (noroeste) presenta precipitaciones excesivas (superiores a 3000 mm), lo que agronómicamente reduce la acumulación de sacarosa por falta de estrés hídrico previo a la cosecha, además de representar un área de conservación. Por el contrario, el interior del noreste (región de la Caatinga) sufre de un déficit hídrico severo (<1000 mm).

Conclusión: La “zona núcleo” agroindustrial se consolida en el Centro-Sur (estado de São Paulo y alrededores) y la franja costera del noreste, donde se cruzan temperaturas cálidas moderadas con un régimen de lluvias suficiente (1500 - 2000 mm) y bien distribuido, permitiendo un secano altamente productivo.

2. India: El contraste topográfico y la dependencia monzónica

El subcontinente indio muestra la mayor variabilidad extrema de los tres países, fuertemente condicionada por el Himalaya y el régimen de los monzones.

Análisis Térmico: A excepción de la extrema cordillera norte (donde las temperaturas medias caen bajo 0°C, imposibilitando cualquier agricultura tropical), la inmensa mayoría de la península india y las llanuras del Ganges ofrecen un perfil térmico ideal y sostenido (tonos cálidos en el mapa, 25°C+).

Análisis Hídrico: El mapa de precipitaciones muestra extremos dramáticos. La costa oeste (Ghats occidentales) y el noreste reciben lluvias torrenciales (más de 4000-6000 mm), lo cual genera anegamiento, pudrición de raíces y baja luminosidad. En contraste, el noroeste (Rajastán) es desértico.

Conclusión: Las áreas de mayor potencial, como las llanuras del norte (Uttar Pradesh) y el centro-oeste (Maharashtra), gozan de excelente temperatura pero presentan precipitaciones moderadas a bajas en la línea base anual. Esto indica que, a diferencia de Brasil, el éxito masivo de la caña en la India depende de una fuerte infraestructura de riego para compensar la estacionalidad del monzón y sostener el cultivo durante los meses secos.

3. Colombia: El control orográfico y los microclimas de alto rendimiento

La distribución climática en Colombia no responde a la latitud, sino casi exclusivamente a la altitud impuesta por las tres cordilleras de los Andes.

Análisis Térmico: El mapa es una calca de la topografía. Las zonas andinas (tonos oscuros/morados) tienen temperaturas medias por debajo de los 15°C, descartándolas totalmente para la caña azucarera comercial (aunque aptas para la producción tradicional de panela a media ladera). Las grandes extensiones planas (Llanos orientales, costa Caribe y valles interandinos) ofrecen el calor necesario (amarillo intenso, ~24-28°C).

Análisis Hídrico: Colombia posee áreas de precipitación extrema, como la costa Pacífica (Chocó) y el piedemonte amazónico (más de 6000 mm), que superan ampliamente la tolerancia de la caña. El norte extremo (La Guajira) es demasiado árido.

Conclusión: La aptitud óptima en Colombia se encuentra encapsulada en microclimas específicos. El Valle del río Cauca (el cual evaluarás en tu siguiente paso) se destaca como un nicho geográfico perfecto: al estar encañonado entre cordilleras, está protegido de las lluvias extremas del Pacífico y la Amazonía, manteniendo un régimen pluviométrico moderado-adecuado con una temperatura cálida constante y alta radiación solar. Esto es lo que permite que en esta región se pueda sembrar y cosechar caña los 365 días del año, un fenómeno agroindustrial único en el mundo.

El Periodo Base: 1970 - 2000

Los datos extraídos provienen de la base global WorldClim (versión 2.1). Esta plataforma utiliza como estándar meteorológico un periodo ininterrumpido de 30 años, abarcando exactamente desde el 1 de enero de 1970 hasta el 31 de diciembre de 2000.

¿Por qué 30 años? La Organización Meteorológica Mundial (OMM) establece que para definir el “clima” de una región (y no solo su “estado del tiempo” o clima a corto plazo), se requiere un promedio mínimo de tres décadas. Esto se conoce como una Normal Climatológica.

Filtrado de Anomalías: Al promediar el clima entre 1970 y 2000, el modelo absorbe y suaviza eventos climáticos extremos. Si tomáramos solo un año reciente (por ejemplo, un año bajo un fenómeno de El Niño severo), las gráficas mostrarían una falsa sequía que no representa la realidad constante del Valle del Cauca.

¿Cómo se refleja este periodo en la gráfica? Cada marcador en tus gráficas es un concentrado matemático de esos 30 años:

En la Temperatura: El punto correspondiente al mes de “Ene” (Enero) no muestra el calor que hizo en enero de 1990 o de 1998. Muestra la temperatura media calculada al sumar todos los eneros de esos 30 años y dividirlos entre 30.

En la Precipitación: La barra del mes de “Abr” (Abril) representa el volumen de agua promedio que históricamente cae en ese mes, calculado a partir de los 30 abriles contenidos entre 1970 y el año 2000.

library(terra)
library(geodata)
library(ggplot2)
library(tidyterra)
library(patchwork)

# ---------------------------------------------------------
# 1. Descarga de datos de ALTA RESOLUCIÓN (0.5 minutos, ~1 km)
# ---------------------------------------------------------
cat("Descargando climatología de alta resolución para Colombia...\n")
## Descargando climatología de alta resolución para Colombia...
# Al descargar por país a máxima resolución, evitamos colapsar la memoria
tavg_col <- worldclim_country(country = "COL", var = "tavg", path = tempdir(), res = 0.5)
## Cached as: C:\Users\ACER\AppData\Local\Temp\Rtmp8WaS6M/climate/wc2.1_country/COL_wc2.1_30s_tavg.tif
prec_col <- worldclim_country(country = "COL", var = "prec", path = tempdir(), res = 0.5)
## Cached as: C:\Users\ACER\AppData\Local\Temp\Rtmp8WaS6M/climate/wc2.1_country/COL_wc2.1_30s_prec.tif
# Calculamos las variables anuales (Media de temp, Suma de lluvia)
tmean_col_annual <- mean(tavg_col)
prec_col_annual <- sum(prec_col)

# ---------------------------------------------------------
# 2. Definición del Polígono y Corte Espacial
# ---------------------------------------------------------
cat("Descargando límites y definiendo el polígono...\n")
## Descargando límites y definiendo el polígono...
# SOLUCIÓN: Definir valle_cauca en este mismo entorno antes de usarlo
col_dept <- gadm(country = "COL", level = 1, path = tempdir())
## Cached as: C:\Users\ACER\AppData\Local\Temp\Rtmp8WaS6M/gadm/gadm41_COL_1_pk.rds
valle_cauca <- col_dept[col_dept$NAME_1 == "Valle del Cauca", ]

cat("Realizando corte de alta definición...\n")
## Realizando corte de alta definición...
temp_valle_hr <- mask(crop(tmean_col_annual, valle_cauca), valle_cauca)
prec_valle_hr <- mask(crop(prec_col_annual, valle_cauca), valle_cauca)

# ---------------------------------------------------------
# 3. Visualización Cartográfica de Alta Resolución
# ---------------------------------------------------------
# Mapa de Temperatura
mapa_temp_hr <- ggplot() +
  geom_spatraster(data = temp_valle_hr) +
  geom_spatvector(data = valle_cauca, fill = NA, color = "black", linewidth = 0.5) +
  scale_fill_viridis_c(
    option = "inferno", 
    name = "Temp (°C)",
    na.value = "transparent"
  ) +
  labs(
    title = "Temperatura Media Anual",
    subtitle = "Valle del Cauca (Resolución ~1 km²)",
    x = "Longitud", y = "Latitud"
  ) +
  theme_minimal() +
  theme(legend.position = "right", plot.title = element_text(face = "bold"))

# Mapa de Precipitación
mapa_prec_hr <- ggplot() +
  geom_spatraster(data = prec_valle_hr) +
  geom_spatvector(data = valle_cauca, fill = NA, color = "black", linewidth = 0.5) +
  scale_fill_viridis_c(
    option = "mako", 
    direction = -1,
    name = "Prec (mm)",
    na.value = "transparent"
  ) +
  labs(
    title = "Precipitación Total Anual",
    subtitle = "Valle del Cauca (Resolución ~1 km²)",
    x = "Longitud", y = "Latitud"
  ) +
  theme_minimal() +
  theme(legend.position = "right", plot.title = element_text(face = "bold"))

# Unir y mostrar gráficos
grafico_hr_combinado <- mapa_temp_hr + mapa_prec_hr +
  plot_annotation(
    title = "Climatología de Alta Resolución - Valle del Cauca",
    theme = theme(plot.title = element_text(size = 16, face = "bold", hjust = 0.5))
  )

print(grafico_hr_combinado)

# ---------------------------------------------------------
# INSTALACIÓN Y CARGA DE LIBRERÍAS
# ---------------------------------------------------------
# install.packages(c("terra", "geodata", "ggplot2", "tidyr", "dplyr", "patchwork"))

library(terra)
library(geodata)
library(ggplot2)
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.5.2
## 
## Adjuntando el paquete: 'tidyr'
## The following object is masked from 'package:terra':
## 
##     extract
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.2
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:terra':
## 
##     intersect, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(patchwork)

# ---------------------------------------------------------
# 1. DESCARGA DE DATOS BASE (Si ya los tienes en memoria, esto será rápido)
# ---------------------------------------------------------
cat("Descargando datos climáticos mensuales globales...\n")
## Descargando datos climáticos mensuales globales...
tavg <- worldclim_global(var = "tavg", res = 10, path = tempdir())
prec <- worldclim_global(var = "prec", res = 10, path = tempdir())

cat("Descargando límites departamentales de Colombia...\n")
## Descargando límites departamentales de Colombia...
col_dept <- gadm(country = "COL", level = 1, path = tempdir())

# ---------------------------------------------------------
# 2. SELECCIÓN DE PUNTOS EN EL VALLE DEL CAUCA
# ---------------------------------------------------------
# Filtrar específicamente el Valle del Cauca
valle_cauca <- col_dept[col_dept$NAME_1 == "Valle del Cauca", ]

# Generar 3 puntos al azar (devuelve SpatVector directamente)
set.seed(42) 
puntos_vec <- spatSample(valle_cauca, size = 3, method = "random")

# Extraer las coordenadas exactas para los nombres
puntos_xy <- as.data.frame(crds(puntos_vec))

# ---------------------------------------------------------
# 3. EXTRACCIÓN DE SERIES DE TIEMPO (CORRECCIÓN DE CONFLICTO)
# ---------------------------------------------------------
cat("Extrayendo datos de los puntos (forzando uso de terra::extract)...\n")
## Extrayendo datos de los puntos (forzando uso de terra::extract)...
# El uso de terra:: evita el conflicto con la librería tidyr
temp_ext <- terra::extract(tavg, puntos_vec)
prec_ext <- terra::extract(prec, puntos_vec)

# ---------------------------------------------------------
# 4. LIMPIEZA Y TRANSFORMACIÓN DE DATOS
# ---------------------------------------------------------
procesar_serie <- function(df_extraido, variable_nombre, coordenadas) {
  # Quitar la columna ID generada por extract
  df_meses <- df_extraido[, -1] 
  
  # Nombres de meses
  meses_esp <- c("Ene", "Feb", "Mar", "Abr", "May", "Jun", 
                 "Jul", "Ago", "Sep", "Oct", "Nov", "Dic")
  colnames(df_meses) <- meses_esp
  
  # Asignar nombres a los puntos basados en sus coordenadas
  df_meses$Sitio <- paste0("Punto ", 1:nrow(df_meses), " (Lat: ", 
                           round(coordenadas$y, 2), ", Lon: ", 
                           round(coordenadas$x, 2), ")")
  
  # Transformar de formato ancho a largo
  df_largo <- df_meses %>%
    pivot_longer(cols = all_of(meses_esp), 
                 names_to = "Mes", 
                 values_to = "Valor") %>%
    mutate(
      Mes = factor(Mes, levels = meses_esp),
      Variable = variable_nombre
    )
  
  return(df_largo)
}

df_temp <- procesar_serie(temp_ext, "Temperatura", puntos_xy)
df_prec <- procesar_serie(prec_ext, "Precipitación", puntos_xy)

# ---------------------------------------------------------
# 5. VISUALIZACIÓN DE TENDENCIA (GGPLOT + PATCHWORK)
# ---------------------------------------------------------
cat("Generando gráficos...\n")
## Generando gráficos...
# Gráfico de Temperatura
grafico_temp <- ggplot(df_temp, aes(x = Mes, y = Valor, color = Sitio, group = Sitio)) +
  geom_line(linewidth = 1.2, alpha = 0.8) +
  geom_point(size = 3.5, shape = 21, fill = "white", stroke = 1.5) +
  scale_color_viridis_d(option = "plasma", end = 0.8) +
  labs(title = "Evolución Térmica Mensual",
       subtitle = "Temperatura media (°C) en el Valle del Cauca",
       y = "Temperatura (°C)", x = NULL) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "none",
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    plot.title = element_text(face = "bold", size = 14)
  )

# Gráfico de Precipitación
grafico_prec <- ggplot(df_prec, aes(x = Mes, y = Valor, fill = Sitio)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7, alpha = 0.9) +
  scale_fill_viridis_d(option = "mako", begin = 0.2, end = 0.8) +
  labs(title = "Régimen Pluviométrico Mensual",
       subtitle = "Precipitación acumulada (mm)",
       y = "Lluvia (mm)", x = "Mes del Año") +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "bottom",
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    plot.title = element_text(face = "bold", size = 14),
    legend.title = element_text(face = "bold"),
    legend.text = element_text(size = 9)
  ) +
  guides(fill = guide_legend(title = "Ubicación del Muestreo", ncol = 1))

# Unir gráficos
panel_series <- grafico_temp / grafico_prec + 
  plot_layout(heights = c(1, 1.2))

# Mostrar gráfico
print(panel_series)

1. Comportamiento de la Temperatura (Evolución Térmica)

Isotermia Ecuatorial (Estabilidad): Las tres líneas de temperatura son prácticamente planas a lo largo de los 12 meses. Esto demuestra la ausencia de estaciones térmicas (verano/invierno) típica de las latitudes ecuatoriales; la variación mes a mes en un mismo punto es casi nula (menos de 1°C de diferencia a lo largo del año).

Fuerte Gradiente Altitudinal: Aunque la temperatura no cambia en el tiempo, sí cambia drásticamente en el espacio. Hay una brecha de aproximadamente 5°C a 6°C entre el punto más cálido (azul oscuro, ~25.5°C) y el más fresco (naranja, ~20°C). Al ser puntos aleatorios dentro de un mismo departamento, esta diferencia confirma que cayeron en distintas elevaciones topográficas (el más cálido en la parte plana del valle y el más fresco hacia la ladera de la cordillera).

2. Comportamiento de la Precipitación (Régimen Pluviométrico)

Régimen Bimodal: Las gráficas de barras evidencian un patrón clásico de la región andina colombiana con dos temporadas de altas precipitaciones (abril-mayo y septiembre-noviembre) y dos temporadas de “secano” relativo o menores lluvias (enero-febrero y julio-agosto).

Alta Variabilidad Microclimática: A pesar de estar en el mismo departamento, los volúmenes de lluvia son drásticamente distintos. El sitio representado por las barras color morado oscuro es extremadamente húmedo (con picos que superan los 400 mm y 500 mm mensuales), mientras que el sitio representado por las barras verde claro se mantiene sistemáticamente por debajo de los 250 mm mensuales.

3. Síntesis Agronómica (Enfoque en Caña de Azúcar)

Contraste de Aptitud: El punto que registra la temperatura más alta (línea azul oscuro, ~25°C) junto con un régimen de precipitación moderado (probablemente las barras verde claro o azul medio, con lluvias entre 100 y 250 mm por mes) representa el nicho ideal para la caña de azúcar, permitiendo crecimiento continuo todo el año. Por el contrario, el punto naranja (20°C) ralentizaría el crecimiento, y el punto de lluvias extremas (morado oscuro, >500 mm/mes) generaría anegamiento y problemas severos de maduración y concentración de sacarosa.

# Asegúrate de tener cargadas: terra, ggplot2, tidyterra
library(terra)
library(ggplot2)
library(tidyterra)

# ---------------------------------------------------------
# 1. Definir el "Clima Ideal" del Valle del Cauca
# ---------------------------------------------------------
cat("Extrayendo perfil climático de los 3 puntos de referencia...\n")
## Extrayendo perfil climático de los 3 puntos de referencia...
# Usamos las coordenadas generadas en el paso anterior (puntos_vec)
ref_temp <- terra::extract(tmean_annual, puntos_vec)[,2]
ref_prec <- terra::extract(prec_annual, puntos_vec)[,2]

# Promediamos los 3 puntos para crear nuestro "Nicho de Referencia"
mean_ref_temp <- mean(ref_temp, na.rm = TRUE)
mean_ref_prec <- mean(ref_prec, na.rm = TRUE)

# ---------------------------------------------------------
# 2. Estandarización Estadística (Z-score)
# ---------------------------------------------------------
cat("Estandarizando variables climáticas a nivel global...\n")
## Estandarizando variables climáticas a nivel global...
# Calculamos medias y desviaciones estándar globales
global_temp_mean <- global(tmean_annual, "mean", na.rm=TRUE)[[1]]
global_temp_sd <- global(tmean_annual, "sd", na.rm=TRUE)[[1]]

global_prec_mean <- global(prec_annual, "mean", na.rm=TRUE)[[1]]
global_prec_sd <- global(prec_annual, "sd", na.rm=TRUE)[[1]]

# Aplicamos la fórmula Z = (X - Media) / Desviación Estándar
tmean_z <- (tmean_annual - global_temp_mean) / global_temp_sd
prec_z <- (prec_annual - global_prec_mean) / global_prec_sd

# Estandarizamos también nuestro punto de referencia
ref_temp_z <- (mean_ref_temp - global_temp_mean) / global_temp_sd
ref_prec_z <- (mean_ref_prec - global_prec_mean) / global_prec_sd

# ---------------------------------------------------------
# 3. Cálculo de Distancia Euclidiana y Similaridad
# ---------------------------------------------------------
cat("Calculando métrica de similaridad espacial...\n")
## Calculando métrica de similaridad espacial...
# Distancia Euclidiana: sqrt( (X_raster - X_ref)^2 + (Y_raster - Y_ref)^2 )
distancia_euclidiana <- sqrt((tmean_z - ref_temp_z)^2 + (prec_z - ref_prec_z)^2)

# Convertimos la Distancia a un Índice de Similaridad (0 a 1)
# Donde 1 es un clima matemáticamente idéntico al Valle del Cauca
max_dist <- global(distancia_euclidiana, "max", na.rm=TRUE)[[1]]
similaridad_global <- 1 - (distancia_euclidiana / max_dist)

# ---------------------------------------------------------
# 4. Visualización del Mapa Cartográfico
# ---------------------------------------------------------
cat("Generando mapa de similaridad...\n")
## Generando mapa de similaridad...
mapa_similaridad <- ggplot() +
  geom_spatraster(data = similaridad_global) +
  # Misma escala de colores: Gris (Diferente), Amarillo (Análogo Medio), Verde (Análogo Exacto)
  scale_fill_gradientn(
    colors = c("#E0E0E0", "#FFC107", "#4CAF50"), 
    name = "Índice de\nSimilaridad",
    na.value = "transparent"
  ) +
  labs(
    title = "Mapa Global de Analogía Climática",
    subtitle = "Zonas climáticamente análogas al clúster muestreado en el Valle del Cauca",
    x = "Longitud",
    y = "Latitud"
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold", size = 14),
    panel.grid.major = element_line(color = "gray90")
  )
## <SpatRaster> resampled to 5e+05 cells.
# Imprimir gráfico
print(mapa_similaridad)

Introducción a la Técnica: Analogía Climática mediante Distancia Euclidiana

El mapa presentado es el resultado de aplicar una técnica espacial empírica conocida como Analogía Climática, utilizando la Distancia Euclidiana como métrica estadística. A diferencia de los modelos deductivos (como el enfoque de la FAO) que evalúan si un territorio cumple con rangos teóricos de tolerancia (máximos y mínimos), esta técnica busca “clonar” un nicho ecológico específico que ya ha demostrado ser exitoso.

El proceso matemático consiste en:

Extraer la “firma climática” (los valores exactos de temperatura y precipitación) de tu sitio de referencia (el Valle del Cauca).

Estandarizar todas las variables climáticas a nivel global mediante el cálculo de Z-scores. Esto es vital para que la lluvia (medida en miles de milímetros) no opaque a la temperatura (medida en decenas de grados) en el cálculo.

Calcular la distancia geométrica en un espacio multivariado (la fórmula de Pitágoras extendida) entre el perfil del Valle del Cauca y cada píxel individual del planeta.

Invertir esta distancia para generar un índice de similaridad continuo, donde el valor máximo representa una coincidencia climática matemáticamente idéntica.

Análisis Breve de la Gráfica Resultante Al interpretar la visualización generada, se extraen las siguientes conclusiones clave sobre el comportamiento espacial del modelo:

Naturaleza Continua del Gradiente:

A diferencia del primer mapa de la FAO, que clasificaba el mundo en categorías rígidas (Apto, Marginal, No Apto), este gráfico proyecta un gradiente continuo. La inmensa capa amarilla que cubre gran parte del globo no significa “aptitud”, sino que representa el valor base o de “alta distancia matemática” (baja similitud). Ecosistemas boreales, templados o desérticos caen en este espectro amarillo porque sus perfiles estandarizados están muy alejados del cálido y húmedo trópico andino.

Alta Especificidad Microclimática (Zonas Verdes):

Las áreas que se tornan verdes (alta similitud) están drásticamente más restringidas que en el modelo de la FAO. Se limitan a franjas muy específicas y profundas del cinturón ecuatorial: núcleos en la cuenca Amazónica, la cuenca del Congo en África, y fragmentos del sudeste asiático insular. Esto demuestra visualmente que el microclima del Valle del Cauca es sumamente particular y estadísticamente “raro” a nivel mundial.

Implicación Agronómica y de Transferencia Tecnológica:

El análisis comparativo revela que una cosa es que una región sea apta para cultivar caña, y otra muy distinta es que tenga el mismo clima que el Valle del Cauca. Mientras que el modelo FAO abría la puerta a sembrar caña en casi todo Brasil o India, este mapa de analogía te indica que solo en las selectas áreas verdes podrías exportar e introducir las variedades de caña específicas (los clones) desarrolladas para el Valle del Cauca, esperando que tengan exactamente la misma tasa de crecimiento, fotosíntesis y respuesta al riego, sin necesidad de adaptar tu paquete tecnológico.

Comparación de Mapas

# Asegúrate de tener cargadas: terra, ggplot2, tidyterra, patchwork
library(terra)
library(ggplot2)
library(tidyterra)
library(patchwork) # Para panelar los mapas

cat("Generando imagen panelada de comparación...\n")
## Generando imagen panelada de comparación...
# --- MAPA 1: APTITUD FAO (Recreación rápida del Paso 1) ---
# Clasificación basada en rangos óptimos de temperatura y precipitación anual
aptitud_temp_opt <- ifel(tmean_annual >= 24 & tmean_annual <= 30, 2, 0)
aptitud_prec_opt <- ifel(prec_annual >= 1500 & prec_annual <= 2500, 2, 0)
# Factor limitante (mínimo)
aptitud_final_plot <- min(aptitud_temp_opt, aptitud_prec_opt)
# Convertir a factor para colores discretos
aptitud_final_plot <- as.factor(aptitud_final_plot)
levels(aptitud_final_plot) <- data.frame(ID=c(0, 2), Aptitud=c("No Apto", "Aptitud Óptima"))

plot1_aptitud <- ggplot() +
  geom_spatraster(data = aptitud_final_plot) +
  scale_fill_manual(
    values = c("No Apto" = "#E0E0E0", "Aptitud Óptima" = "#4CAF50"),
    na.translate = FALSE, name = "Nivel de Aptitud"
  ) +
  labs(title = "A. Modelo Teórico de Aptitud (Rangos FAO)",
       subtitle = "Zonas biológicamente viables para caña de azúcar") +
  theme_minimal() + theme(legend.position = "right", plot.title = element_text(face="bold"))
## <SpatRaster> resampled to 5e+05 cells.
# --- MAPA 2: SIMILARIDAD GLOBAL (Recreación del Paso 4) ---
# Cálculo de similaridad basada en la distancia euclidiana estandarizada (Punto 4)

# similaridad_global_plot <- 1 - (distancia_euclidiana / max_dist) # Calculado en el paso anterior

plot2_similaridad <- ggplot() +
  geom_spatraster(data = similaridad_global) + # Usamos el objeto del Paso 4
  scale_fill_gradientn(
    colors = c("#E0E0E0", "#FFC107", "#4CAF50"), # Gris -> Amarillo -> Verde
    name = "Índice de\nSimilaridad", na.value = "transparent"
  ) +
  labs(title = "B. Enfoque Empírico de Analogía Climática (Distancia Euclidiana)",
       subtitle = "Zonas climáticamente análogas al Valle del Cauca") +
  theme_minimal() + theme(legend.position = "right", plot.title = element_text(face="bold"))
## <SpatRaster> resampled to 5e+05 cells.
# --- ENSAMBLAJE FINAL CON PATCHWORK ---
panel_comparativo <- plot1_aptitud / plot2_similaridad + 
  plot_layout(heights = c(1, 1)) +
  plot_annotation(
    title = "Comparativo de Modelado Agroclimático para Caña de Azúcar",
    subtitle = "Teoría FAO vs. Analogía Empírica del Valle del Cauca",
    theme = theme(plot.title = element_text(size = 16, face = "bold", hjust = 0.5))
  )

print(panel_comparativo)

El contraste visual y analítico entre el Mapa A (Modelo Teórico FAO) y el Mapa B (Analogía Climática Empírica) evidencia la transición fundamental entre evaluar la viabilidad biológica de una especie y evaluar su viabilidad agroindustrial específica.

Filtro de Precisión:

El Mapa A demuestra que el trópico global es inmensamente permisivo para la supervivencia de la caña de azúcar, mostrando vastas extensiones de tierras aptas. Sin embargo, el Mapa B actúa como un filtro de alta resolución matemática. Al aplicar la Distancia Euclidiana, las grandes masas verdes se reducen drásticamente, demostrando que encontrar un clon climático del Valle del Cauca —con su delicado balance de isotermia cálida y lluvias moderadas bien distribuidas— es estadísticamente excepcional a nivel global.

##Toma de Decisiones: Mientras que el primer enfoque sirve para macro-políticas de expansión agrícola (identificar dónde la planta no morirá), el segundo enfoque es la herramienta definitiva para la transferencia tecnológica. Las zonas que permanecen de color verde intenso en el Mapa B indican que allí se podría exportar exactamente el mismo paquete agronómico del Valle del Cauca (mismas variedades genéticas desarrolladas en Cenicaña, mismos esquemas de fertilización y cronogramas de maduración) con una altísima probabilidad de replicar sus excepcionales rendimientos industriales, sin necesidad de costosas curvas de adaptación.

Top 5: Países con Alto Potencial (Análogos al Valle del Cauca)

# Asegúrate de tener cargadas: terra, geodata, ggplot2, tidyterra
library(terra)
library(geodata)
library(ggplot2)
library(tidyterra)

cat("Aislando los clústeres de alta similaridad...\n")
## Aislando los clústeres de alta similaridad...
# 1. Aplicar un filtro (Umbral estricto)
# Usamos el objeto 'similaridad_global' del paso 4.
# La función ifel() conserva los valores >= 0.75 y vuelve NA (nulo) todo lo demás.
cluster_nicho <- ifel(similaridad_global >= 0.75, similaridad_global, NA)

# 2. Descargar un mapa base mundial (polígonos) para dar contexto de fondo
mundo <- world(resolution = 5, path = tempdir())
## Cached as: C:\Users\ACER\AppData\Local\Temp\Rtmp8WaS6M/gadm/gadm36_adm0_r5_pk.rds
# 3. Visualización Cartográfica del Clúster
mapa_clusters_aislados <- ggplot() +
  # Capa 1: Continentes de fondo en gris claro (para ubicar visualmente los clústeres)
  geom_spatvector(data = mundo, fill = "#F5F5F5", color = "#D6D6D6", linewidth = 0.3) +
  
  # Capa 2: Superponer solo los píxeles del clúster (verde intenso)
  geom_spatraster(data = cluster_nicho) +
  scale_fill_gradient(
    low = "#8BC34A",   # Verde claro (Similaridad ~0.75)
    high = "#1B5E20",  # Verde muy oscuro (Similaridad casi 1.0)
    name = "Índice de\nAlta Similaridad",
    na.value = "transparent" # Hace invisible el resto del mundo
  ) +
  
  labs(
    title = "Clústeres Globales de Alta Similaridad (Nicho Analógico)",
    subtitle = "Aislamiento estricto (Índice > 0.75) respecto al Valle del Cauca",
    x = "Longitud",
    y = "Latitud"
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold", size = 14),
    panel.grid.major = element_line(color = "gray95"),
    panel.background = element_rect(fill = "white", color = NA)
  )
## <SpatRaster> resampled to 5e+05 cells.
# Imprimir el gráfico
print(mapa_clusters_aislados)

Observando estrictamente los clústeres de alta similaridad (tonos verdes con índice > 0.75 en el Mapa B) cruzados con la viabilidad del Mapa A, se identifican los siguientes 5 países con un potencial agroindustrial extraordinario para replicar este modelo:

Indonesia:

Presenta algunas de las manchas verdes más consistentes en el sudeste asiático (particularmente en Sumatra, Java y el sur de Borneo). Su ubicación ecuatorial garantiza la estabilidad térmica, y regiones específicas bajo la sombra orográfica de sus cadenas montañosas replican el régimen pluviométrico moderado necesario.

Filipinas:

Fuertemente iluminado en el mapa de similaridad, especialmente en las islas centrales y del sur (como Mindanao). Este país ya es un productor histórico, lo que valida empíricamente el modelo matemático; sin embargo, el mapa sugiere que la inserción de tecnología colombiana en estas zonas verdes específicas tendría un acople casi perfecto.

República Democrática del Congo (RDC):

Aunque el Mapa A pinta casi toda África central de verde, el Mapa B restringe la similaridad a zonas específicas en la periferia de la cuenca del Congo y hacia los Grandes Lagos. Aquí, la elevación topográfica modera el calor extremo de la selva baja, creando un clima cálido-moderado idéntico al de los valles interandinos colombianos.

Madagascar:

En el Mapa B se aprecia claramente una franja verde longitudinal en la costa oriental y norte de la isla. La interacción de los vientos alisios del Océano Índico con su relieve central crea microclimas de secano y temperaturas sumamente estables a lo largo del año, ideales para el metabolismo C4 de la caña.

Guatemala (y la vertiente del Pacífico Centroamericano):

El istmo centroamericano resalta en ambos mapas. Específicamente, las zonas costeras y piedemontes del lado del Pacífico en Guatemala (y extendiéndose hacia El Salvador y Costa Rica) ofrecen una radiación solar excepcionalmente alta combinada con suelos volcánicos ricos. Las señales de similaridad climática aquí indican un entorno donde las prácticas de zafra continua (o semi-continua) podrían prosperar con adaptación mínima.