1 Introducción y planteamiento del caso

Un agricultor quiere saber en qué zonas del mundo puede establecer caña de azúcar y obtener rendimientos potenciales. El clima no es el único factor limitante (también influyen suelo, riego, manejo, variedad), pero sí es el primer filtro. Según la información del caso, los rangos climáticos óptimos son:

Variable Rango óptimo
Temperatura media 22,5 – 28 °C
Precipitación anual 1.500 – 3.500 mm
Precipitación mensual 125 – 290 mm

Además, se sabe que el Valle del Cauca (Colombia) tiene rendimientos altos. Esto permite dos aproximaciones complementarias, que se desarrollan en este informe:

  1. Aproximación por rangos óptimos (basada en conocimiento): se marcan como aptas las celdas del planeta cuyo clima cae dentro de los rangos de la tabla.
  2. Aproximación por similaridad (basada en un sitio de referencia): se buscan las celdas del planeta cuyo clima se parece al de sitios del Valle del Cauca, usando la distancia euclidiana sobre variables climáticas estandarizadas.

Al final se comparan ambos resultados.

2 Preparación: paquetes y datos

2.1 Paquetes

Se usa terra para el manejo de ráster y vectores (es el sucesor de raster y es mucho más rápido), y geodata para descargar los datos de clima de línea base (WorldClim) y los límites de países.

# Instala los paquetes que falten (solo la primera vez)
paquetes <- c("terra", "geodata", "knitr")
faltantes <- paquetes[!paquetes %in% rownames(installed.packages())]
if (length(faltantes) > 0) install.packages(faltantes)

library(terra)    # rásters y vectores espaciales
library(geodata)  # descarga de WorldClim y límites administrativos
library(knitr)    # tablas (kable)

2.2 Datos de clima de línea base (WorldClim 2.1)

Se usa WorldClim 2.1, que es la línea base climática global (promedios del periodo 1970–2000). Cada variable viene como un ráster de 12 capas, una por mes:

  • tavg: temperatura media mensual (°C).
  • prec: precipitación total mensual (mm).

Se trabaja a una resolución de 10 minutos de arco (~18 km en el ecuador), suficiente para un análisis global y liviana para procesar. La descarga se hace una sola vez: geodata guarda los archivos en la carpeta datos/ y en las siguientes ejecuciones los lee de allí.

dir.create("datos", showWarnings = FALSE)

# Temperatura media mensual (12 capas) y precipitación mensual (12 capas)
tavg <- worldclim_global(var = "tavg", res = 10, path = "datos")
prec <- worldclim_global(var = "prec", res = 10, path = "datos")

# Nombres de capas más legibles: ene, feb, ..., dic
meses <- c("Ene", "Feb", "Mar", "Abr", "May", "Jun",
           "Jul", "Ago", "Sep", "Oct", "Nov", "Dic")
names(tavg) <- meses
names(prec) <- meses

tavg
## 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       :     Ene,        Feb,       Mar,        Abr,        May,        Jun, ...
## 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, ...

2.3 Límites de países (shape global)

# Límites de países (nivel 0 = país). resolution = 5 es la versión simplificada,
# adecuada para mapas globales.
mundo <- world(resolution = 5, level = 0, path = "datos")
mundo
## class       : SpatVector
## geometry    : polygons
## dimensions  : 231, 2  (geometries, attributes)
## extent      : -180, 180, -90, 83.65625  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names       : GID_0      NAME_0
## type        : <chr>       <chr>
## values      :   ABW       Aruba
##                 AFG Afghanistan
##                 AGO      Angola
##               ...

Nota: si el curso entregó sus propios archivos (por ejemplo un shape de países en .shp), basta con reemplazar la línea anterior por mundo <- vect("ruta/al/archivo.shp"). El resto del código no cambia.

3 Paso 1 — Mapas de aptitud climática por rangos óptimos

3.1 Variables anuales

Los rangos del caso se expresan en términos anuales (temperatura media y precipitación acumulada), así que primero se resumen las 12 capas mensuales:

  • Temperatura media anual = promedio de las 12 temperaturas mensuales.
  • Precipitación anual = suma de las 12 precipitaciones mensuales.
tmed_anual <- mean(tavg)   # °C
prec_anual <- sum(prec)    # mm/año
names(tmed_anual) <- "tmed_anual"
names(prec_anual) <- "prec_anual"

# Resumen rápido de los valores globales
summary(c(tmed_anual, prec_anual))
##    tmed_anual         prec_anual    
##  Min.   :-54.7195   Min.   :   0.0  
##  1st Qu.:-22.8111   1st Qu.: 110.0  
##  Median : -0.5266   Median : 337.0  
##  Mean   : -4.0459   Mean   : 550.3  
##  3rd Qu.: 19.0177   3rd Qu.: 690.0  
##  Max.   : 30.7095   Max.   :7011.0  
##  NAs    :65591      NAs    :65591
par(mfrow = c(2, 1))

# Temperatura: paleta divergente (azul = frío, rojo = calor)
plot(tmed_anual,
     col = hcl.colors(100, "RdYlBu", rev = TRUE),
     main = "Temperatura media anual (°C) — WorldClim 1970-2000",
     axes = FALSE, mar = c(1, 1, 2, 5))
lines(mundo, lwd = 0.2, col = "grey30")

# Precipitación: paleta secuencial (amarillo = seco, azul = húmedo).
# Se recorta en 5.000 mm para que la escala no la dominen unos pocos píxeles extremos.
plot(clamp(prec_anual, upper = 5000),
     col = hcl.colors(100, "YlGnBu", rev = TRUE),
     main = "Precipitación anual (mm) — WorldClim 1970-2000 (valores > 5.000 mm saturados)",
     axes = FALSE, mar = c(1, 1, 2, 5))
lines(mundo, lwd = 0.2, col = "grey30")

3.2 Aptitud por cada variable

Se aplica una reclasificación booleana: una celda es apta (1) si su valor está dentro del rango óptimo, y no apta (0) si no.

# Rangos óptimos del caso
T_MIN <- 22.5; T_MAX <- 28
P_MIN <- 1500; P_MAX <- 3500

apt_temp <- (tmed_anual >= T_MIN) & (tmed_anual <= T_MAX)
apt_prec <- (prec_anual >= P_MIN) & (prec_anual <= P_MAX)

col_bin <- c("grey85", "#1a9850")  # no apto / apto

par(mfrow = c(2, 1))
plot(apt_temp, col = col_bin, type = "classes",
     levels = c("No apto", "Apto"),
     main = "Aptitud térmica: temperatura media entre 22,5 y 28 °C",
     axes = FALSE, mar = c(1, 1, 2, 7))
lines(mundo, lwd = 0.2, col = "grey40")

