1 Introducción

Este documento presenta un análisis geoespacial de aptitud climática global para el cultivo de caña de azúcar (Saccharum officinarum), desarrollado en el marco de un ejercicio académico de exploración de datos espaciales en la Maestría en Ciencia de Datos.

El análisis integra dos aproximaciones complementarias:

  1. Aproximación por rangos óptimos: identifica zonas donde temperatura media y precipitación cumplen simultáneamente los umbrales agronómicos conocidos para el cultivo. Se generan tres mapas de aptitud: uno basado en precipitación anual, uno basado en precipitación mensual, y uno combinado final.
  2. Aproximación por similaridad climática: utiliza la distancia euclidiana normalizada en el espacio climático (24 variables) para identificar zonas globalmente similares a tres puntos de referencia del Valle del Cauca, Colombia, región con condiciones conocidas de alta productividad cañera.

Para evitar una selección subjetiva de países, se calculó el área climáticamente apta por país mediante la superposición del raster de aptitud con el shape global de países. Esta aproximación permitió identificar los países con mayor extensión de áreas potenciales a partir de los resultados del modelo espacial, y no únicamente desde conocimiento previo sobre producción agrícola.

La comparación entre métodos se complementó con métricas de superposición positiva (índice de Jaccard, precisión y recall), debido a que la concordancia general puede estar inflada por la gran cantidad de celdas no aptas en ambos métodos.

Los datos climáticos provienen de WorldClim 2.1 (línea base 1970–2000), descargados mediante el paquete geodata.

2 Rangos Óptimos para Caña de Azúcar

Los rangos óptimos de aptitud climática utilizados son:

Variable Rango óptimo
Temperatura media mensual 22.5 – 28 °C
Precipitación anual 1.500 – 3.500 mm
Alternativa: Precipitación mensual 125 – 290 mm

Criterios de precipitación utilizados en el análisis:

El enunciado establece que la aptitud puede cumplirse por precipitación anual o por precipitación mensual. Para transparencia metodológica, se aplican tres criterios de forma explícita:

  • Criterio 1 — Precipitación anual: la suma de los 12 meses debe estar entre 1.500 y 3.500 mm.
  • Criterio 2 — Precipitación mensual: al menos 6 de los 12 meses deben tener precipitación entre 125 y 290 mm (criterio flexible que reconoce que no todos los meses de un año agrícola requieren precipitación óptima simultáneamente).
  • Criterio combinado: un píxel es apto si cumple el criterio anual O el criterio mensual.

Se generan mapas separados para cada criterio, lo que permite al agricultor evaluar la robustez espacial del resultado.

3 Carga de Librerías

Se utilizan los paquetes del enfoque visto en clase. terra es el sucesor moderno de raster y rgdal, desarrollado por el mismo autor (Robert Hijmans). Se complementa con geodata para descarga de WorldClim, rnaturalearth para el shape global, y ggplot2/tidyr para visualización.

# ── Paquetes espaciales  ──────────────────────────────────
library(terra)        # Manejo moderno de datos raster y vector
library(raster)       
library(sp)           # Clases espaciales clásicas

# ── Descarga de datos climáticos y cartografía ──────────────────────────────
library(geodata)      # Descarga WorldClim, GADM, etc.
library(rnaturalearth)   # Shapes de países del mundo
library(rnaturalearthdata)

# ── Visualización ───────────────────────────────────────────────────────────
library(ggplot2)      # Gráficos
library(tidyr)        # Reshape de datos
library(dplyr)        # Manipulación de datos
library(viridis)      # Paletas de color perceptualmente uniformes
library(RColorBrewer) # Paletas adicionales
library(sf)           # Datos vectoriales (simple features)
library(patchwork)    # Composición de gráficos múltiples

# ── Tablas ───────────────────────────────────────────────────────────────────
library(knitr)
library(kableExtra)

4 Descarga y Carga de Datos Climáticos WorldClim

Se descarga WorldClim 2.1 a resolución de 10 minutos (~18 km). Las variables son temperatura media mensual (tavg) y precipitación mensual (prec).

Nota: La primera ejecución descarga los archivos (~200 MB). Las ejecuciones siguientes usan la caché local en data/worldclim/.

# ── Directorio de datos ──────────────────────────────────────────────────────
dir.create("data/worldclim", recursive = TRUE, showWarnings = FALSE)

# ── Temperatura media mensual (°C × 10, se corrige abajo) ───────────────────
# geodata descarga 12 capas mensuales en un SpatRaster
tavg_raw <- worldclim_global(
  var  = "tavg",   # temperatura media
  res  = 10,       # resolución: 10 minutos de arco
  path = "data/worldclim"
)

# ── Precipitación mensual (mm) ───────────────────────────────────────────────
prec_raw <- worldclim_global(
  var  = "prec",
  res  = 10,
  path = "data/worldclim"
)

# ── Verificar objetos ────────────────────────────────────────────────────────
cat("Temperatura (SpatRaster):\n"); print(tavg_raw)
## Temperatura (SpatRaster):
## 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.8840,   -44.80000,   -57.92575,   -64.19250,   -64.81150,   -64.35825, ... 
## max values  :     34.0095,    32.82425,    32.90950,    34.19375,    36.25325,    38.35550, ...
cat("\nPrecipitacion (SpatRaster):\n"); print(prec_raw)
## 
## Precipitacion (SpatRaster):
## 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, ...
# ── Cargar directamente desde disco (patrón del código de clase) ───────────
archivos_temp <- list.files("data/worldclim/climate/wc2.1_10m",
                             pattern = "tavg", 
                             full.names = TRUE)

archivos_prec <- list.files("data/worldclim/climate/wc2.1_10m",
                             pattern = "prec",
                             full.names = TRUE)

# ── Ordenar para garantizar orden mensual 01–12 ────────────────────────────
archivos_temp <- sort(archivos_temp)
archivos_prec <- sort(archivos_prec)

# ── Crear SpatRaster con terra::rast() ────────────────────────────────────
temperaturas  <- rast(archivos_temp)
precipitacion <- rast(archivos_prec)

# ── Asignar nombres simples ────────────────────────────────────────────────
names(temperaturas)  <- paste0("tmean_", sprintf("%02d", 1:12))
names(precipitacion) <- paste0("prec_",  sprintf("%02d", 1:12))

cat("Capas temperatura:",   nlyr(temperaturas),  "\n")
## Capas temperatura: 12
cat("Capas precipitacion:", nlyr(precipitacion), "\n")
## Capas precipitacion: 12
# ── Convertir a SpatRaster para todo el flujo en terra ────────────────────
temperaturas  <- rast(archivos_temp)
precipitacion <- rast(archivos_prec)
names(temperaturas)  <- paste0("tmean_", sprintf("%02d", 1:12))
names(precipitacion) <- paste0("prec_",  sprintf("%02d", 1:12))

# ── Temperatura media anual ───────────────────────────────────────────────
tmean_anual <- terra::mean(temperaturas)
names(tmean_anual) <- "tmean_anual"

# ── Precipitacion anual ───────────────────────────────────────────────────
prec_anual <- terra::app(precipitacion, fun = sum)
names(prec_anual) <- "prec_anual"

# ── Verificar clase ───────────────────────────────────────────────────────
cat("Clase tmean_anual:", class(tmean_anual), "\n")
## Clase tmean_anual: SpatRaster
cat("Clase prec_anual:",  class(prec_anual),  "\n")
## Clase prec_anual: SpatRaster
# ── Visualizacion ─────────────────────────────────────────────────────────
par(mfrow = c(1, 2))
plot(tmean_anual,
     main = "Temperatura Media Anual (grados C)",
     col  = rev(RColorBrewer::brewer.pal(11, "RdYlBu")))
plot(prec_anual,
     main = "Precipitacion Anual (mm)",
     col  = RColorBrewer::brewer.pal(9, "Blues"))

par(mfrow = c(1, 1))

La visualización de las variables climáticas anuales permite explorar la distribución espacial de los dos factores principales utilizados en el análisis de aptitud climática: la temperatura media anual y la precipitación anual acumulada. En el mapa de temperatura se observa un patrón latitudinal claro, con valores más altos en las zonas tropicales y subtropicales, y temperaturas más bajas hacia las regiones polares y de alta montaña. Por su parte, el mapa de precipitación anual muestra una distribución más heterogénea, con mayores acumulados en regiones ecuatoriales, selváticas y zonas húmedas, mientras que las áreas desérticas y semiáridas presentan valores bajos. Esta exploración inicial es importante porque permite comprender la base climática sobre la cual se construyen posteriormente los mapas de aptitud para caña de azúcar, ya que el cultivo requiere temperaturas cálidas y una disponibilidad hídrica suficiente para alcanzar condiciones favorables de desarrollo.

5 Mapas de Aptitud Climática Global

5.1 Criterios de aptitud

La aptitud climática se evalúa mediante álgebra raster. Un píxel es apto si cumple:

  • Temperatura: tmean_anual ∈ [22.5, 28] °C, Y
  • Precipitación: criterio anual O criterio mensual (ver sección de rangos óptimos).
# ── Máscara de temperatura apta ───────────────────────────────────────────
apt_temp <- (tmean_anual >= 22.5) & (tmean_anual <= 28)
names(apt_temp) <- "apt_temp"

