library(terra)    
## terra 1.9.50
library(sf)         
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
library(geodata)  
library(ggplot2)   
library(tidyterra)
## 
## Adjuntando el paquete: 'tidyterra'
## The following object is masked from 'package:stats':
## 
##     filter
ruta_datos <- tempdir() 

prec_global <- worldclim_global(var = "prec", res = 10, path = ruta_datos)
## Cached as: C:\Users\RESPUE~1.COL\AppData\Local\Temp\RtmpA1Svam/climate/wc2.1_10m/wc2.1_10m_prec.zip
tavg_global <- worldclim_global(var = "tavg", res = 10, path = ruta_datos)
## Cached as: C:\Users\RESPUE~1.COL\AppData\Local\Temp\RtmpA1Svam/climate/wc2.1_10m/wc2.1_10m_tavg.zip
# Dimensiones y nombres de las capas (12 capas por raster, una por mes)
print(prec_global)
## class       : SpatRaster
## size        : 1080, 2160, 12  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## sources     : wc2.1_10m_prec_01.tif
##               wc2.1_10m_prec_02.tif
##               wc2.1_10m_prec_03.tif
##               ... and 9 more sources
## names       : wc2.1~ec_01, wc2.1~ec_02, wc2.1~ec_03, wc2.1~ec_04, wc2.1~ec_05, wc2.1~ec_06, ...
## min values  :           0,           0,           0,           0,           0,           0, ...
## max values  :         908,         793,         720,        1004,        2068,        2210, ...
print(tavg_global)
## class       : SpatRaster
## size        : 1080, 2160, 12  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## sources     : wc2.1_10m_tavg_01.tif
##               wc2.1_10m_tavg_02.tif
##               wc2.1_10m_tavg_03.tif
##               ... and 9 more sources
## names       : wc2.1~vg_01, wc2.1~vg_02, wc2.1~vg_03, wc2.1~vg_04, wc2.1~vg_05, wc2.1~vg_06, ...
## min values  :     -45.884,  -44.799999,   -57.92575,  -64.192497,  -64.811502,  -64.358252, ...
## max values  :     34.0095,    32.82425,     32.9095,    34.19375,    36.25325,     38.3555, ...

1.2 Cálculo de la aptitud climática

# Calcular variables anuales a partir de los datos mensuales
prec_anual <- sum(prec_global)
tavg_anual <- mean(tavg_global)

# Definir las condiciones óptimas (Reglas lógicas: 1 = Cumple, 0 = No cumple)
aptitud_prec <- prec_anual >= 1500 & prec_anual <= 3500
aptitud_temp <- tavg_anual >= 22.5 & tavg_anual <= 28

# Cruzar ambas condiciones (Aptitud Global)
aptitud_total <- aptitud_prec & aptitud_temp

# Limpiar el raster: Convertimos los valores 0 en NA (nulos) 
# Esto sirve para que al graficar, las zonas no aptas queden transparentes
aptitud_total[aptitud_total == 0] <- NA

1.3 Mapeo de la aptitud climática global

mundo <- world(path = ruta_datos)
## Cached as: C:\Users\RESPUE~1.COL\AppData\Local\Temp\RtmpA1Svam/gadm/gadm36_adm0_r5_pk.rds
mapa_global <- ggplot() +
  # Capa base vectorial (Países)
  geom_spatvector(data = mundo, fill = "gray90", color = "white", size = 0.2) +
  # Capa raster (Zonas aptas) convertida a factor para clasificación categórica
  geom_spatraster(data = as.factor(aptitud_total)) +
  # Escala de colores adecuada
  scale_fill_manual(values = c("1" = "#2ca25f"), # Verde oscuro para caña
                    na.value = "transparent",
                    labels = c("Área con potencial"),
                    name = "Aptitud") +
  theme_minimal() +
  labs(title = "Zonas con potencial climático global para Caña de Azúcar",
       subtitle = "T. Media: 22.5-28 °C | Precipitación anual: 1500-3500 mm",
       x = "Longitud",
       y = "Latitud") +
  theme(legend.position = "bottom")
## <SpatRaster> resampled to 5e+05 cells.
mapa_global

2. Recorte por países con alto potencial

Primero, calculamos con los datos cuáles son los 3 países con mayor superficie apta para la caña de azúcar según nuestros rangos.

# Extraer y sumar los píxeles de aptitud (valor 1) para cada país
aptitud_por_pais <- terra::extract(aptitud_total, mundo, fun = sum, na.rm = TRUE)

# Añadir el nombre del país al resultado
aptitud_por_pais$Pais <- mundo$NAME_0