plot(apt_prec, col = col_bin, type = "classes",
     levels = c("No apto", "Apto"),
     main = "Aptitud hídrica: precipitación anual entre 1.500 y 3.500 mm",
     axes = FALSE, mar = c(1, 1, 2, 7))
lines(mundo, lwd = 0.2, col = "grey40")

3.3 Aptitud climática combinada

Para combinar ambas condiciones sin perder información se construye un mapa categórico con una suma ponderada (temp × 1 + prec × 2), que da 4 clases únicas:

Valor Clase
0 No apto
1 Solo temperatura apta (limitado por lluvia)
2 Solo precipitación apta (limitado por temperatura)
3 Apto (temperatura y precipitación)

Las clases 1 y 2 son útiles para el agricultor: la clase 1 indica zonas donde el clima térmico es ideal y la lluvia podría suplirse con riego.

# Versión numérica (0-3), útil para cálculos posteriores
aptitud_num <- apt_temp * 1 + apt_prec * 2
names(aptitud_num) <- "aptitud"

# Versión categórica (con etiquetas y colores fijos), para los mapas
aptitud <- aptitud_num * 1   # copia independiente
clases_apt <- data.frame(
  id    = 0:3,
  clase = c("No apto", "Solo temperatura", "Solo precipitación", "Apto (T y P)")
)
levels(aptitud) <- clases_apt

col_apt <- c("grey88", "#fdae61", "#74add1", "#006837")
# Se fija la tabla de colores al ráster: así cada clase conserva su color
# aunque, al recortar un país, alguna clase no aparezca.
coltab(aptitud) <- data.frame(value = 0:3, col = col_apt)

plot(aptitud,
     main = "Aptitud climática para caña de azúcar (rangos óptimos)",
     axes = FALSE, mar = c(1, 1, 2, 11))
lines(mundo, lwd = 0.2, col = "grey40")

Para cuantificar el resultado se calcula la superficie de cada clase. Como las celdas en coordenadas geográficas no tienen todas el mismo tamaño (se achican hacia los polos), se usa cellSize() para obtener el área real de cada celda en km².

area_celda <- cellSize(aptitud, unit = "km")
area_clase <- zonal(area_celda, aptitud, fun = "sum", na.rm = TRUE)
names(area_clase) <- c("Clase", "Area_km2")
area_clase$Porcentaje <- round(100 * area_clase$Area_km2 / sum(area_clase$Area_km2), 2)
area_clase$Area_km2 <- format(round(area_clase$Area_km2), big.mark = ".")

kable(area_clase, caption = "Superficie terrestre por clase de aptitud climática",
      col.names = c("Clase", "Área (km²)", "% de tierra firme"))
Superficie terrestre por clase de aptitud climática
Clase Área (km²) % de tierra firme
No apto 103.335.958 66.95
Solo temperatura 28.916.692 18.73
Solo precipitación 5.266.034 3.41
Apto (T y P) 16.830.487 10.90

3.4 Aptitud con criterio mensual (complemento)

El caso también da un rango mensual de precipitación (125–290 mm). El promedio anual puede ocultar estaciones secas marcadas: un sitio con 1.800 mm concentrados en 4 meses cumple el criterio anual pero sufre déficit hídrico el resto del año. Por eso se cuenta, para cada celda, cuántos meses del año cumplen simultáneamente los dos rangos (temperatura mensual entre 22,5 y 28 °C y precipitación mensual entre 125 y 290 mm).

# Cada comparación devuelve 12 capas TRUE/FALSE; la suma cuenta los meses que cumplen
mes_apto <- (tavg >= T_MIN & tavg <= T_MAX) & (prec >= 125 & prec <= 290)
n_meses_aptos <- sum(mes_apto)
names(n_meses_aptos) <- "meses_aptos"

plot(n_meses_aptos,
     col = c("grey88", hcl.colors(12, "YlGn", rev = TRUE)),
     breaks = -0.5:12.5, type = "interval",
     plg = list(legend = 0:12, title = "Meses"),
     main = "Meses al año con temperatura y precipitación mensual óptimas",
     axes = FALSE, mar = c(1, 1, 2, 6))
lines(mundo, lwd = 0.2, col = "grey40")

3.5 Interpretación del paso 1

  • Las zonas aptas en ambas variables se concentran en la franja tropical (aprox. entre 20° N y 20° S): norte de Sudamérica y Brasil, Centroamérica y el Caribe, África occidental y central, sur y sureste asiático (India, Bangladesh, Indochina, Indonesia, Filipinas) y el extremo norte de Australia.
  • La temperatura es el filtro más restrictivo en latitud: fuera del trópico prácticamente ninguna zona alcanza una media anual de 22,5 °C. La precipitación es el filtro dentro del trópico: excluye los desiertos cálidos (Sahara, Arabia, Australia central), que aparecen como “solo temperatura”, y las zonas extremadamente lluviosas (> 3.500 mm), como el Chocó biogeográfico colombiano o partes de la Amazonía occidental.
  • La clase “solo temperatura” agrupa dos situaciones distintas: zonas secas (Sahel, India, norte de Australia, noreste de Brasil), donde el déficit de lluvia podría compensarse con riego, y zonas demasiado lluviosas (> 3.500 mm, como la costa Pacífica colombiana), donde el exceso de agua, el drenaje y las enfermedades son el problema.
  • El mapa mensual es más exigente: muestra que muchas zonas aptas por el criterio anual tienen una estación seca marcada (pocos meses óptimos), mientras que las zonas con 9–12 meses óptimos tienen una lluvia bien distribuida a lo largo del año.

4 Paso 2 — Recorte para países con alto potencial

4.1 Selección objetiva de países: ranking por área apta

Visualmente, el mapa global muestra que las zonas aptas se concentran cerca de la línea ecuatorial. Para no elegir los países “a ojo”, se calcula para cada país:

  • Área apta (T y P) en km²: suma del área de las celdas de clase 3 dentro del polígono del país.
  • % del territorio apto: área apta / área total del país.
  • Área “solo temperatura”: zonas térmicamente ideales que requerirían riego.

Técnicamente se hace una extracción zonal (extract() con fun = sum) del ráster de áreas sobre los polígonos del shape global.

# Ráster de área (km²) de las celdas aptas; las demás valen 0
area_apta      <- ifel(aptitud_num == 3, area_celda, 0)
area_solo_temp <- ifel(aptitud_num == 1, area_celda, 0)
area_tierra    <- mask(area_celda, aptitud_num)   # solo celdas con dato (tierra)

# Suma por país (una fila por polígono del shape)
sumar <- function(r) extract(r, mundo, fun = sum, na.rm = TRUE, ID = FALSE)[, 1]