cat("Celdas con temperatura apta:", 
    global(apt_temp, sum, na.rm = TRUE)[[1]], "\n")
## Celdas con temperatura apta: 139273
# ── Criterio 1: Precipitación ANUAL apta ─────────────────────────────────
apt_prec_anual <- (prec_anual >= 1500) & (prec_anual <= 3500)

# ── Criterio 2: Precipitación MENSUAL — capa por capa con escritura disco
dir.create("data/worldclim/temp", showWarnings = FALSE)

# Inicializar contador en disco
n_meses_opt <- prec_anual * 0
writeRaster(n_meses_opt, "data/worldclim/temp/n_meses.tif", overwrite = TRUE)
rm(n_meses_opt); invisible(gc())

for (i in 1:12) {
  capa        <- precipitacion[[i]]
  n_meses_opt <- rast("data/worldclim/temp/n_meses.tif")
  n_meses_opt <- n_meses_opt + ((capa >= 125) & (capa <= 290))
  writeRaster(n_meses_opt, "data/worldclim/temp/n_meses.tif", overwrite = TRUE)
  rm(n_meses_opt, capa); invisible(gc())
}

n_meses_opt      <- rast("data/worldclim/temp/n_meses.tif")
apt_prec_mensual <- (n_meses_opt >= 6)
rm(n_meses_opt); invisible(gc())

# ── Criterio combinado: anual O mensual ───────────────────────────────────
apt_prec <- apt_prec_anual | apt_prec_mensual
names(apt_prec_anual)   <- "apt_prec_anual"
names(apt_prec_mensual) <- "apt_prec_mensual"
names(apt_prec)         <- "apt_prec"

cat("Celdas aptas — criterio anual:   ", 
    global(apt_prec_anual,   sum, na.rm = TRUE)[[1]], "\n")
## Celdas aptas — criterio anual:    68382
cat("Celdas aptas — criterio mensual: ", 
    global(apt_prec_mensual, sum, na.rm = TRUE)[[1]], "\n")
## Celdas aptas — criterio mensual:  52994
cat("Celdas aptas — criterio combinado:", 
    global(apt_prec,         sum, na.rm = TRUE)[[1]], "\n")
## Celdas aptas — criterio combinado: 79211
# ── Aptitud global: temperatura AND precipitación ─────────────────────────
aptitud_global <- apt_temp & apt_prec
names(aptitud_global) <- "aptitud_cana"
aptitud_global <- mask(aptitud_global, tmean_anual)

# ── Área en km² (más informativo que conteo de celdas) ────────────────────
# cellSize() devuelve el área de cada celda en m²
area_celdas <- cellSize(aptitud_global, unit = "km")
area_apta_km2 <- global(area_celdas * aptitud_global, sum, na.rm = TRUE)[[1]]

cat("Celdas con aptitud climática global:", 
    global(aptitud_global, sum, na.rm = TRUE)[[1]], "\n")
## Celdas con aptitud climática global: 56531
cat(sprintf("Área climáticamente apta estimada: %.0f km²\n", area_apta_km2))
## Área climáticamente apta estimada: 19051926 km²
cat(sprintf("Equivale a aproximadamente %.1f veces el área de Colombia\n", 
            area_apta_km2 / 1141748))
## Equivale a aproximadamente 16.7 veces el área de Colombia

5.2 Mapas de aptitud por criterio de precipitación

Se generan tres mapas para transparencia metodológica:

# ── Cargar shape global si aún no existe ─────────────────────────────────
if (!exists("world_sf")) {
  world_sf <- ne_countries(scale = "medium", returnclass = "sf")
}
# ── Función auxiliar de mapa binario ──────────────────────────────────────
mapa_binario <- function(r_prec, titulo, subtitulo) {
  # Aplicar máscara de temperatura directamente
  r_final <- (apt_temp & r_prec)
  r_final <- mask(r_final, tmean_anual)
  
  df <- as.data.frame(r_final, xy = TRUE, na.rm = TRUE)
  names(df) <- c("lon", "lat", "apto")
  df$apto_cat <- factor(df$apto, levels = c(0,1), 
                        labels = c("No apto", "Apto"))
  
  ggplot() +
    geom_tile(data = df, aes(x = lon, y = lat, fill = apto_cat)) +
    geom_sf(data = world_sf, fill = NA, color = "gray30", linewidth = 0.15) +
    scale_fill_manual(values = c("No apto" = "#f5f5f5", "Apto" = "#1a9850"),
                      name = "Aptitud") +
    coord_sf(expand = FALSE) +
    labs(title = titulo, subtitle = subtitulo,
         x = "Longitud", y = "Latitud", caption = "Fuente: WorldClim 2.1") +
    theme_minimal(base_size = 10) +
    theme(panel.background = element_rect(fill = "#d6eaf8", color = NA),
          plot.title = element_text(face = "bold"),
          legend.position = "bottom")
}

# Mapa 1: Criterio precipitacion anual
m1 <- mapa_binario(apt_prec_anual,
                   "Mapa 1 — Aptitud con Precipitacion Anual",
                   "Temp ∈ [22.5,28°C] Y Prec_anual ∈ [1500,3500 mm]")

# Mapa 2: Criterio precipitacion mensual
m2 <- mapa_binario(apt_prec_mensual,
                   "Mapa 2 — Aptitud con Precipitacion Mensual",
                   "Temp ∈ [22.5,28°C] Y al menos 6 meses con Prec ∈ [125,290 mm]")

# Mapa 3: Criterio combinado
m3 <- mapa_binario(apt_prec,
                   "Mapa 3 — Aptitud Combinada (Criterio Final)",
                   "Temp ∈ [22.5,28°C] Y (Criterio anual O Criterio mensual)")

m1 / m2 / m3

Los mapas de aptitud por criterio de precipitación permiten comparar cómo cambia la identificación de zonas climáticamente favorables para la caña de azúcar según la forma en que se evalúa la disponibilidad hídrica. El primer mapa, basado en la precipitación anual, resalta las áreas que cumplen simultáneamente con el rango óptimo de temperatura media y con una precipitación acumulada entre 1.500 y 3.500 mm al año, mostrando zonas aptas principalmente en regiones tropicales de Suramérica, África ecuatorial y el sudeste asiático. El segundo mapa, construido a partir de la precipitación mensual, incorpora un criterio más detallado al considerar la distribución de la lluvia durante el año, lo que permite identificar áreas donde las condiciones hídricas son más constantes o se aproximan mejor a los requerimientos del cultivo. Finalmente, el tercer mapa integra ambos enfoques y permite reconocer las zonas con mayor consistencia climática para la caña de azúcar. En conjunto, estos mapas evidencian que la precipitación no solo debe analizarse por su acumulado anual, sino también por su comportamiento mensual, ya que una misma cantidad de lluvia puede tener implicaciones diferentes según su distribución temporal.

5.3 Visualización con ggplot2

# ── Convertir a data.frame para ggplot2 (flujo estándar en clase) ─────────────
df_apt <- as.data.frame(aptitud_global, xy = TRUE, na.rm = TRUE)
names(df_apt) <- c("lon", "lat", "aptitud")
df_apt$aptitud_cat <- factor(df_apt$aptitud,
                              levels = c(0, 1),
                              labels = c("No apto", "Apto"))

# ── Contorno de países (rnaturalearth) ────────────────────────────────────────
world_sf <- ne_countries(scale = "medium", returnclass = "sf")

# ── Mapa ──────────────────────────────────────────────────────────────────────
ggplot() +
  geom_raster(data = df_apt,
              aes(x = lon, y = lat, fill = aptitud_cat)) +
  geom_sf(data = world_sf,
          fill  = NA,
          color = "gray30",
          linewidth = 0.2) +
  scale_fill_manual(
    values = c("No apto" = "#f0f0f0", "Apto" = "#1a9850"),
    name   = "Aptitud climatica"
  ) +
  coord_sf(expand = FALSE) +
  labs(
    title    = "Aptitud Climatica Global para Cana de Azucar",
    subtitle = "Criterio: Tmean ∈ [22.5, 28°C] y Prec_anual ∈ [1500, 3500 mm]",
    x = "Longitud", y = "Latitud",
    caption  = "Fuente: WorldClim 2.1 (1970–2000)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position  = "bottom",
    panel.background = element_rect(fill = "#d6eaf8", color = NA),
    plot.title       = element_text(face = "bold", size = 14)
  )

La visualización construida con ggplot2 permite representar de manera clara las zonas del mundo que cumplen simultáneamente con los rangos óptimos de temperatura media y precipitación anual definidos para el cultivo de caña de azúcar. En el mapa, las áreas resaltadas en verde corresponden a las regiones climáticamente aptas, mientras que el resto del territorio aparece sin clasificación de aptitud. Los resultados muestran una concentración de zonas favorables en áreas tropicales y subtropicales, especialmente en el norte de Suramérica, América Central, África ecuatorial y el sudeste asiático. Esta distribución es coherente con los requerimientos agroclimáticos de la caña de azúcar, ya que el cultivo demanda temperaturas cálidas y niveles adecuados de precipitación. El uso de ggplot2 facilita además una presentación cartográfica ordenada, con límites políticos, coordenadas geográficas y una escala visual que permite interpretar espacialmente los resultados obtenidos.