# La columna con la suma toma el nombre de la capa base, la identificamos por su índice
col_suma <- colnames(aptitud_por_pais)[2]

# Ordenar los países de mayor a menor potencial
datos_ordenados <- aptitud_por_pais[order(-aptitud_por_pais[[col_suma]]), ]

# Extraer los nombres de los 3 primeros
top3_paises <- head(datos_ordenados$Pais, 3)

cat("Según los datos espaciales, los 3 países con mayor potencial son:\n")
## Según los datos espaciales, los 3 países con mayor potencial son:
print(top3_paises)
## [1] "Brazil"                           "Democratic Republic of the Congo"
## [3] "Indonesia"
# Mapa individual con contexto regional (Mantenemos la misma función)
crear_mapa_pais <- function(nombre_pais) {
  
  pais <- mundo[mundo$NAME_0 == nombre_pais, ]
  ext_ampliada <- ext(pais) + 15
  mundo_regional <- crop(mundo, ext_ampliada)
  
  aptitud_regional <- crop(aptitud_total, ext_ampliada)
  aptitud_pais_solo <- mask(aptitud_regional, pais)
  
  mapa <- ggplot() +
    geom_spatvector(data = mundo_regional, fill = "gray90", color = "white", size = 0.3) +
    geom_spatvector(data = pais, fill = NA, color = "black", linewidth = 0.8) +
    geom_spatraster(data = as.factor(aptitud_pais_solo)) +
    scale_fill_manual(values = c("1" = "#2ca25f"), 
                      na.value = "transparent",
                      labels = c("Área con potencial"),
                      name = "Aptitud") +
    theme_minimal() +
    labs(title = paste("Potencial Climático para Caña de Azúcar"),
         subtitle = paste("Enfoque continental:", nombre_pais),
         x = "Longitud",
         y = "Latitud") +
    theme(legend.position = "bottom")
  
  return(mapa)
}

# Aunque el top 3 por área es Brasil, Congo e Indonesia, 
# forzamos la lista para incluir a Colombia por requerimiento del caso de estudio.
paises_analisis <- c("Colombia", "Brazil", "Indonesia")

mapa_col <- crear_mapa_pais(paises_analisis[1])
mapa_bra <- crear_mapa_pais(paises_analisis[2])
mapa_ind <- crear_mapa_pais(paises_analisis[3])

mapa_col

mapa_bra

mapa_ind

3. Extracción de series de tiempo climáticas en el Valle del Cauca

Primero, cargamos las librerías para manipulación de datos y definimos nuestros puntos de interés.

Punto 1 (Sur): Cerca a Candelaria/Palmira (Lat: 3.45, Lon: -76.35)

Punto 2 (Centro): Cerca a Buga (Lat: 3.88, Lon: -76.30)

Punto 3 (Norte): Cerca a Zarzal (Lat: 4.35, Lon: -76.05)

library(dplyr)
## 
## 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(tidyr)
## 
## Adjuntando el paquete: 'tidyr'
## The following object is masked from 'package:terra':
## 
##     extract
puntos_valle <- data.frame(
  id = c("Punto 1 (Sur)", "Punto 2 (Centro)", "Punto 3 (Norte)"),
  lon = c(-76.35, -76.30, -76.05),
  lat = c(3.45, 3.88, 4.35)
)

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

# Extraer los datos mensuales de temperatura y precipitación

temp_ext <- terra::extract(tavg_global, puntos_sf)
prec_ext <- terra::extract(prec_global, puntos_sf)

# Preparar los datos para graficar (formato largo)
temp_ext$Punto <- puntos_valle$id
prec_ext$Punto <- puntos_valle$id

meses <- c("Ene", "Feb", "Mar", "Abr", "May", "Jun", "Jul", "Ago", "Sep", "Oct", "Nov", "Dic")

# Transformar Temperatura a formato largo
temp_long <- temp_ext %>%
  select(-ID) %>% # Quitamos la columna ID que genera extract()
  pivot_longer(cols = -Punto, names_to = "Capa", values_to = "Temperatura") %>%
  mutate(Mes = rep(1:12, times = 3),
         Mes_nombre = factor(meses[Mes], levels = meses))

# Transformar Precipitación a formato largo
prec_long <- prec_ext %>%
  select(-ID) %>%
  pivot_longer(cols = -Punto, names_to = "Capa", values_to = "Precipitacion") %>%
  mutate(Mes = rep(1:12, times = 3),
         Mes_nombre = factor(meses[Mes], levels = meses))