ranking <- data.frame(
  Pais           = mundo$NAME_0,
  Codigo         = mundo$GID_0,
  Area_apta      = sumar(area_apta),
  Area_solo_temp = sumar(area_solo_temp),
  Area_total     = sumar(area_tierra)
)
ranking$Pct_apto <- 100 * ranking$Area_apta / ranking$Area_total
ranking <- ranking[order(-ranking$Area_apta), ]

# Tabla de los 15 países con mayor área apta
top15 <- head(ranking, 15)
kable(
  data.frame(
    Pais = top15$Pais,
    `Área apta (km²)` = format(round(top15$Area_apta), big.mark = "."),
    `% del país apto` = round(top15$Pct_apto, 1),
    `Área solo temp. (km²)` = format(round(top15$Area_solo_temp), big.mark = "."),
    check.names = FALSE
  ),
  caption = "Países con mayor superficie climáticamente apta (temperatura y precipitación)"
)
Países con mayor superficie climáticamente apta (temperatura y precipitación)
Pais Área apta (km²) % del país apto Área solo temp. (km²)
Brazil 4.949.701 58.2 1.986.635
Democratic Republic of the Congo 1.470.451 63.1 479.628
Indonesia 1.375.329 73.8 242.578
Colombia 760.789 66.9 119.201
Peru 625.184 48.4 64.665
Venezuela 534.487 58.6 286.797
Bolivia 346.183 32.0 361.048
Cameroon 318.027 68.5 106.734
India 313.937 10.0 2.284.972
Congo 288.394 85.2 48.773
Gabon 257.386 96.7 5.803
Philippines 251.720 88.1 17.547
Malaysia 250.602 76.5 58.416
Papua New Guinea 246.838 53.6 85.626
Myanmar 213.069 32.0 212.293 
par(mar = c(4, 10, 2, 1))
barplot(rev(top15$Area_apta) / 1e3, names.arg = rev(top15$Pais),
        horiz = TRUE, las = 1, col = "#006837", border = NA,
        xlab = "Área apta (miles de km²)",
        main = "Top 15 de países por área climáticamente apta")

4.2 Países seleccionados

Con base en el ranking se seleccionan tres países, uno por continente, para representar las tres grandes regiones aptas del trópico:

  • Brasil (BRA): el país con mayor área apta del mundo (Amazonía) y, además, el primer productor mundial de caña de azúcar.
  • Colombia (COL): país de referencia del caso (Valle del Cauca) y con amplias zonas aptas en la Orinoquía y la Amazonía.
  • Indonesia (IDN): la mayor área apta del sureste asiático, con clima ecuatorial húmedo y lluvia bien distribuida en el año.
# Para analizar otros países basta con cambiar los códigos ISO3 (ej. "IND", "COD")
paises_sel <- c("BRA", "COL", "IDN")
ranking[ranking$Codigo %in% paises_sel, c("Pais", "Area_apta", "Pct_apto")]
##         Pais Area_apta Pct_apto
## 30    Brazil 4949700.8 58.22146
## 95 Indonesia 1375328.8 73.81590
## 45  Colombia  760789.1 66.89459

4.3 Recorte de los mapas con el shape global

Para cada país se hace un corte en dos pasos:

  1. crop(): recorta el ráster al rectángulo (extensión) del país.
  2. mask(): deja en NA las celdas que quedan fuera del polígono, para que solo se vea el territorio del país.

Se añaden los límites departamentales/provinciales (GADM nivel 1) como referencia geográfica para interpretar dónde quedan las zonas aptas.

# Recorta y enmascara un ráster con el polígono de un país
recortar_pais <- function(r, cod) {
  pol <- mundo[mundo$GID_0 == cod, ]
  mask(crop(r, pol), pol)
}

# Grafica, para un país, el mapa de aptitud combinada y el de meses óptimos
graficar_pais <- function(cod) {
  pol   <- mundo[mundo$GID_0 == cod, ]
  nom   <- pol$NAME_0
  div1  <- gadm(country = cod, level = 1, path = "datos", resolution = 2)

  apt_p   <- recortar_pais(aptitud, cod)
  meses_p <- recortar_pais(n_meses_aptos, cod)

  par(mfrow = c(1, 2))
  plot(apt_p, all_levels = TRUE, axes = FALSE, mar = c(1, 1, 2, 9),
       main = paste0(nom, ": aptitud climática"))
  lines(div1, lwd = 0.3, col = "grey35")
  lines(pol, lwd = 1)

  plot(meses_p,
       col = c("grey88", hcl.colors(12, "YlGn", rev = TRUE)),
       breaks = -0.5:12.5, type = "interval",
       plg = list(legend = 0:12, title = "Meses"),
       axes = FALSE, mar = c(1, 1, 2, 5),
       main = paste0(nom, ": meses óptimos al año"))
  lines(div1, lwd = 0.3, col = "grey35")
  lines(pol, lwd = 1)
}

4.3.1 Brasil

graficar_pais("BRA")

4.3.2 Colombia

graficar_pais("COL")

4.3.3 Indonesia

graficar_pais("IDN")

4.4 Interpretación del paso 2

  • El ranking confirma que el potencial se concentra en países ecuatoriales: Brasil (≈ 4,95 millones de km² aptos), RD del Congo (≈ 1,47 M), Indonesia (≈ 1,38 M) y Colombia (≈ 0,76 M) encabezan la lista. Los tres países seleccionados están entre los cuatro primeros.
  • India es un caso aparte: solo el 10 % de su territorio es apto en ambas variables, pero tiene la mayor área “solo temperatura” del mundo (≈ 2,28 M km²) y es el segundo productor mundial de caña. Ahí el riego reemplaza la lluvia que falta.
  • Brasil concentra la mayor superficie apta del planeta (58 % del país), casi toda en la cuenca amazónica. Es un resultado importante para el agricultor: la zona “climáticamente óptima” coincide con selva tropical, donde existen restricciones ambientales y legales, suelos pobres y poca infraestructura. En cambio, la principal región cañera real de Brasil (São Paulo, en el centro-sur) queda en el límite o por fuera de los rangos óptimos. Esto confirma que el clima es una condición necesaria, pero no suficiente.
  • Colombia (67 % del país apto) muestra zonas aptas en la Orinoquía, la Amazonía y el Magdalena Medio. La costa Pacífica aparece como “solo temperatura” por exceso de lluvia, y la región andina como no apta o “solo precipitación” por temperaturas bajas debidas a la altitud. Lo más llamativo es el valle geográfico del río Cauca: aunque es la región de mayores rendimientos del país, no aparece como apta. A esta resolución (~18 km), cada celda del valle (~25 km de ancho) se mezcla con las laderas de las cordilleras Occidental y Central, y la temperatura promedio de la celda cae por debajo de 22,5 °C. Además, la parte plana del valle tiene una precipitación cercana al límite inferior del rango, y allí la caña se cultiva con riego. Esto justifica la segunda aproximación (pasos 3 y 4), que parte de sitios donde se sabe que la caña rinde bien.
  • Indonesia (74 % del país apto) es apta en la mayor parte de su territorio de baja altitud (Sumatra, Borneo/Kalimantan, Sulawesi, Papúa), con muchos meses óptimos gracias a su lluvia bien distribuida. Las zonas montañosas quedan excluidas por temperatura.