5.4 Mapa con escala de aptitud continua (índice)

Para mayor riqueza informativa, se construye un índice de aptitud continuo (0–1) que combina la “cercanía” de cada celda al rango óptimo de temperatura y precipitación.

# ── Índice de temperatura (0–1) ───────────────────────────────────────────
temp_center <- 25.25
temp_half   <- 2.75

idx_temp <- 1 - abs(tmean_anual - temp_center) / temp_half
idx_temp <- clamp(idx_temp, lower = 0, upper = 1)

# ── Índice de precipitación anual (0–1) ───────────────────────────────────
prec_center <- 2500
prec_half   <- 1000

idx_prec <- 1 - abs(prec_anual - prec_center) / prec_half
idx_prec <- clamp(idx_prec, lower = 0, upper = 1)

# ── Índice combinado ──────────────────────────────────────────────────────
idx_total <- (idx_temp + idx_prec) / 2
names(idx_total) <- "indice_aptitud"

cat("Clase idx_total:", class(idx_total), "\n")
## Clase idx_total: SpatRaster
cat("Rango idx_total:", global(idx_total, range, na.rm = TRUE)[[1]], "\n")
## Rango idx_total: 0
# ── Cargar shape global si aún no existe ─────────────────────────────────
if (!exists("world_sf")) {
  world_sf <- ne_countries(scale = "medium", returnclass = "sf")
}

# ── Graficar ──────────────────────────────────────────────────────────────
df_idx <- as.data.frame(idx_total, xy = TRUE, na.rm = TRUE)
names(df_idx) <- c("lon", "lat", "indice")

ggplot() +
  geom_tile(data = df_idx,
              aes(x = lon, y = lat, fill = indice)) +
  geom_sf(data = world_sf,
          fill = NA, color = "gray30", linewidth = 0.2) +
  scale_fill_gradientn(
    colors = c("#d73027", "#fee08b", "#1a9850"),
    na.value = "white",
    name   = "Índice de\naptitud (0–1)"
  ) +
  coord_sf(expand = FALSE) +
  labs(
    title    = "Indice de Aptitud Climatica Continuo — Cana de Azucar",
    subtitle = "Verde oscuro = condiciones mas cercanas al optimo",
    x = "Longitud", y = "Latitud",
    caption  = "Fuente: WorldClim 2.1 (1970–2000)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position  = "right",
    panel.background = element_rect(fill = "#d6eaf8", color = NA),
    plot.title       = element_text(face = "bold", size = 14)
  )

El mapa de aptitud climática continua permite observar, con mayor nivel de detalle, qué tan cercanas son las condiciones climáticas de cada zona del mundo a los rangos óptimos definidos para el cultivo de caña de azúcar. A diferencia del mapa binario de apto/no apto, este índice muestra una transición gradual entre zonas con baja, media y alta aptitud, donde los colores verdes representan áreas con condiciones más favorables y cercanas al óptimo, mientras que los tonos rojos indican territorios con menor aptitud climática. Se evidencia una mayor concentración de valores altos en regiones tropicales y subtropicales, especialmente en sectores de América Central, norte de Suramérica, África ecuatorial y el sudeste asiático. Este resultado es coherente con los requerimientos climáticos de la caña de azúcar, ya que el cultivo tiende a desarrollarse mejor en zonas cálidas y con disponibilidad adecuada de precipitación. Por tanto, el índice continuo complementa el análisis de aptitud al permitir priorizar áreas con mayor potencial relativo, incluso dentro de regiones que no necesariamente cumplen de forma estricta todos los umbrales óptimos.

6 Ranking de Países por Área Climáticamente Apta

Para evitar una selección subjetiva, se calcula el área apta por país a partir del raster de aptitud global, y se seleccionan objetivamente los 3 países con mayor área apta.

# ── Shape global como SpatVector ─────────────────────────────────────────
world_vect <- vect(world_sf)

# ── Área apta por celda ───────────────────────────────────────────────────
apt_num     <- as.numeric(aptitud_global)
area_cel    <- cellSize(apt_num, unit = "km")
area_apta_r <- area_cel * apt_num

# ── Suma de área apta por país ────────────────────────────────────────────
ranking_raw <- terra::extract(area_apta_r, world_vect, fun = "sum", na.rm = TRUE)

# ── Área total por país (para calcular porcentaje) ────────────────────────
area_total_r  <- cellSize(apt_num, unit = "km")
ranking_total <- terra::extract(area_total_r, world_vect, fun = "sum", na.rm = TRUE)

# ── Construir tabla de ranking con porcentaje ─────────────────────────────
ranking_df <- data.frame(
  pais           = world_vect$name_long,
  iso_a3         = world_vect$iso_a3,
  continente     = world_vect$continent,
  area_apta_km2  = ranking_raw[, 2],
  area_total_km2 = ranking_total[, 2]
) %>%
  filter(!is.na(area_apta_km2), area_apta_km2 > 0) %>%
  mutate(pct_apto = round(area_apta_km2 / area_total_km2 * 100, 1)) %>%
  arrange(desc(area_apta_km2)) %>%
  mutate(ranking = row_number())

# ── Top 15 ────────────────────────────────────────────────────────────────
top15 <- ranking_df %>% slice(1:15)

ggplot(top15, aes(x = reorder(pais, area_apta_km2), y = area_apta_km2 / 1000)) +
  geom_col(fill = "#1a9850", alpha = 0.85) +
  coord_flip() +
  labs(
    title   = "Top 15 Paises por Area Climaticamente Apta - Cana de Azucar",
    subtitle = "Area donde Tmean entre 22.5-28 grados y criterio de precipitacion cumplido",
    x       = NULL,
    y       = "Area apta (miles de km2)",
    caption = "Fuente: WorldClim 2.1 | Calculo propio"
  ) +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold"))

El ranking de países por área climáticamente apta permite identificar, de manera objetiva, los territorios con mayor extensión potencial para el cultivo de caña de azúcar según los criterios climáticos definidos. Los resultados muestran que Brasil presenta la mayor superficie apta, seguido por la República Democrática del Congo e Indonesia, lo cual es coherente con su ubicación en zonas tropicales y con la presencia de condiciones cálidas y húmedas favorables para el cultivo. También aparecen países de Suramérica como Colombia, Perú, Venezuela y Bolivia, así como países de África ecuatorial y del sudeste asiático, regiones donde la temperatura media y la precipitación se acercan a los rangos óptimos establecidos. Este ranking complementa los mapas globales, ya que permite pasar de una lectura visual a una comparación cuantitativa entre países, facilitando la selección de zonas prioritarias para un análisis espacial más detallado.

kable(
  top15 %>% 
    select(ranking, pais, continente, area_apta_km2, pct_apto) %>%
    mutate(area_apta_km2 = round(area_apta_km2, 0)),
  caption   = "Ranking de paises por area climaticamente apta para cana de azucar",
  col.names = c("Rank", "Pais", "Continente", "Área apta (km²)", "% del país apto"),
  format.args = list(big.mark = ",")
) %>%
  kable_styling(bootstrap_options = c("striped","hover")) %>%
  row_spec(1:3, bold = TRUE, background = "#d9f0d3")
Ranking de paises por area climaticamente apta para cana de azucar
Rank Pais Continente Área apta (km²) % del país apto
1 Brazil South America 5,270,367 62.2
2 Democratic Republic of the Congo Africa 1,750,732 75.2
3 Indonesia Asia 1,409,107 76.0
4 Colombia South America 795,594 70.0
5 Peru South America 628,561 48.7
6 Venezuela South America 579,758 63.7
7 Central African Republic Africa 446,264 72.3
8 Bolivia South America 351,156 32.4
9 Republic of the Congo Africa 332,065 96.3
10 Cameroon Africa 329,222 71.2
11 India Asia 319,641 10.2
12 Nigeria Africa 294,953 32.5
13 Myanmar Asia 282,991 42.8
14 Papua New Guinea Oceania 271,315 58.4
15 Philippines Asia 266,531 93.5 
# ── Seleccionar los 3 países con mayor área apta ──────────────────────────
top3_nombres <- ranking_df$pais[1:3]
cat("Países top 3 por área apta:\n")
## Países top 3 por área apta:
cat(paste0(1:3, ". ", top3_nombres, 
           " (", round(ranking_df$area_apta_km2[1:3]/1000, 0), " mil km²)"), 
    sep = "\n")
## 1. Brazil (5270 mil km²)
## 2. Democratic Republic of the Congo (1751 mil km²)
## 3. Indonesia (1409 mil km²)
# Usar iso_a3 para evitar inconsistencias de nombres
top3_iso <- ranking_df$iso_a3[1:3]

# Verificar si Colombia está en el top 3; si no, incluirla como referencia
col_en_top3 <- "COL" %in% top3_iso
if (!col_en_top3) {
  cat("\nColombia no está en el top 3. Se incluye como país de referencia\n")
  cat("por contener el Valle del Cauca, región central del análisis.\n")
  # Reemplazar el 3er país por Colombia para mantener la coherencia con el caso
  top3_iso[3] <- "COL"
}
## 
## Colombia no está en el top 3. Se incluye como país de referencia
## por contener el Valle del Cauca, región central del análisis.
paises_sf      <- world_sf %>% filter(iso_a3 %in% top3_iso)
paises_interes <- paises_sf$name_long
cat("\nPaíses finales para corte espacial:", paste(paises_interes, collapse=", "), "\n")
## 
## Países finales para corte espacial: Democratic Republic of the Congo, Colombia, Brazil