grafico_temp <- ggplot(temp_long, aes(x = Mes_nombre, y = Temperatura, group = Punto, color = Punto)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  theme_minimal() +
  labs(title = "Comportamiento Mensual de la Temperatura",
       subtitle = "Puntos de muestreo en el Valle del Cauca",
       x = "Mes", y = "Temperatura (°C)") +
  theme(legend.position = "bottom")


grafico_prec <- ggplot(prec_long, aes(x = Mes_nombre, y = Precipitacion, group = Punto, fill = Punto)) +
  geom_col(position = "dodge") +
  theme_minimal() +
  labs(title = "Comportamiento Mensual de la Precipitación",
       subtitle = "Puntos de muestreo en el Valle del Cauca",
       x = "Mes", y = "Precipitación (mm)") +
  theme(legend.position = "bottom")


grafico_temp

grafico_prec

4. Mapa de similaridad climática global (Distancia Euclidiana)

clima_global <- c(tavg_global, prec_global)

clima_norm <- scale(clima_global)

# Extraer los valores climáticos estandarizados de los 3 puntos del Valle del Cauca
perfil_puntos <- terra::extract(clima_norm, puntos_sf)

# Calcular el "perfil climático promedio" del Valle (promediamos los 3 puntos)
perfil_promedio <- colMeans(perfil_puntos[, -1])

# Calcular la Distancia Euclidiana entre todo el mundo y el perfil promedio del Valle
distancia_euclidiana <- sqrt(sum((clima_norm - perfil_promedio)^2))

# Graficar el mapa de similaridad global
mapa_similaridad <- ggplot() +
  geom_spatvector(data = mundo, fill = "gray90", color = "white", size = 0.2) +
  geom_spatraster(data = distancia_euclidiana) +
  # Usamos una paleta de colores invertida (direction = -1) 
  # para que los colores intensos/oscuros representen las zonas más similares (menor distancia)
  scale_fill_viridis_c(direction = -1, 
                       option = "magma", 
                       na.value = "transparent",
                       name = "Distancia\nEuclidiana\n(Menor = Más similar)") +
  theme_minimal() +
  labs(title = "Similaridad climática global con el Valle del Cauca, Colombia",
       subtitle = "Métrica: Distancia Euclidiana (24 variables estandarizadas de temp y prec)",
       x = "Longitud",
       y = "Latitud") +
  theme(legend.position = "bottom")
## <SpatRaster> resampled to 5e+05 cells.
# Mostrar el mapa
mapa_similaridad

5. Comparación de enfoques y conclusiones

Al contrastar los mapas generados por ambas aproximaciones metodológicas, se evidencian diferencias sustanciales en la delimitación de zonas con potencial para el cultivo de caña de azúcar:

1. Aproximación por Rangos Óptimos (Filtro Binario) * Resultados: El primer mapa muestra enormes extensiones de tierra catalogadas como aptas (color verde), concentrándose fuertemente en las cuencas del Amazonas en Brasil, el Congo en África y el archipiélago de Indonesia[cite: 5, 8]. * Análisis: Esta metodología utiliza umbrales estrictos anuales (temperatura entre 22.5-28 °C y precipitación entre 1500-3500 mm)[cite: 5]. Al ser un enfoque binario (apto/no apto), generaliza demasiado. Un país podría cumplir con los 3000 mm de lluvia anual, pero concentrados en solo tres meses (provocando inundaciones), y el modelo lo seguiría clasificando como ideal.

2. Aproximación por Similaridad Climática (Distancia Euclidiana) * Resultados: El mapa de calor global muestra zonas de alta similaridad (colores oscuros) mucho más restringidas y específicas[cite: 11]. * Análisis: A diferencia del primer enfoque, este método evalúa 24 variables simultáneas (12 meses de temperatura y 12 de precipitación)[cite: 11]. Al tomar como referencia el Valle del Cauca, el modelo no solo busca volúmenes totales, sino que exige que se replique la estacionalidad: una temperatura estable todo el año (21-24 °C)[cite: 9] y un régimen de lluvias bimodal (con picos alrededor de abril y octubre)[cite: 10].

Conclusión Final La aproximación por rangos teóricos es un excelente primer filtro macro-regional, pero tiende a sobreestimar las áreas viables[cite: 5]. Por otro lado, la aproximación por similaridad es mucho más robusta para transferir paquetes tecnológicos, ya que encontrar un “gemelo climático” del Valle del Cauca garantiza que las prácticas agronómicas, el calendario de siembra y las variedades genéticas de la caña de azúcar tendrán un comportamiento casi idéntico al de la zona de referencia. Lo ideal en un estudio de zonificación agrícola es combinar ambas: usar los rangos para delimitar la región, y la distancia euclidiana para priorizar los polígonos de siembra.