5 Paso 3 — Clima de puntos en el Valle del Cauca

5.1 Selección de los puntos al azar

La caña en el Valle del Cauca se cultiva en el valle geográfico del río Cauca (la parte plana, ~950–1.050 m s. n. m., entre Jamundí y La Victoria), no en las cordilleras ni en la costa Pacífica del departamento. Si los puntos se tomaran al azar en todo el departamento, podrían caer en zonas de montaña sin caña. Por eso:

  1. Se delimitó la parte plana del valle con un polígono cuyos vértices se tomaron de Google Maps (bordes de las cordilleras Occidental y Central).
  2. Se dividió el valle en tres tramos (sur, centro y norte) y se tomó un punto al azar en cada tramo (muestreo aleatorio estratificado). Así los puntos quedan repartidos a lo largo del valle y no amontonados.
  3. Se fijó la semilla (set.seed(42)) para que el resultado sea reproducible.
# Polígono aproximado de la parte plana del valle del río Cauca
# (vértices lon/lat tomados de Google Maps, de Jamundí a La Victoria)
coords_valle <- rbind(
  c(-76.52, 3.25), c(-76.45, 3.55), c(-76.38, 3.85), c(-76.28, 4.15), c(-76.14, 4.45),
  c(-76.02, 4.45), c(-76.15, 4.15), c(-76.28, 3.85), c(-76.26, 3.55), c(-76.32, 3.25),
  c(-76.52, 3.25)
)
valle_plano <- vect(coords_valle, type = "polygons", crs = "EPSG:4326")

# Tramos del valle (franjas de latitud)
tramos <- data.frame(
  tramo   = c("Sur", "Centro", "Norte"),
  lat_min = c(3.25, 3.65, 4.05),
  lat_max = c(3.65, 4.05, 4.45)
)

# Genera un punto aleatorio uniforme dentro del polígono y la franja de latitud
# (se sortean coordenadas hasta que una caiga dentro del polígono)
punto_al_azar <- function(pol, lat_min, lat_max) {
  e <- ext(pol)
  repeat {
    p <- vect(cbind(runif(1, e$xmin, e$xmax), runif(1, lat_min, lat_max)),
              crs = "EPSG:4326")
    if (is.related(p, pol, "intersects")) return(p)
  }
}

set.seed(42)
puntos <- do.call(rbind, lapply(seq_len(nrow(tramos)), function(i)
  punto_al_azar(valle_plano, tramos$lat_min[i], tramos$lat_max[i])))

xy <- crds(puntos)
puntos$id         <- paste0("P", 1:3)
puntos$tramo      <- tramos$tramo
puntos$referencia <- c("Rozo – Palmira", "Buga", "Roldanillo")  # municipio más cercano
puntos$lon        <- round(xy[, 1], 4)
puntos$lat        <- round(xy[, 2], 4)
puntos$google_maps <- sprintf("https://www.google.com/maps?q=%.4f,%.4f",
                              puntos$lat, puntos$lon)

kable(as.data.frame(puntos)[, c("id", "tramo", "referencia", "lat", "lon", "google_maps")],
      caption = "Puntos seleccionados al azar en la zona cañera del Valle del Cauca",
      col.names = c("Punto", "Tramo", "Referencia", "Latitud", "Longitud", "Ver en Google Maps"))
Puntos seleccionados al azar en la zona cañera del Valle del Cauca
Punto Tramo Referencia Latitud Longitud Ver en Google Maps
P1 Sur Rozo – Palmira 3.5822 -76.3769 https://www.google.com/maps?q=3.5822,-76.3769
P2 Centro Buga 3.9376 -76.2911 https://www.google.com/maps?q=3.9376,-76.2911
P3 Norte Roldanillo 4.3744 -76.1512 https://www.google.com/maps?q=4.3744,-76.1512 
# Límite del departamento del Valle del Cauca (GADM nivel 1)
col_div1   <- gadm(country = "COL", level = 1, path = "datos", resolution = 2)
valle_dpto <- col_div1[col_div1$NAME_1 == "Valle del Cauca", ]

# Aptitud por rangos (paso 1) recortada al departamento
apt_valle <- mask(crop(aptitud, valle_dpto), valle_dpto)

col_pts <- c("#1b9e77", "#d95f02", "#7570b3")  # un color por punto

plot(apt_valle, all_levels = TRUE, axes = TRUE, mar = c(2, 2, 2, 10),
     main = "Valle del Cauca: aptitud por rangos (10') y puntos de referencia")
lines(valle_dpto, lwd = 1.2)
lines(valle_plano, lwd = 1.5, lty = 2, col = "black")
points(puntos, pch = 21, bg = col_pts, cex = 2)
text(puntos, labels = puntos$id, pos = 4, font = 2, cex = 0.9)

5.2 Extracción de la información de clima

Para cada punto se extraen los 12 valores mensuales de temperatura media y precipitación con extract(). Se usan dos resoluciones:

  • 10’ (~18 km): la misma de los mapas globales de los pasos 1 y 2.
  • 30” (~1 km): WorldClim a resolución fina para Colombia, que representa mejor el clima real de la parte plana del valle, porque no se mezcla con las laderas de las cordilleras.
# Clima a 30" (~1 km) solo para Colombia
tavg_col <- worldclim_country("COL", var = "tavg", path = "datos")
prec_col <- worldclim_country("COL", var = "prec", path = "datos")
names(tavg_col) <- meses
names(prec_col) <- meses

# Extracción: una fila por punto, una columna por mes
t_10  <- extract(tavg,     puntos, ID = FALSE)
p_10  <- extract(prec,     puntos, ID = FALSE)
t_30s <- extract(tavg_col, puntos, ID = FALSE)
p_30s <- extract(prec_col, puntos, ID = FALSE)

# Resumen anual a ambas resoluciones
resumen_pts <- data.frame(
  Punto        = puntos$id,
  Referencia   = puntos$referencia,
  Tmed_10min   = round(rowMeans(t_10), 1),
  Tmed_30s     = round(rowMeans(t_30s), 1),
  Prec_10min   = round(rowSums(p_10)),
  Prec_30s     = round(rowSums(p_30s)),
  Meses_opt_30s = rowSums(t_30s >= T_MIN & t_30s <= T_MAX &
                          p_30s >= 125 & p_30s <= 290)
)
kable(resumen_pts,
      caption = "Clima anual de los puntos a 10' y 30\" de resolución",
      col.names = c("Punto", "Referencia", "T media 10' (°C)", "T media 30\" (°C)",
                    "Prec. anual 10' (mm)", "Prec. anual 30\" (mm)",
                    "Meses óptimos (30\")"))