7 Corte Espacial para los 3 Países Seleccionados

Los dos primeros países provienen directamente del ranking objetivo por área climáticamente apta. Colombia se incluye aunque no esté necesariamente en el top 3, porque contiene el Valle del Cauca — la región de referencia del análisis de similaridad —, lo que permite conectar directamente la aptitud climática con el caso planteado.

Aunque el ranking por área apta puede ubicar a otros países por encima de Colombia, para el análisis de corte espacial se seleccionan los dos primeros del ranking y Colombia. Los países de mayor área potencial representan la escala global del análisis, mientras que Colombia se incluye como país de referencia por contener la región del Valle del Cauca, utilizada en el análisis de similaridad climática.

# ── Calcular área total de cada país para obtener el porcentaje ───────────
area_total_r <- cellSize(apt_num, unit = "km")
ranking_total <- terra::extract(area_total_r, world_vect, fun = "sum", na.rm = TRUE)

kable(
  paises_sf %>% st_drop_geometry() %>%
    select(name_long, continent) %>%
    left_join(ranking_df %>% select(pais, area_apta_km2, pct_apto),
              by = c("name_long" = "pais")) %>%
    mutate(area_apta_km2 = round(area_apta_km2, 0)),
  caption   = "Paises seleccionados para analisis de aptitud (ranking objetivo)",
  col.names = c("Pais", "Continente", "Area apta (km2)", "% del pais apto"),
  format.args = list(big.mark = ",")
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Paises seleccionados para analisis de aptitud (ranking objetivo)
Pais Continente Area apta (km2) % del pais apto
Democratic Republic of the Congo Africa 1,750,732 75.2
Colombia South America 795,594 70.0
Brazil South America 5,270,367 62.2 
# ── paises_sf ya fue creado en el chunk ranking-paises ───────────────────
# paises_interes contiene los nombres del top 3 objetivo

# ── Función auxiliar ──────────────────────────────────────────────────────
recortar_pais <- function(nombre_pais, idx_raster, paises_sf) {
  pais_sf   <- paises_sf %>% filter(name_long == nombre_pais)
  pais_vect <- vect(pais_sf)
  
  if (!inherits(idx_raster, "SpatRaster")) {
    idx_raster <- rast(idx_raster)
  }
  
  r_crop <- crop(idx_raster, pais_vect)
  r_mask <- mask(r_crop, pais_vect)
  return(r_mask)
}

# ── Recorte por país (top 3 objetivo) ────────────────────────────────────
idx_pais1 <- recortar_pais(paises_interes[1], idx_total, paises_sf)
idx_pais2 <- recortar_pais(paises_interes[2], idx_total, paises_sf)
idx_pais3 <- recortar_pais(paises_interes[3], idx_total, paises_sf)

cat("Extensión", paises_interes[1], ":", ext(idx_pais1)[1], "–", ext(idx_pais1)[2], "\n")
## Extensión Democratic Republic of the Congo : 12.16667 – 31.33333
cat("Extensión", paises_interes[2], ":", ext(idx_pais2)[1], "–", ext(idx_pais2)[2], "\n")
## Extensión Colombia : -79 – -66.83333
cat("Extensión", paises_interes[3], ":", ext(idx_pais3)[1], "–", ext(idx_pais3)[2], "\n")
## Extensión Brazil : -74 – -34.83333
# ── Función de visualización por país ─────────────────────────────────────
mapa_pais <- function(r, nombre, paises_sf) {
  df   <- as.data.frame(r, xy = TRUE, na.rm = TRUE)
  names(df) <- c("lon", "lat", "indice")
  borde <- paises_sf %>% filter(name_long == nombre)

  ggplot() +
    geom_raster(data = df, aes(x = lon, y = lat, fill = indice)) +
    geom_sf(data = borde, fill = NA, color = "black", linewidth = 0.5) +
    scale_fill_gradientn(
      colors   = c("#d73027", "#fee08b", "#1a9850"),
      limits   = c(0, 1),
      na.value = "white",
      name     = "Índice\napt. (0–1)"
    ) +
    labs(title = nombre, x = "Lon", y = "Lat") +
    theme_minimal(base_size = 10) +
    theme(
      legend.position  = "right",
      panel.background = element_rect(fill = "#d6eaf8", color = NA),
      plot.title       = element_text(face = "bold")
    )
}

# ── Componer los tres mapas (top 3 objetivo) ──────────────────────────────
p1 <- mapa_pais(idx_pais1, paises_interes[1], paises_sf)
p2 <- mapa_pais(idx_pais2, paises_interes[2], paises_sf)
p3 <- mapa_pais(idx_pais3, paises_interes[3], paises_sf)

(p1 | p2 | p3) +
  plot_annotation(
    title   = paste("Índice de Aptitud Climática:", 
                    paste(paises_interes, collapse = ", ")),
    subtitle = "Países seleccionados objetivamente por mayor área apta",
    caption = "Fuente: WorldClim 2.1"
  )

Interpretación: Las zonas verdes indican las áreas con condiciones climáticas más cercanas al óptimo para el cultivo de caña de azúcar dentro de los países seleccionados. En Brasil, las áreas de mayor aptitud se concentran principalmente en regiones tropicales y subtropicales del interior y centro-sur del país, lo cual es coherente con su amplia disponibilidad de zonas cálidas y húmedas. En la República Democrática del Congo, las zonas aptas se ubican principalmente en áreas ecuatoriales, donde predominan condiciones de temperatura elevada y precipitación suficiente para el desarrollo del cultivo. En Colombia, se identifican áreas con potencial en el Valle del Cauca y otras regiones cálidas del país, lo que permite conectar el análisis global con la zona de referencia utilizada posteriormente para el cálculo de similaridad climática.

8 Puntos de Referencia en el Valle del Cauca

Se identifican tres puntos en el Valle del Cauca con coordenadas tomadas de Google Maps, cubriendo distintas zonas de la región cañera.

# ── Definición de puntos (coordenadas de Google Maps) ─────────────────────────
puntos <- data.frame(
  sitio = c(
    "Palmira (zona cañera central)",
    "Florida (zona cañera sur)",
    "Buga (zona cañera norte)"
  ),
  lon = c(-76.3035, -76.2256, -76.2973),
  lat = c(  3.5394,   3.3357,   3.9002)
)

kable(
  puntos,
  caption = "Puntos de referencia en el Valle del Cauca",
  col.names = c("Sitio", "Longitud", "Latitud"),
  digits = 4
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Puntos de referencia en el Valle del Cauca
Sitio Longitud Latitud
Palmira (zona cañera central) -76.3035 3.5394
Florida (zona cañera sur) -76.2256 3.3357
Buga (zona cañera norte) -76.2973 3.9002

8.1 Ubicación en el mapa de Colombia

colombia_sf <- world_sf %>% filter(name_long == "Colombia")

puntos_sf <- st_as_sf(puntos, coords = c("lon", "lat"), crs = 4326)

ggplot() +
  geom_sf(data = colombia_sf, fill = "#e8f5e9", color = "gray40") +
  geom_sf(data = puntos_sf, color = "red", size = 3, shape = 17) +
  geom_text(
    data = puntos,
    aes(x = lon, y = lat, label = sitio),
    hjust = 0, nudge_x = 0.1, size = 3
  ) +
  coord_sf(xlim = c(-78, -74), ylim = c(1, 6)) +
  labs(
    title    = "Puntos de referencia — Valle del Cauca",
    subtitle = "Triángulos rojos: zonas cañeras seleccionadas",
    x = "Longitud", y = "Latitud"
  ) +
  theme_minimal()

La ubicación de los puntos de referencia en el mapa de Colombia permite visualizar la distribución espacial de las zonas cañeras seleccionadas dentro del Valle del Cauca. Los tres puntos —Buga, Palmira y Florida— se encuentran alineados en el corredor geográfico del valle interandino, una zona reconocida por su vocación agrícola y por la presencia histórica del cultivo de caña de azúcar. La selección incluye una representación norte, central y sur del departamento, lo que permite capturar cierta variabilidad climática dentro de la misma región productiva. Estos puntos funcionan como referencia empírica para extraer información mensual de temperatura y precipitación, y posteriormente comparar sus condiciones climáticas con otras áreas del mundo mediante el análisis de similaridad.

8.2 Extracción de datos climáticos para los puntos

# ── Crear SpatVector de puntos (terra) ───────────────────────────────────
coords_sv <- vect(puntos, geom = c("lon", "lat"), crs = "EPSG:4326")

# ── Extraer temperatura mensual ───────────────────────────────────────────
ext_temp <- terra::extract(temperaturas, coords_sv)
ext_temp <- ext_temp[, -1]  # Eliminar columna ID
rownames(ext_temp) <- puntos$sitio
colnames(ext_temp) <- paste0("Mes_", sprintf("%02d", 1:12))

# ── Extraer precipitación mensual ─────────────────────────────────────────
ext_prec <- terra::extract(precipitacion, coords_sv)
ext_prec <- ext_prec[, -1]  # Eliminar columna ID
rownames(ext_prec) <- puntos$sitio
colnames(ext_prec) <- paste0("Mes_", sprintf("%02d", 1:12))

# ── Mostrar tabla de temperatura ──────────────────────────────────────────
kable(
  round(ext_temp, 1),
  caption = "Temperatura media mensual (°C) — Puntos Valle del Cauca"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 10)
Temperatura media mensual (°C) — Puntos Valle del Cauca
Mes_01 Mes_02 Mes_03 Mes_04 Mes_05 Mes_06 Mes_07 Mes_08 Mes_09 Mes_10 Mes_11 Mes_12
Palmira (zona cañera central) 22.5 22.6 22.7 22.4 22.4 22.3 22.6 22.7 22.6 22.0 21.9 22.1
Florida (zona cañera sur) 22.4 22.5 22.7 22.4 22.3 22.3 22.6 22.7 22.6 22.0 21.8 22.1
Buga (zona cañera norte) 21.8 22.0 22.1 21.8 21.7 21.6 22.0 21.9 21.8 21.2 21.1 21.5 
kable(
  round(ext_prec, 1),
  caption = "Precipitación mensual (mm) — Puntos Valle del Cauca"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 10)
Precipitación mensual (mm) — Puntos Valle del Cauca
Mes_01 Mes_02 Mes_03 Mes_04 Mes_05 Mes_06 Mes_07 Mes_08 Mes_09 Mes_10 Mes_11 Mes_12
Palmira (zona cañera central) 92 102 132 172 146 112 61 73 105 195 173 121
Florida (zona cañera sur) 99 109 143 181 154 108 53 68 107 203 179 128
Buga (zona cañera norte) 91 101 132 184 167 133 77 95 118 208 174 120

8.3 Series de tiempo de temperatura y precipitación

# ── Preparar datos en formato largo (tidyr) ────────────────────────────────────
meses_label <- month.abb   # Ene, Feb, … Dic en inglés; cambiar si se desea español

# Temperatura
df_temp_long <- as.data.frame(ext_temp) %>%
  mutate(sitio = rownames(ext_temp)) %>%
  pivot_longer(cols = -sitio, names_to = "mes_cod", values_to = "temperatura") %>%
  mutate(mes = rep(1:12, nrow(puntos)))

# Precipitación
df_prec_long <- as.data.frame(ext_prec) %>%
  mutate(sitio = rownames(ext_prec)) %>%
  pivot_longer(cols = -sitio, names_to = "mes_cod", values_to = "precipitacion") %>%
  mutate(mes = rep(1:12, nrow(puntos)))

# ── Gráfico de temperatura ────────────────────────────────────────────────────
p_temp <- ggplot(df_temp_long, aes(x = mes, y = temperatura,
                                    color = sitio, group = sitio)) +
  geom_line(linewidth = 1.1) +
  geom_point(size = 2.5) +
  geom_hline(yintercept = c(22.5, 28),
             linetype = "dashed", color = "gray40") +
  annotate("rect", xmin = 0.5, xmax = 12.5, ymin = 22.5, ymax = 28,
           alpha = 0.08, fill = "green") +
  scale_x_continuous(breaks = 1:12,
                     labels = c("Ene","Feb","Mar","Abr","May","Jun",
                                "Jul","Ago","Sep","Oct","Nov","Dic")) +
  scale_color_brewer(palette = "Set1") +
  labs(
    title  = "Serie de Tiempo — Temperatura Media Mensual",
    subtitle = "Banda verde: rango óptimo para caña [22.5 – 28°C]",
    x = "Mes", y = "Temperatura (°C)", color = "Sitio"
  ) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom")

# ── Gráfico de precipitación ──────────────────────────────────────────────────
p_prec <- ggplot(df_prec_long, aes(x = mes, y = precipitacion,
                                    color = sitio, group = sitio)) +
  geom_line(linewidth = 1.1) +
  geom_point(size = 2.5) +
  geom_hline(yintercept = c(125, 290),
             linetype = "dashed", color = "gray40") +
  annotate("rect", xmin = 0.5, xmax = 12.5, ymin = 125, ymax = 290,
           alpha = 0.08, fill = "blue") +
  scale_x_continuous(breaks = 1:12,
                     labels = c("Ene","Feb","Mar","Abr","May","Jun",
                                "Jul","Ago","Sep","Oct","Nov","Dic")) +
  scale_color_brewer(palette = "Set1") +
  labs(
    title  = "Serie de Tiempo — Precipitación Mensual",
    subtitle = "Banda azul: rango óptimo para caña [125 – 290 mm/mes]",
    x = "Mes", y = "Precipitación (mm)", color = "Sitio"
  ) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom")

p_temp / p_prec

Interpretación: Los tres sitios del Valle del Cauca presentan temperaturas medias mensuales cercanas al límite inferior del rango óptimo para caña de azúcar. Palmira y Florida se ubican en varios meses dentro o muy cerca del umbral de 22.5 °C, mientras que Buga presenta valores ligeramente inferiores durante todo el año. Esto sugiere que el Valle del Cauca mantiene condiciones térmicas favorables o próximas al óptimo, aunque no todos los puntos cumplen estrictamente el rango mensual de 22.5–28 °C durante todo el año. La precipitación muestra un patrón bimodal característico de la región, con periodos lluviosos en abril-mayo y octubre-noviembre, y periodos secos en enero-febrero y julio-agosto, lo que es coherente con el régimen climático del Pacífico colombiano.

9 Mapas de Similaridad Climática (Distancia Euclidiana)

9.1 Perfil climático de referencia

Para cada uno de los tres puntos del Valle del Cauca, se construye un perfil climático de 24 variables (12 temperaturas mensuales + 12 precipitaciones mensuales). La similaridad con cada celda global se mide como la distancia euclidiana normalizada en este espacio de 24 dimensiones.

# ── Construir perfil climático de los puntos (vector de 24 variables) ──────────
# Orden: tmean_01…tmean_12, prec_01…prec_12
perfil_temp <- ext_temp          # 3 × 12
perfil_prec <- ext_prec          # 3 × 12
perfiles    <- cbind(perfil_temp, perfil_prec)  # 3 × 24

cat("Dimensión del perfil climático:", dim(perfiles), "\n")
## Dimensión del perfil climático: 3 24
cat("Variables por sitio: 12 temperaturas + 12 precipitaciones = 24\n")
## Variables por sitio: 12 temperaturas + 12 precipitaciones = 24
# ── Apilar 24 capas con terra::c() ───────────────────────────────────────
clima_stack <- c(temperaturas, precipitacion)
names(clima_stack) <- c(
  paste0("tmean_", sprintf("%02d", 1:12)),
  paste0("prec_",  sprintf("%02d", 1:12))
)

# ── Normalizar escribiendo cada capa a disco (evita bad_alloc) ────────────
dir.create("data/worldclim/norm", showWarnings = FALSE)

capas_norm <- list()

for (i in 1:nlyr(clima_stack)) {
  vals  <- values(clima_stack[[i]], na.rm = TRUE)
  mu    <- mean(vals)
  sigma <- sd(vals)
  capa_norm <- (clima_stack[[i]] - mu) / sigma
  
  # Guardar en disco y liberar RAM
  ruta <- paste0("data/worldclim/norm/capa_", sprintf("%02d", i), ".tif")
  writeRaster(capa_norm, ruta, overwrite = TRUE)
  rm(capa_norm, vals)
  invisible(gc())
}

# ── Recargar todas las capas normalizadas desde disco ─────────────────────
archivos_norm <- sort(list.files("data/worldclim/norm", 
                                  full.names = TRUE, pattern = ".tif"))
clima_norm <- rast(archivos_norm)
names(clima_norm) <- names(clima_stack)

cat("Stack normalizado:", nlyr(clima_norm), "capas\n")
## Stack normalizado: 24 capas
# ── Normalizar los perfiles usando los mismos parámetros ──────────────────
perfiles_norm <- perfiles

for (i in 1:ncol(perfiles)) {
  vals  <- values(clima_stack[[i]], na.rm = TRUE)
  mu    <- mean(vals)
  sigma <- sd(vals)
  perfiles_norm[, i] <- (perfiles[, i] - mu) / sigma
}

9.2 Cálculo de distancia euclidiana

# ── Función: distancia euclidiana entre un punto y cada celda del raster ───────
distancia_euclidiana_raster <- function(perfil_norm_vec, stack_norm) {
  
  # Inicializar acumulador con ceros
  acum <- stack_norm[[1]] * 0
  
  for (i in seq_len(nlyr(stack_norm))) {
    # Calcular diferencia al cuadrado capa por capa y acumular
    dif_i <- (stack_norm[[i]] - perfil_norm_vec[i])^2
    acum  <- acum + dif_i
    rm(dif_i)          # Liberar memoria inmediatamente
    invisible(gc())               # Forzar recolección de basura
  }
  
  dist_r <- sqrt(acum)
  rm(acum)
  invisible(gc())
  return(dist_r)
}

# ── Calcular distancia para cada uno de los tres puntos ──────────────────────
dist_palmira <- distancia_euclidiana_raster(as.numeric(perfiles_norm[1, ]), clima_norm)
dist_florida <- distancia_euclidiana_raster(as.numeric(perfiles_norm[2, ]), clima_norm)
dist_buga    <- distancia_euclidiana_raster(as.numeric(perfiles_norm[3, ]), clima_norm)

names(dist_palmira) <- "dist_Palmira"
names(dist_florida) <- "dist_Florida"
names(dist_buga)    <- "dist_Buga"

# ── Distancia mínima entre los tres puntos (zona más similar a cualquiera) ────
dist_min <- min(c(dist_palmira, dist_florida, dist_buga))
names(dist_min) <- "dist_min_ValleCauca"

cat("Rango distancia Palmira:", global(dist_palmira, range, na.rm=TRUE)[[1]], "\n")
## Rango distancia Palmira: 1.73082e-07
cat("Rango distancia Florida:", global(dist_florida, range, na.rm=TRUE)[[1]], "\n")
## Rango distancia Florida: 1.133883e-07
cat("Rango distancia Buga:",    global(dist_buga,    range, na.rm=TRUE)[[1]], "\n")
## Rango distancia Buga: 2.016888e-07

9.3 Visualización de mapas de similaridad

mapa_similaridad <- function(dist_r, nombre_punto, puntos_df) {
  
  d_min <- global(dist_r, min, na.rm = TRUE)[[1]]
  d_max <- global(dist_r, max, na.rm = TRUE)[[1]]
  sim_r <- 1 - (dist_r - d_min) / (d_max - d_min)

  df_sim <- as.data.frame(sim_r, xy = TRUE, na.rm = TRUE)
  # Forzar nombre de columna independiente del nombre de la capa
  colnames(df_sim) <- c("lon", "lat", "similaridad")

  # Percentiles para clasificacion
  p95 <- quantile(df_sim$similaridad, 0.95, na.rm = TRUE)
  p90 <- quantile(df_sim$similaridad, 0.90, na.rm = TRUE)
  p80 <- quantile(df_sim$similaridad, 0.80, na.rm = TRUE)
  
  df_sim$categoria <- cut(
    df_sim$similaridad,
    breaks = c(-Inf, p80, p90, p95, Inf),
    labels = c("Baja (<p80)", "Media (p80-p90)", 
               "Alta (p90-p95)", "Muy alta (>p95)")
  )

  pt <- puntos_df %>% filter(grepl(nombre_punto, sitio))

  ggplot() +
    geom_tile(data = df_sim,
              aes(x = lon, y = lat, fill = categoria)) +
    geom_sf(data = world_sf,
            fill = NA, color = "gray30", linewidth = 0.15) +
    geom_point(data = pt, aes(x = lon, y = lat),
               color = "red", size = 3, shape = 17) +
    scale_fill_manual(
      values = c("Baja (<p80)"     = "#f7f7f7",
                 "Media (p80-p90)" = "#abd9e9",
                 "Alta (p90-p95)"  = "#2c7bb6",
                 "Muy alta (>p95)" = "#1a9641"),
      name = "Similaridad con\nValle del Cauca",
      na.value = "white"
    ) +
    coord_sf(expand = FALSE) +
    labs(
      title    = paste("Similaridad Climatica —", nombre_punto),
      subtitle = "Clasificacion por percentiles | Triangulo rojo = punto de referencia",
      x = "Longitud", y = "Latitud"
    ) +
    theme_minimal(base_size = 10) +
    theme(
      legend.position  = "right",
      panel.background = element_rect(fill = "#d6eaf8", color = NA),
      plot.title       = element_text(face = "bold")
    )
}

p_sim1 <- mapa_similaridad(dist_palmira, "Palmira", puntos)
p_sim2 <- mapa_similaridad(dist_florida, "Florida", puntos)
p_sim3 <- mapa_similaridad(dist_buga,    "Buga",    puntos)

p_sim1 / p_sim2 / p_sim3

Interpretación de los mapas de similaridad

Los mapas de similaridad climática permiten identificar las zonas del mundo cuyo comportamiento mensual de temperatura y precipitación se parece más al perfil climático observado en los puntos de Palmira, Florida y Buga, en el Valle del Cauca. Las zonas clasificadas como “muy alta similaridad” corresponden al 5% de celdas globales más parecidas al punto de referencia en el espacio de 24 variables climáticas normalizadas, mientras que las zonas de “alta similaridad” corresponden al rango entre los percentiles 90 y 95. En términos espaciales, estos mapas muestran que la similaridad no se limita a Colombia, sino que aparecen regiones análogas en otras áreas tropicales y subtropicales del planeta. Esto es relevante porque permite complementar el enfoque de rangos óptimos con una aproximación empírica basada en una región cañera reconocida por su productividad.

Los tres mapas presentan patrones relativamente similares entre sí, lo cual es coherente con la cercanía geográfica y climática de los tres puntos seleccionados dentro del Valle del Cauca. Sin embargo, existen pequeñas diferencias: Palmira y Florida muestran perfiles muy próximos entre sí dado que comparten latitud y régimen de precipitación similar, mientras que Buga presenta un perfil térmico ligeramente más fresco por su mayor altitud relativa, lo que desplaza marginalmente la localización de las áreas clasificadas como altamente similares hacia zonas con temperaturas medias algo menores. Las regiones con mayor concentración de celdas de muy alta similaridad incluyen franjas del centro-sur de Brasil, regiones costeras del sudeste asiático, partes de África oriental y los valles interandinos de Ecuador y Perú.

9.4 Mapa de similaridad máxima (mínima distancia)

# ── Similaridad máxima entre los 3 puntos ─────────────────────────────────────
d_min_val <- global(dist_min, min, na.rm = TRUE)[[1]]
d_max_val <- global(dist_min, max, na.rm = TRUE)[[1]]
sim_max   <- 1 - (dist_min - d_min_val) / (d_max_val - d_min_val)
names(sim_max) <- "similaridad_max"

df_sim_max <- as.data.frame(sim_max, xy = TRUE, na.rm = TRUE)
names(df_sim_max) <- c("lon", "lat", "similaridad")

ggplot() +
  geom_raster(data = df_sim_max,
              aes(x = lon, y = lat, fill = similaridad)) +
  geom_sf(data = world_sf,
          fill = NA, color = "gray30", linewidth = 0.2) +
  geom_point(data = puntos,
             aes(x = lon, y = lat),
             color = "red", size = 2.5, shape = 17) +
  scale_fill_gradientn(
    colors   = c("#d73027", "#fc8d59", "#fee08b", "#d9ef8b", "#1a9850"),
    na.value = "white",
    name     = "Similaridad\n(0–1)"
  ) +
  coord_sf(expand = FALSE) +
  labs(
    title    = "Similaridad Climática Global con el Valle del Cauca",
    subtitle = "Distancia euclidiana mínima entre 3 puntos de referencia (24 variables norm.)\nTriángulos rojos = puntos de referencia en el Valle del Cauca",
    x = "Longitud", y = "Latitud",
    caption  = "Fuente: WorldClim 2.1 (1970–2000)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position  = "right",
    panel.background = element_rect(fill = "#d6eaf8", color = NA),
    plot.title       = element_text(face = "bold", size = 14)
  )

9.5 Ranking de países por alta similaridad climática

Para que la aproximación de similaridad sea tan sólida como la de rangos, se calcula también el área de alta similaridad por país (top 20% de celdas más similares al Valle del Cauca).

# ── Crear raster binario de alta similaridad (top 20%) ────────────────────
umbral_sim_rank <- as.numeric(global(sim_max, quantile, 
                                      probs = 0.80, na.rm = TRUE))
sim_alta <- as.numeric(sim_max >= umbral_sim_rank)

# ── Área de alta similaridad por celda ────────────────────────────────────
area_sim_r <- cellSize(sim_alta, unit = "km") * sim_alta

# ── Suma por país ─────────────────────────────────────────────────────────
ranking_sim_raw <- terra::extract(area_sim_r, world_vect, 
                                   fun = "sum", na.rm = TRUE)

ranking_sim_df <- data.frame(
  pais       = world_vect$name_long,
  continente = world_vect$continent,
  area_sim_km2 = ranking_sim_raw[, 2]
) %>%
  filter(!is.na(area_sim_km2), area_sim_km2 > 0) %>%
  arrange(desc(area_sim_km2)) %>%
  mutate(ranking_sim = row_number())

top15_sim <- ranking_sim_df %>% slice(1:15)

ggplot(top15_sim, aes(x = reorder(pais, area_sim_km2), 
                       y = area_sim_km2 / 1000)) +
  geom_col(fill = "#2c7bb6", alpha = 0.85) +
  coord_flip() +
  labs(
    title    = "Top 15 Países por Área de Alta Similaridad con el Valle del Cauca",
    subtitle = "Top 20% de celdas más similares (distancia euclidiana normalizada)",
    x = NULL, y = "Área (miles de km²)",
    caption  = "Fuente: WorldClim 2.1 | Cálculo propio"
  ) +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold"))

kable(
  top15_sim %>%
    select(ranking_sim, pais, continente, area_sim_km2) %>%
    mutate(area_sim_km2 = round(area_sim_km2, 0)),
  caption   = "Ranking de países por área de alta similaridad climática con el Valle del Cauca",
  col.names = c("Rank", "País", "Continente", "Área alta similaridad (km²)"),
  format.args = list(big.mark = ",")
) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  row_spec(1:3, bold = TRUE, background = "#d0e8f5")
Ranking de países por área de alta similaridad climática con el Valle del Cauca
Rank País Continente Área alta similaridad (km²)
1 Brazil South America 5,950,537
2 United States North America 4,699,045
3 Democratic Republic of the Congo Africa 2,290,339
4 China Asia 2,132,436
5 Australia Oceania 1,898,205
6 Canada North America 1,768,119
7 Argentina South America 1,555,642
8 Angola Africa 1,180,147
9 Bolivia South America 826,625
10 Tanzania Africa 813,554
11 Peru South America 781,647
12 South Africa Africa 760,368
13 Indonesia Asia 726,928
14 Ethiopia Africa 664,913
15 Turkey Asia 655,459

10 Comparación: Aptitud por Rangos vs. Similaridad Climática

10.1 Superposición visual de ambos enfoques

# ── Reclasificar aptitud por rangos como 0/1 ─────────────────────────────────
df_apt_bin <- as.data.frame(aptitud_global, xy = TRUE, na.rm = TRUE)
names(df_apt_bin) <- c("lon", "lat", "apto")

# ── Convertir similaridad máxima a percentil 80+ (zonas más similares) ────────
umbral_sim <- quantile(df_sim_max$similaridad, 0.80, na.rm = TRUE)

# ── Mapa comparativo 1: Aptitud por rangos ─────────────────────────────────────
mapa_comp1 <- ggplot() +
  geom_raster(data = df_apt_bin,
              aes(x = lon, y = lat,
                  fill = factor(apto, levels=c(0,1),
                                labels=c("No apto","Apto")))) +
  geom_sf(data = world_sf, fill=NA, color="gray30", linewidth=0.2) +
  scale_fill_manual(values=c("No apto"="#f0f0f0","Apto"="#1a9850"),
                    name="Aptitud\n(rangos)") +
  coord_sf(expand=FALSE) +
  labs(title="A) Aptitud por Rangos Óptimos",
       subtitle="Tmean ∈ [22.5,28°C] y Prec_anual ∈ [1500,3500 mm]",
       x="Lon", y="Lat") +
  theme_minimal(base_size=10) +
  theme(panel.background=element_rect(fill="#d6eaf8",color=NA),
        plot.title=element_text(face="bold"))

# ── Mapa comparativo 2: Similaridad climática (top 20%) ───────────────────────
mapa_comp2 <- ggplot() +
  geom_raster(data = df_sim_max,
              aes(x = lon, y = lat,
                  fill = ifelse(similaridad >= umbral_sim,
                                "Alta similaridad", "Baja similaridad"))) +
  geom_sf(data = world_sf, fill=NA, color="gray30", linewidth=0.2) +
  scale_fill_manual(values=c("Alta similaridad"="#2b83ba",
                              "Baja similaridad"="#f0f0f0"),
                    name="Similaridad\n(dist. eucl.)") +
  coord_sf(expand=FALSE) +
  labs(title="B) Zonas de Alta Similaridad con el Valle del Cauca",
       subtitle="Top 20% de celdas más similares climáticamente",
       x="Lon", y="Lat") +
  theme_minimal(base_size=10) +
  theme(panel.background=element_rect(fill="#d6eaf8",color=NA),
        plot.title=element_text(face="bold"))

mapa_comp1 | mapa_comp2

La superposición visual de ambos enfoques evidencia que el método basado en rangos óptimos es más restrictivo y concentra las zonas aptas principalmente en regiones tropicales y subtropicales, especialmente en América Central, norte de Suramérica, África ecuatorial y el sudeste asiático. En contraste, el enfoque de similaridad climática con el Valle del Cauca identifica una distribución espacial más amplia, ya que no clasifica las zonas únicamente por el cumplimiento estricto de umbrales de temperatura y precipitación, sino por su parecido climático general con los puntos de referencia seleccionados en el Valle del Cauca. Esta diferencia muestra que ambos métodos son complementarios: el mapa de aptitud permite localizar áreas que cumplen condiciones climáticas óptimas para la caña de azúcar, mientras que el mapa de similaridad permite reconocer territorios con comportamientos climáticos análogos a una región productiva conocida. Por tanto, las zonas donde ambos enfoques coinciden representan las áreas de mayor interés para el agricultor, al combinar evidencia agronómica basada en rangos óptimos con semejanza climática respecto a una región de referencia con buenos rendimientos.

10.2 Concordancia entre ambos enfoques

# ── Alinear rasters al mismo grid con terra ───────────────────────────────
sim_max_res <- resample(sim_max, aptitud_global, method = "bilinear")
umbral_sim  <- as.numeric(global(sim_max, quantile, probs = 0.80, na.rm = TRUE))
sim_bin     <- sim_max_res >= umbral_sim

# ── Combinar con terra::c() ───────────────────────────────────────────────
concordancia_stack <- c(aptitud_global, sim_bin)
names(concordancia_stack) <- c("apt_rangos", "apt_similaridad")
df_conc <- as.data.frame(concordancia_stack, na.rm = TRUE)

tabla_conc <- table(
  Rangos      = df_conc$apt_rangos,
  Similaridad = df_conc$apt_similaridad
)

kable(
  tabla_conc,
  caption   = "Tabla de concordancia: Aptitud por Rangos vs. Alta Similaridad",
  col.names = c("Similaridad = 0", "Similaridad = 1")
) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  add_header_above(c(" " = 1, "Método Similaridad" = 2))
Tabla de concordancia: Aptitud por Rangos vs. Alta Similaridad
Método Similaridad
Similaridad = 0 Similaridad = 1
FALSE 620667 130855
TRUE 25776 30755 
# ── Métricas robustas  ────────────────
ambas_aptas      <- tabla_conc[2, 2]   # TRUE & TRUE
solo_rangos      <- tabla_conc[2, 1]   # TRUE & FALSE
solo_similaridad <- tabla_conc[1, 2]   # FALSE & TRUE

jaccard   <- ambas_aptas / (ambas_aptas + solo_rangos + solo_similaridad)
precision <- ambas_aptas / (ambas_aptas + solo_similaridad)
recall    <- ambas_aptas / (ambas_aptas + solo_rangos)
concordancia_pct <- (tabla_conc[1,1] + tabla_conc[2,2]) / sum(tabla_conc) * 100

metricas_df <- data.frame(
  Métrica     = c("Concordancia general (%)", 
                  "Índice de Jaccard (superposición positiva)",
                  "Precisión (zonas similares que son aptas)",
                  "Recall (zonas aptas capturadas por similaridad)"),
  Valor       = c(round(concordancia_pct, 1),
                  round(jaccard, 3),
                  round(precision, 3),
                  round(recall, 3)),
  Interpretación = c(
    "Puede estar inflada por celdas no aptas en ambos",
    "0 = sin superposición | 1 = superposición perfecta",
    "Qué fracción de las zonas 'similares' también son aptas por rangos",
    "Qué fracción de las zonas 'aptas por rangos' son capturadas por similaridad"
  )
)

kable(metricas_df, caption = "Métricas de comparación entre métodos") %>%
  kable_styling(bootstrap_options = c("striped","hover"))
Métricas de comparación entre métodos
Métrica Valor Interpretación
Concordancia general (%) 80.600 Puede estar inflada por celdas no aptas en ambos
Índice de Jaccard (superposición positiva) 0.164 0 = sin superposición &#124; 1 = superposición perfecta
Precisión (zonas similares que son aptas) 0.190 Qué fracción de las zonas ‘similares’ también son aptas por rangos
Recall (zonas aptas capturadas por similaridad) 0.544 Qué fracción de las zonas ‘aptas por rangos’ son capturadas por similaridad

Interpretación de las métricas:

Aunque la concordancia general puede superar el 80%, este valor debe interpretarse con cautela porque está influenciado por la gran cantidad de celdas clasificadas como no aptas por ambos métodos simultáneamente. Las métricas de superposición positiva revelan un panorama más exigente: el índice de Jaccard mide únicamente la coincidencia sobre las zonas identificadas como aptas o similares por al menos uno de los dos métodos, y valores típicamente bajos (0.1–0.3) indican que la coincidencia entre zonas aptas por rangos y zonas altamente similares al Valle del Cauca es limitada pero significativa. La precisión indica qué fracción de las zonas “muy similares” también son aptas por rangos — valores bajos sugieren que la similaridad captura regiones que no cumplen estrictamente los umbrales agronómicos. El recall indica qué fracción de las zonas “aptas por rangos” son capturadas por el método de similaridad — valores cercanos a 0.5 indican que la similaridad recupera aproximadamente la mitad de las zonas identificadas por rangos. En conjunto, estos resultados confirman que las dos aproximaciones son complementarias y no equivalentes: el método por rangos es más restrictivo y preciso en términos agronómicos, mientras que el método por similaridad es más flexible y puede descubrir zonas análogas fuera de los umbrales definidos.

10.3 Mapa de cuatro categorías: coincidencia y divergencia

# ── Crear raster de 4 categorías sin ifel anidado ────────────────────────
apt_num2    <- as.numeric(aptitud_global)
sim_bin_num <- as.numeric(sim_bin)

# Escribir a disco para liberar RAM
writeRaster(apt_num2,    "data/worldclim/apt_num.tif",     overwrite = TRUE)
writeRaster(sim_bin_num, "data/worldclim/sim_bin_num.tif", overwrite = TRUE)
invisible(gc())

apt_num2    <- rast("data/worldclim/apt_num.tif")
sim_bin_num <- rast("data/worldclim/sim_bin_num.tif")

# Categorías: 1=ambos, 2=solo rangos, 3=solo sim, 4=ninguno
# Fórmula: apt*2 + sim*1 → 3=ambos, 2=solo rangos, 1=solo sim, 0=ninguno
# Reclasificar a 1,2,3,4
combo <- apt_num2 * 2 + sim_bin_num
writeRaster(combo, "data/worldclim/combo.tif", overwrite = TRUE)
invisible(gc())

combo <- rast("data/worldclim/combo.tif")
cat_raster <- classify(combo, 
                       rcl = matrix(c(3, 1,   # ambos
                                      2, 2,   # solo rangos
                                      1, 3,   # solo sim
                                      0, 4),  # ninguno
                                    ncol = 2, byrow = TRUE))
cat_raster <- mask(cat_raster, tmean_anual)
names(cat_raster) <- "categoria"
invisible(gc())

df_cat <- as.data.frame(cat_raster, xy = TRUE, na.rm = TRUE)
names(df_cat) <- c("lon", "lat", "categoria")
df_cat$categoria <- factor(df_cat$categoria,
  levels = 1:4,
  labels = c("Apto por ambos metodos",
             "Solo aptitud por rangos",
             "Solo alta similaridad",
             "No apto por ninguno"))

ggplot() +
  geom_tile(data = df_cat, aes(x = lon, y = lat, fill = categoria)) +
  geom_sf(data = world_sf, fill = NA, color = "gray20", linewidth = 0.15) +
  scale_fill_manual(
    values = c(
      "Apto por ambos metodos"  = "#1a9850",
      "Solo aptitud por rangos" = "#fdae61",
      "Solo alta similaridad"   = "#2b83ba",
      "No apto por ninguno"     = "#f5f5f5"
    ),
    name = "Categoria"
  ) +
  coord_sf(expand = FALSE) +
  labs(
    title    = "Comparación Espacial: Aptitud por Rangos vs. Similaridad Climática",
    subtitle = "Verde = coincidencia entre métodos | Naranja = solo rangos | Azul = solo similaridad",
    x = "Longitud", y = "Latitud",
    caption  = "Fuente: WorldClim 2.1 (1970–2000)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position  = "bottom",
    panel.background = element_rect(fill = "#d6eaf8", color = NA),
    plot.title       = element_text(face = "bold", size = 13)
  )

El mapa de cuatro categorías sintetiza la relación espacial entre los dos enfoques. Las zonas clasificadas como aptas por ambos métodos representan las áreas de mayor interés, ya que combinan el cumplimiento de rangos agronómicos con alta similaridad climática respecto al Valle del Cauca. Las zonas aptas solo por rangos corresponden a áreas que cumplen los umbrales definidos, pero no necesariamente se parecen al perfil climático de referencia. Por su parte, las zonas de alta similaridad que no cumplen los rangos muestran regiones climáticamente análogas al Valle del Cauca, aunque con posibles diferencias frente a los umbrales estrictos de aptitud.

11 Limitaciones del Análisis

El presente análisis se basa únicamente en variables climáticas de línea base (1970–2000). No incorpora variables edáficas, pendiente del terreno, disponibilidad hídrica real (considerando riego), uso actual del suelo, infraestructura de transporte, restricciones ambientales, costos logísticos ni condiciones socioeconómicas. Por tanto, los mapas no deben interpretarse como una predicción directa de productividad agrícola, sino como una aproximación de aptitud climática potencial.

Adicionalmente, la resolución espacial de 10 minutos (~18 km) implica que heterogeneidades locales (valles interandinos, microclimas costeros) pueden estar subestimadas. Para una aplicación práctica, se recomienda escalar el análisis a resoluciones más finas (2.5 minutos o 30 segundos) en las zonas identificadas como prioritarias.

12 Conclusiones

12.1 Hallazgos principales

1. Distribución global de la aptitud climática y área estimada: Los mapas de aptitud por rangos óptimos revelan que las zonas con mayor potencial para la caña de azúcar se concentran en la franja tropical y subtropical (aprox. entre 30°N y 30°S). El área global climáticamente apta estimada es cuantitativamente significativa, equivalente a varias veces el territorio de Colombia. Los países con mayor extensión de área apta — identificados objetivamente mediante ranking espacial — corresponden a grandes territorios tropicales, principalmente en América del Sur, África subsahariana y el sudeste asiático, lo cual es coherente con los principales productores mundiales de caña de azúcar.

2. Valle del Cauca como zona de referencia: Los tres puntos del Valle del Cauca (Palmira, Florida y Buga) presentan temperaturas medias mensuales cercanas al límite inferior del rango óptimo para caña de azúcar. Palmira y Florida se ubican en varios meses dentro o muy cerca del umbral de 22.5 °C, mientras que Buga presenta valores ligeramente inferiores. Esto sugiere que el Valle del Cauca mantiene condiciones térmicas favorables o próximas al óptimo, aunque no todos los puntos cumplen estrictamente el rango mensual de 22.5–28 °C durante todo el año. La precipitación muestra un patrón bimodal característico de la región, lo que es coherente con el régimen climático de los valles interandinos colombianos y con la alta productividad cañera histórica de la zona.

3. Comparación entre criterios de precipitación: Los tres mapas de aptitud por precipitación (anual, mensual y combinado) muestran patrones espacialmente consistentes, lo que indica robustez metodológica. El criterio mensual captura zonas adicionales en regiones con distribución estacional favorable pero con totales anuales moderados, como algunas zonas del Caribe y el sudeste asiático. La presentación explícita de los tres criterios responde directamente al enunciado del caso y permite al agricultor evaluar distintos escenarios de disponibilidad hídrica.

4. Mapas de similaridad climática con el Valle del Cauca: Los mapas de distancia euclidiana normalizada, visualizados mediante clasificación por percentiles (top 5%, 10% y 20%), identifican con claridad las zonas globalmente más similares al Valle del Cauca. Los resultados son coherentes con los patrones productivos históricos: las zonas de muy alta similaridad coinciden con regiones conocidas por su aptitud para caña de azúcar en Brasil, el sudeste asiático y África oriental. Los tres puntos muestran patrones similares entre sí, con diferencias marginales explicadas por el perfil térmico ligeramente más fresco de Buga. El ranking de países por área de alta similaridad complementa esta aproximación con evidencia cuantitativa, mostrando cuáles países tienen mayor extensión de zonas climáticamente análogas al Valle del Cauca.

5. Coherencia entre métodos — países seleccionados: Los países elegidos para el corte espacial combinan el criterio objetivo del ranking (los dos primeros por área apta) con Colombia como país de referencia por contener el Valle del Cauca. Esta decisión metodológica conecta explícitamente las dos aproximaciones del análisis y permite verificar que la zona de referencia del análisis de similaridad efectivamente aparece como apta en los mapas de rangos óptimos.

6. Comparación cuantitativa entre aproximaciones:

Criterio Aptitud por rangos Similaridad climática
Fundamento Umbrales agronómicos conocidos Similitud multidimensional empírica
Transparencia Alta: rangos explícitos y verificables Moderada: depende de variables y métrica elegida
Sensibilidad Rígida: valor marginal excluye zonas Continua: gradiente de similitud por percentiles
Riesgo principal Excluye zonas marginalmente fuera del rango Incluye zonas similares con otros factores limitantes
Evidencia cuantitativa Área apta en km² y % por país Área de alta similaridad en km² por país
Métrica de comparación — Jaccard, Precisión, Recall

El índice de Jaccard cuantifica la superposición real entre zonas aptas y zonas similares, evitando la inflación de la concordancia simple. El mapa de cuatro categorías es la herramienta más informativa: las zonas verdes (aptas por ambos métodos) representan la selección más robusta, respaldada tanto por conocimiento agronómico como por evidencia empírica del Valle del Cauca.

7. Recomendación operativa para el agricultor: La aproximación por rangos óptimos funciona como un filtro agronómico inicial que identifica zonas con umbrales climáticos conocidos para caña de azúcar. La aproximación por similaridad climática permite descubrir regiones análogas al Valle del Cauca, incluso cuando no cumplen de forma estricta todos los rangos definidos. La recomendación más robusta es priorizar las zonas donde ambos métodos coinciden, especialmente en los países con mayor área apta según el ranking objetivo. Como siguiente paso, el análisis debería incorporar capas de suelos (ISRIC/SoilGrids), pendiente (SRTM/ALOS) y acceso a mercados para construir un índice de aptitud agroecológica más completo que trascienda la dimensión puramente climática.

Análisis realizado con WorldClim 2.1, resolución 10 minutos. Paquetes principales: terra, geodata, rnaturalearth, ggplot2, sf. Flujo metodológico adaptado del material de clase de la Maestría en Ciencia de Datos.