Clima anual de los puntos a 10’ y 30” de resolución
Punto Referencia T media 10’ (°C) T media 30” (°C) Prec. anual 10’ (mm) Prec. anual 30” (mm) Meses óptimos (30”)
P1 Rozo – Palmira 23.5 23.6 1485 1446 5
P2 Buga 21.7 23.6 1600 1507 6
P3 Roldanillo 23.4 23.6 1632 1535 6

5.3 Series de tiempo de temperatura y precipitación

Los datos de línea base de WorldClim son una climatología mensual: cada mes es el promedio de 1970–2000. Por eso la “serie de tiempo” de cada punto es el ciclo anual (enero a diciembre). Se grafican los datos a 30”, que representan el clima real del cultivo. La franja verde marca el rango óptimo de la caña.

par(mfrow = c(2, 1), mar = c(4, 4.5, 3, 1))

# --- Temperatura media mensual ---
matplot(1:12, t(as.matrix(t_30s)), type = "n", xaxt = "n",
        ylim = range(c(18, 30, unlist(t_30s))),
        xlab = "", ylab = "Temperatura media (°C)",
        main = "Ciclo anual de temperatura media — puntos del Valle del Cauca (30\")")
rect(0.5, T_MIN, 12.5, T_MAX, col = adjustcolor("#1a9850", 0.15), border = NA)
matlines(1:12, t(as.matrix(t_30s)), lty = 1, lwd = 2.5, col = col_pts)
matpoints(1:12, t(as.matrix(t_30s)), pch = 19, col = col_pts)
axis(1, at = 1:12, labels = meses)
legend("bottomright", legend = paste(puntos$id, puntos$referencia),
       col = col_pts, lwd = 2.5, bty = "n", cex = 0.85)
legend("topleft", legend = "Rango óptimo (22,5–28 °C)",
       fill = adjustcolor("#1a9850", 0.15), border = NA, bty = "n", cex = 0.85)

# --- Precipitación mensual ---
matplot(1:12, t(as.matrix(p_30s)), type = "n", xaxt = "n",
        ylim = c(0, max(300, unlist(p_30s))),
        xlab = "Mes", ylab = "Precipitación (mm)",
        main = "Ciclo anual de precipitación — puntos del Valle del Cauca (30\")")
rect(0.5, 125, 12.5, 290, col = adjustcolor("#1a9850", 0.15), border = NA)
matlines(1:12, t(as.matrix(p_30s)), lty = 1, lwd = 2.5, col = col_pts)
matpoints(1:12, t(as.matrix(p_30s)), pch = 19, col = col_pts)
axis(1, at = 1:12, labels = meses)
legend("topright", legend = "Rango óptimo (125–290 mm/mes)",
       fill = adjustcolor("#1a9850", 0.15), border = NA, bty = "n", cex = 0.85)

5.4 Interpretación del paso 3

  • Temperatura: a 30” los tres puntos tienen una temperatura media de 23,6 °C, dentro del rango óptimo (22,5–28 °C) y casi constante durante el año, como corresponde a un clima ecuatorial. Es una temperatura ideal, pero cercana al límite inferior del rango.
  • Efecto de la resolución: en P1 y P3 los valores a 10’ y a 30” son casi iguales, pero en P2 (Buga) la temperatura baja de 23,6 °C a 21,7 °C a 10’. Su celda de ~18 km incluye laderas de la cordillera Central, y el promedio cae por debajo de 22,5 °C. Por eso parte del valle aparece como “no apto” en los mapas globales, aunque en el terreno sí sea apto. Es un efecto de la escala, no del clima real del cultivo.
  • Precipitación: la lluvia anual a 30” está en el límite inferior del rango óptimo: 1.446 mm en P1 (por debajo de 1.500) y 1.507–1.535 mm en P2 y P3 (apenas dentro). El régimen es bimodal, con dos temporadas de lluvia y dos secas, y solo 5–6 meses al año cumplen a la vez los rangos mensuales de temperatura y precipitación.
  • Conclusión del paso 3: la región con los mayores rendimientos de caña del país está al borde de los rangos “óptimos”. Con un criterio binario basta una pequeña diferencia (P1: 54 mm menos de lluvia; P2 a 10’: la mezcla con la montaña) para que una zona excelente quede excluida. Los altos rendimientos del Valle se explican por la combinación de temperatura ideal y estable, alta radiación, suelos fértiles y riego en las épocas secas. Esto muestra que la aproximación por rangos puede excluir zonas productivas y motiva usar el clima real de estos sitios como referencia (paso 4).

6 Paso 4 — Mapas de similaridad climática (distancia de Mahalanobis)

6.1 ¿Por qué Mahalanobis y no la distancia euclidiana?

La distancia euclidiana mide la separación “en línea recta” entre el clima de una celda y el clima de referencia. Tiene dos problemas con variables climáticas:

  1. Escalas distintas: una diferencia de 1 °C y una de 1 mm pesarían lo mismo. Se puede corregir estandarizando, pero queda el segundo problema.
  2. Correlación entre variables: la lluvia anual y la lluvia del mes más seco, por ejemplo, están muy correlacionadas. Con la euclidiana esa información se cuenta dos veces.

La distancia de Mahalanobis corrige ambos problemas usando la matriz de covarianza (\(\Sigma\)) del clima global:

\[ D^2(\mathbf{x}) = (\mathbf{x} - \boldsymbol{\mu}_{ref})^\top \, \Sigma^{-1} \, (\mathbf{x} - \boldsymbol{\mu}_{ref}) \]

donde \(\mathbf{x}\) es el clima de una celda y \(\boldsymbol{\mu}_{ref}\) el clima del punto de referencia. Si las variables no estuvieran correlacionadas y tuvieran varianza 1, \(D\) sería igual a la distancia euclidiana.

6.2 Variables climáticas para comparar

En vez de comparar los 24 valores mensuales mes a mes, se usan 5 variables resumen. Esto evita un problema importante: en el hemisferio sur las estaciones están invertidas, así que un sitio con el mismo clima que el Valle pero con las lluvias desfasadas seis meses parecería “distinto” si se compara enero con enero.

Variable Qué representa para la caña
Temperatura media anual (°C) Régimen térmico general
Rango térmico anual (°C) = T máx. mensual − T mín. mensual Estabilidad térmica (el Valle casi no tiene estaciones)
Precipitación anual (mm) Disponibilidad total de agua
Precipitación del mes más seco (mm) Severidad de la época seca (necesidad de riego)
Coeficiente de variación de la precipitación mensual (%) Qué tan concentrada está la lluvia en el año 
# --- Variables resumen para la grilla global (10') ---
rango_t  <- max(tavg) - min(tavg)
prec_min <- min(prec)
prec_med <- mean(prec)
prec_cv  <- stdev(prec, pop = TRUE) / ifel(prec_med < 1, 1, prec_med) * 100

vars_glob <- c(tmed_anual, rango_t, prec_anual, prec_min, prec_cv)
names(vars_glob) <- c("tmed", "rango_t", "prec_anual", "prec_min", "prec_cv")

# --- Las mismas variables para los puntos del Valle (30", clima real del cultivo) ---
resumir_puntos <- function(tm, pm) {
  tm <- as.matrix(tm); pm <- as.matrix(pm)
  sd_pob <- function(x) sqrt(mean((x - mean(x))^2))   # desviación poblacional
  data.frame(
    tmed       = rowMeans(tm),
    rango_t    = apply(tm, 1, max) - apply(tm, 1, min),
    prec_anual = rowSums(pm),
    prec_min   = apply(pm, 1, min),
    prec_cv    = apply(pm, 1, sd_pob) / pmax(rowMeans(pm), 1) * 100
  )
}
ref_pts <- resumir_puntos(t_30s, p_30s)
rownames(ref_pts) <- paste(puntos$id, puntos$referencia)

kable(round(ref_pts, 1),
      caption = "Clima de referencia de los puntos del Valle del Cauca (30\")",
      col.names = c("T media (°C)", "Rango térmico (°C)", "Prec. anual (mm)",
                    "Prec. mes más seco (mm)", "CV prec. (%)"))
Clima de referencia de los puntos del Valle del Cauca (30”)
T media (°C) Rango térmico (°C) Prec. anual (mm) Prec. mes más seco (mm) CV prec. (%)
P1 Rozo – Palmira 23.6 0.8 1446 50 37.0
P2 Buga 23.6 1.0 1507 59 36.7
P3 Roldanillo 23.6 1.2 1535 73 35.9

6.3 Matriz de covarianza: ¿sobre qué dominio?

La matriz de covarianza \(\Sigma\) define la “escala” con la que se mide cada diferencia. Por eso es clave decidir sobre qué celdas se estima:

  • Dominio global (toda la tierra firme): incluye climas polares y desérticos, así que la variabilidad de la temperatura es enorme (de −50 °C a +30 °C). Con esa escala, una diferencia de 14 °C (por ejemplo, entre el Valle y Europa occidental) equivale a menos de una desviación estándar y casi no se penaliza. Además, a escala global “más frío” va casi siempre con “más estacional”, y Mahalanobis da poco peso a las diferencias que siguen esa correlación. El resultado (que se muestra más adelante como comparación) asigna alta similaridad a Europa, Australia e incluso a zonas de la Antártida, lo cual no tiene sentido agronómico.
  • Dominio de crecimiento potencial de la caña: celdas cuyo mes más frío tiene una temperatura media ≥ 10 °C (climas tropicales y subtropicales sin heladas, donde la caña puede crecer físicamente). Dentro de este dominio la variabilidad es la relevante para el cultivo, así que una diferencia de unos pocos grados ya se considera importante.

Se usa el dominio de crecimiento potencial. En ambos casos la covarianza se estima con una muestra regular de celdas de tierra firme, que es suficiente para estimarla con precisión sin cargar el ráster completo en memoria.

# Dominio de crecimiento potencial: mes más frío con T media >= 10 °C
T_FRIO_MIN <- 10
dominio    <- min(tavg) >= T_FRIO_MIN
vars_dom   <- mask(vars_glob, dominio, maskvalues = 0)

# Muestra regular de celdas; na.omit() descarta explícitamente océanos y celdas
# fuera del dominio (con method = "regular", na.rm no siempre las elimina)
muestra_dom  <- na.omit(spatSample(vars_dom,  size = 3e5, method = "regular",
                                   as.df = TRUE))
muestra_glob <- na.omit(spatSample(vars_glob, size = 1e5, method = "regular",
                                   as.df = TRUE))
c(celdas_dominio = nrow(muestra_dom), celdas_global = nrow(muestra_glob))
## celdas_dominio  celdas_global 
##          27521          34761
S      <- cov(muestra_dom)    # covarianza usada en el análisis
S_glob <- cov(muestra_glob)   # solo para comparación

# Desviación estándar de cada variable en ambos dominios:
# muestra cuánto se "diluyen" las diferencias con el dominio global
kable(round(rbind(`Dominio caña (≥ 10 °C)` = sqrt(diag(S)),
                  `Global` = sqrt(diag(S_glob))), 1),
      caption = "Desviación estándar de las variables según el dominio")
Desviación estándar de las variables según el dominio
tmed rango_t prec_anual prec_min prec_cv
Dominio caña (≥ 10 °C) 3.2 6.4 927.6 44.6 39.0
Global 24.5 12.5 652.1 28.7 62.8 
# Correlaciones entre variables en el dominio usado
kable(round(cor(muestra_dom), 2),
      caption = "Correlación entre las variables climáticas (dominio de crecimiento potencial)")
Correlación entre las variables climáticas (dominio de crecimiento potencial)
tmed rango_t prec_anual prec_min prec_cv
tmed 1.00 -0.12 0.11 0.03 0.33
rango_t -0.12 1.00 -0.71 -0.43 0.30
prec_anual 0.11 -0.71 1.00 0.73 -0.47
prec_min 0.03 -0.43 0.73 1.00 -0.59
prec_cv 0.33 0.30 -0.47 -0.59 1.00

6.4 Cálculo de la distancia y de la similaridad

Para calcular \(D^2\) en todo el planeta de forma eficiente se usa la descomposición de Cholesky \(\Sigma^{-1} = U^\top U\). Así, \(D^2 = \lVert U(\mathbf{x}-\boldsymbol{\mu}_{ref}) \rVert^2\), que se calcula con álgebra de rásters (sumas y productos de capas) sin recorrer celda por celda.

Para convertir la distancia en una similaridad entre 0 y 1 se usa la distribución chi-cuadrado: si el clima fuera normal multivariado, \(D^2\) seguiría una \(\chi^2\) con 5 grados de libertad (una por variable). Entonces:

\[ \text{Similaridad} = 1 - F_{\chi^2_5}(D^2) \]

Una similaridad de 1 indica un clima idéntico al de referencia, y valores cercanos a 0 indican climas muy diferentes.

U      <- chol(solve(S))        # matriz triangular superior: S^-1 = t(U) %*% U
U_glob <- chol(solve(S_glob))   # versión con covarianza global (comparación)

# Distancia de Mahalanobis al cuadrado entre cada celda y un vector de referencia
mahalanobis_raster <- function(r, ref, U) {
  d <- r - as.numeric(ref)                  # diferencia capa por capa
  z <- lapply(seq_len(nrow(U)), function(k) sum(d * U[k, ]))  # z = U %*% d
  sum(rast(z)^2)                            # D² = suma de z²
}

# Similaridad (0-1) a partir de D²
similaridad_raster <- function(r, ref, U) {
  d2 <- mahalanobis_raster(r, ref, U)
  # lapp() pasa la capa completa como vector: pchisq() se aplica vectorizado (rápido)
  lapp(d2, function(x) 1 - pchisq(x, df = 5))
}

# Comprobación: en una celda cualquiera el resultado coincide con stats::mahalanobis()
celda <- as.numeric(muestra_dom[1, ])
c(raster = sum((U %*% (celda - as.numeric(ref_pts[1, ])))^2),
  stats  = mahalanobis(celda, as.numeric(ref_pts[1, ]), S))
##   raster    stats 
## 11.05267 11.05267
# Un mapa de similaridad por cada punto de referencia
sim_pts <- rast(lapply(seq_len(nrow(ref_pts)), function(i)
  similaridad_raster(vars_glob, ref_pts[i, ], U)))
names(sim_pts) <- puntos$id

# Similaridad combinada: la celda se parece a al menos uno de los tres puntos
sim_max <- max(sim_pts)
names(sim_max) <- "similaridad"

6.5 Mapas de similaridad

Se usa una escala por intervalos para que las clases sean fáciles de leer. Los tonos oscuros indican los climas más parecidos al del Valle del Cauca.

cortes_sim <- c(0, 0.05, 0.10, 0.25, 0.50, 0.75, 1)
col_sim    <- c("grey88", hcl.colors(5, "YlGnBu", rev = TRUE))
etq_sim    <- c("< 0,05", "0,05–0,10", "0,10–0,25", "0,25–0,50", "0,50–0,75", "> 0,75")

graficar_sim <- function(r, titulo) {
  plot(r, breaks = cortes_sim, col = col_sim, type = "interval",
       plg = list(legend = etq_sim, title = "Similaridad"),
       axes = FALSE, mar = c(1, 1, 2, 8), main = titulo)
  lines(mundo, lwd = 0.2, col = "grey40")
}

par(mfrow = c(3, 1))
for (i in 1:3) {
  graficar_sim(sim_pts[[i]],
               paste0("Similaridad climática con ", puntos$id[i], " (",
                      puntos$referencia[i], ")"))
  points(puntos[i], pch = 21, bg = "red", cex = 1.2)
}

graficar_sim(sim_max, "Similaridad climática con el Valle del Cauca (máximo de P1, P2 y P3)")
points(puntos, pch = 21, bg = "red", cex = 1)

6.5.1 Comparación: covarianza global vs. dominio de crecimiento potencial

Para mostrar la importancia de la elección del dominio, se repite el mapa combinado usando la covarianza global.

sim_max_glob <- max(rast(lapply(seq_len(nrow(ref_pts)), function(i)
  similaridad_raster(vars_glob, ref_pts[i, ], U_glob))))

par(mfrow = c(2, 1))
graficar_sim(sim_max_glob, "Similaridad con covarianza GLOBAL (no recomendada)")
graficar_sim(sim_max, "Similaridad con covarianza del dominio de la caña (usada)")

Con la covarianza global, Europa occidental, el sur de Australia y partes de la Antártida aparecen como “similares” al Valle del Cauca, aunque su temperatura media es mucho más baja (de 10 °C a más de 40 °C menos). Con la covarianza del dominio de crecimiento potencial, esas zonas quedan descartadas y la similaridad alta se concentra en climas tropicales. El resultado de Mahalanobis depende de sobre qué población se mide la variabilidad. Ese dominio debe definirse con criterio agronómico.

6.6 Países con mayor área climáticamente similar

Igual que en el paso 2, para cada país se suma el área de las celdas con similaridad ≥ 0,5. Este umbral equivale a \(D^2 \le 4{,}35\), la mediana de la \(\chi^2_5\), y selecciona climas muy cercanos al de referencia.

UMBRAL_SIM <- 0.5
area_similar <- ifel(sim_max >= UMBRAL_SIM, area_celda, 0)

ranking_sim <- data.frame(
  Pais         = mundo$NAME_0,
  Codigo       = mundo$GID_0,
  Area_similar = sumar(area_similar),
  Area_total   = sumar(area_tierra)
)
ranking_sim$Pct_similar <- 100 * ranking_sim$Area_similar / ranking_sim$Area_total
ranking_sim <- ranking_sim[order(-ranking_sim$Area_similar), ]

top15_sim <- head(ranking_sim, 15)
kable(
  data.frame(
    Pais = top15_sim$Pais,
    `Área similar (km²)` = format(round(top15_sim$Area_similar), big.mark = "."),
    `% del país` = round(top15_sim$Pct_similar, 1),
    check.names = FALSE
  ),
  caption = "Países con mayor superficie de clima similar al Valle del Cauca (similaridad ≥ 0,5)"
)
Países con mayor superficie de clima similar al Valle del Cauca (similaridad ≥ 0,5)
Pais Área similar (km²) % del país
Brazil 5.058.133 59.5
Democratic Republic of the Congo 2.119.571 90.9
Indonesia 783.602 42.1
Bolivia 703.811 65.1
Tanzania 630.101 66.9
South Sudan 609.458 97.2
Central African Republic 598.484 96.7
Venezuela 528.892 58.0
Ethiopia 486.386 43.1
Peru 480.281 37.1
Mozambique 468.791 59.7
Angola 394.554 31.7
Paraguay 352.055 88.0
Mexico 325.262 16.6
Kenya 318.376 54.3

6.7 Interpretación del paso 4

  • Al usar la covarianza del dominio de la caña, la similaridad alta queda solo en la franja tropical, y las zonas templadas quedan descartadas.
  • Las zonas más parecidas al Valle están en África central y oriental (cuenca del Congo, Grandes Lagos), el norte de Sudamérica y Centroamérica y el Caribe: climas de 23–25 °C, estables en el año y con lluvias moderadas en dos temporadas.
  • La similaridad es menor en climas mucho más lluviosos (Amazonía occidental) o con una estación seca fuerte (Sahel, India, norte de Australia).
  • Los tres puntos dan mapas casi iguales porque tienen prácticamente el mismo clima: la zona cañera del valle es homogénea.

7 Paso 5 — Comparación de aproximaciones y conclusiones

7.1 Mapa de consenso entre ambas aproximaciones

Para comparar formalmente los resultados se cruzan los dos mapas binarios:

  • Apto por rangos: clase “Apto (T y P)” del paso 1.
  • Similar al Valle: similaridad combinada ≥ 0,5 del paso 4.

El cruce produce 4 clases, que se codifican con la misma técnica del paso 1 (rangos × 1 + similaridad × 2):

Valor Clase Significado para el agricultor
0 Ninguna Clima no recomendado
1 Solo rangos Clima “de libro” óptimo, pero distinto al de una zona de éxito comprobado
2 Solo similaridad Clima análogo al Valle del Cauca, aunque no cumpla los rangos estrictos
3 Ambas (consenso) Zonas prioritarias: cumplen los rangos y además se parecen al Valle 
apto_rangos <- aptitud_num == 3
similar     <- sim_max >= UMBRAL_SIM

consenso_num <- apto_rangos * 1 + similar * 2
consenso     <- consenso_num * 1
levels(consenso) <- data.frame(
  id    = 0:3,
  clase = c("Ninguna", "Solo rangos", "Solo similaridad", "Ambas (consenso)")
)
col_cons <- c("grey88", "#f1a340", "#998ec3", "#1b7837")
coltab(consenso) <- data.frame(value = 0:3, col = col_cons)

plot(consenso, all_levels = TRUE, axes = FALSE, mar = c(1, 1, 2, 10),
     main = "Comparación: aptitud por rangos vs. similaridad con el Valle del Cauca")
lines(mundo, lwd = 0.2, col = "grey40")
points(puntos, pch = 21, bg = "red", cex = 0.9)

# Zoom en Colombia, el país del agricultor
graficar_consenso_pais <- function(cod) {
  pol <- mundo[mundo$GID_0 == cod, ]
  plot(mask(crop(consenso, pol), pol), all_levels = TRUE, axes = FALSE,
       mar = c(1, 1, 2, 10), main = paste0(pol$NAME_0, ": comparación de aproximaciones"))
  lines(gadm(country = cod, level = 1, path = "datos", resolution = 2),
        lwd = 0.3, col = "grey35")
  lines(pol, lwd = 1)
}
graficar_consenso_pais("COL")
points(puntos, pch = 21, bg = "red", cex = 1.2)

7.2 Cuantificación del acuerdo

Se calcula el área de cada clase y el índice de Jaccard, que mide el grado de acuerdo entre las dos aproximaciones:

\[ J = \frac{\text{Área(ambas)}}{\text{Área(ambas)} + \text{Área(solo rangos)} + \text{Área(solo similaridad)}} \]

\(J = 1\) indica que ambos métodos señalan exactamente las mismas zonas, y \(J = 0\) que no coinciden en nada.

area_cons <- zonal(area_celda, consenso, fun = "sum", na.rm = TRUE)
names(area_cons) <- c("Clase", "Area_km2")
a <- setNames(area_cons$Area_km2, area_cons$Clase)

jaccard <- a["Ambas (consenso)"] /
  (a["Ambas (consenso)"] + a["Solo rangos"] + a["Solo similaridad"])

kable(
  data.frame(Clase = area_cons$Clase,
             `Área (km²)` = format(round(area_cons$Area_km2), big.mark = "."),
             check.names = FALSE),
  caption = "Superficie por clase de comparación"
)
Superficie por clase de comparación
Clase Área (km²)
Ninguna 124.796.424
Solo rangos 8.616.828
Solo similaridad 12.722.259
Ambas (consenso) 8.213.660 
cat(sprintf("Índice de Jaccard entre ambas aproximaciones: %.2f\n", jaccard))
## Índice de Jaccard entre ambas aproximaciones: 0.28

7.3 Países prioritarios según cada aproximación

area_ambas <- ifel(consenso_num == 3, area_celda, 0)
ranking_cons <- data.frame(Pais = mundo$NAME_0, Area = sumar(area_ambas))
ranking_cons <- ranking_cons[order(-ranking_cons$Area), ]

n_top <- min(10, nrow(ranking_cons))
kable(
  data.frame(
    Posicion            = seq_len(n_top),
    `Por rangos`        = head(ranking$Pais, n_top),
    `Por similaridad`   = head(ranking_sim$Pais, n_top),
    `Consenso (ambas)`  = head(ranking_cons$Pais, n_top),
    `Área consenso (km²)` = format(round(head(ranking_cons$Area, n_top)), big.mark = "."),
    check.names = FALSE
  ),
  caption = "Top 10 de países según cada aproximación"
)
Top 10 de países según cada aproximación
Posicion Por rangos Por similaridad Consenso (ambas) Área consenso (km²)
1 Brazil Brazil Brazil 2.533.984
2 Democratic Republic of the Congo Democratic Republic of the Congo Democratic Republic of the Congo 1.470.110
3 Indonesia Indonesia Indonesia 671.605
4 Colombia Bolivia Peru 385.540
5 Peru Tanzania Bolivia 301.088
6 Venezuela South Sudan Congo 239.569
7 Bolivia Central African Republic Guyana 204.696
8 Cameroon Venezuela Cameroon 204.670
9 India Ethiopia Venezuela 188.362
10 Congo Peru Central African Republic 161.828

7.4 Comparación e interpretación

  • Coinciden en lo esencial: ambas aproximaciones ubican el potencial en la franja tropical, y los tres primeros países son los mismos (Brasil, RD del Congo e Indonesia).
  • Difieren en el tipo de clima que favorecen: los rangos aceptan climas muy lluviosos (Amazonía, Colombia, Camerún), mientras que la similaridad favorece climas moderados y con dos temporadas de lluvia, como el del Valle. Por eso suben en el ranking países de África central y oriental (Tanzania, Sudán del Sur, República Centroafricana, Etiopía) y Bolivia.
  • Colombia sale del top 10 por similaridad y por consenso: su área apta está sobre todo en la Amazonía y la Orinoquía, más lluviosas que el Valle.
  • El propio Valle del Cauca queda al borde de los rangos óptimos, aunque tiene rendimientos altos. Los rangos fijos son una simplificación, y el riego y el manejo compensan las diferencias moderadas.

8 Conclusiones

  1. El clima apto para la caña está en el trópico. Los países con mayor potencial son Brasil, RD del Congo e Indonesia, por cualquiera de las dos aproximaciones.
  2. Las dos aproximaciones se complementan. Los rangos muestran dónde el clima es teóricamente óptimo, y la similaridad dónde se parece a una zona que ya produce bien. Las zonas de consenso son la mejor recomendación para el agricultor.
  3. Los rangos no deben usarse como límites rígidos: el Valle del Cauca no los cumple del todo y es muy productivo gracias al riego.
  4. Los métodos requieren criterio: la resolución de los datos y el dominio de la covarianza cambian bastante los resultados.
  5. El clima no es suficiente: antes de invertir, hay que evaluar suelos, agua para riego, restricciones ambientales (gran parte del área apta es selva) e infraestructura.