Aptitud climática de la caña de azúcar a escala global

Módulo 1 - Unidad 2 · Análisis de datos geográficos y espaciales

1 Introducción y objetivo

La caña de azúcar (Saccharum officinarum) es el cultivo más importante del valle geográfico del río Cauca. Para saber en qué otras partes del mundo se podría cultivar con un clima parecido, se usan datos de clima de línea base a nivel global: temperatura media mensual y precipitación mensual en una grilla regular (12 capas por variable, una por mes).

Fuente de datos. Se usa WorldClim 2.1 (Fick y Hijmans, 2017). Representa el clima promedio del periodo 1970-2000 con una resolución de 10 minutos de arco (unos 18,5 km en el ecuador), que es más fina que la de 0,5° del ejemplo.

Objetivo: identificar las zonas del mundo con condiciones climáticas aptas para la caña de azúcar con dos aproximaciones diferentes y compararlas:

1. Aptitud por rangos óptimos (enfoque ecofisiológico): a cada celda se le asigna un índice entre 0 y 1 según qué tan cerca están su temperatura y su precipitación de los rangos que necesita el cultivo.
2. Similaridad climática (enfoque de análogos): a cada celda se le mide la distancia euclidiana entre su clima mensual y el de sitios del Valle del Cauca donde la caña ya se cultiva con éxito.

Se siguen los pasos del ejemplo de la unidad (cargar las imágenes raster mes a mes, apilarlas, renombrar las capas y visualizarlas con una escala de color común). Como las librerías raster, rgdal y rasterVis del ejemplo fueron retiradas o dependen de paquetes retirados de CRAN, se usan sus reemplazos modernos: terra para los raster (rast() ≈ raster()/stack(), crop()/mask() iguales, terra::extract() igual) y sf para los vectores (st_read() ≈ readOGR()).

2 Librerías y datos

Los datos deben estar en una carpeta data_global/ junto a este archivo .Rmd, con una subcarpeta por variable (tmean/ y prec/). En tmean/ va el contenido de wc2.1_10m_tavg.zip y en prec/ el de wc2.1_10m_prec.zip, descargados de https://www.worldclim.org/data/worldclim21.html.

library(terra)     # lectura y álgebra de imágenes raster (reemplaza a raster)
library(sf)        # datos vectoriales: shape global de países (reemplaza a rgdal)
library(dplyr)     # manipulación de tablas
library(tidyr)     # reestructuración de tablas (ojo: tidyr::extract enmascara terra::extract)
library(ggplot2)   # mapas y gráficos
library(scales)    # formatos de ejes y leyendas
library(patchwork) # unir varios gráficos en una figura
library(knitr)     # tablas

ruta_datos <- "data_global"
ruta_tmean <- file.path(ruta_datos, "tmean")
ruta_prec  <- file.path(ruta_datos, "prec")
meses <- c("Ene", "Feb", "Mar", "Abr", "May", "Jun",
           "Jul", "Ago", "Sep", "Oct", "Nov", "Dic")

2.1 Imágenes raster de clima

Igual que en el ejemplo, se listan los 12 archivos de cada variable con list.files() y se cargan todos juntos. En terra, rast() con un vector de archivos hace lo mismo que stack().

# Lista los archivos raster (asc, tif o bil) y los ordena por el número de mes
listar_meses <- function(carpeta) {
  archivos <- list.files(carpeta, pattern = "\\.(asc|tif|bil)$", full.names = TRUE)
  mes <- as.integer(sub(".*_(\\d+)\\.[a-z]+$", "\\1", basename(archivos)))
  archivos[order(mes)]
}

tmean <- rast(listar_meses(ruta_tmean))
prec  <- rast(listar_meses(ruta_prec))

# Se asegura el sistema de coordenadas geográficas WGS84 (los .asc del ejemplo no lo traen)
crs(tmean) <- "EPSG:4326"
crs(prec)  <- "EPSG:4326"

# Se cambian los nombres de las capas por los meses (como names(...) = month.name)
names(tmean) <- meses
names(prec)  <- meses

# Se excluye la Antártida (al sur de 60° S), igual que la extensión de los datos del
# ejemplo: no tiene interés agrícola y sus temperaturas extremas distorsionan las estadísticas
tmean <- crop(tmean, ext(-180, 180, -60, 90))
prec  <- crop(prec,  ext(-180, 180, -60, 90))
tmean

## class       : SpatRaster
## size        : 900, 2160, 12  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s)   : memory
## names       :        Ene,        Feb,        Mar,        Abr,        May,       Jun, ...
## min values  : -45.883999, -43.172001, -42.453999, -33.623749, -21.366499,     -14.2, ...
## max values  :  34.009499,  32.824249,    32.9095,  34.193748,   36.25325, 38.355499, ...

En la versión 1.4 de WorldClim (la del ejemplo) la temperatura viene multiplicada por 10 para guardarla como entero; en la 2.1 ya viene en °C. Se revisa el rango y, si hace falta, se divide por 10.

rango_t <- global(tmean[[1]], range, na.rm = TRUE)
rango_t

##          X1      X2
## Ene -45.884 34.0095

if (max(rango_t) > 100) tmean <- tmean / 10   # de décimas de °C a °C
global(tmean, range, na.rm = TRUE)

##            X1       X2
## Ene -45.88400 34.00950
## Feb -43.17200 32.82425
## Mar -42.45400 32.90950
## Abr -33.62375 34.19375
## May -21.36650 36.25325
## Jun -14.20000 38.35550
## Jul -13.42875 39.54950
## Ago -13.99375 38.43275
## Sep -20.59575 35.79000
## Oct -29.69400 32.65125
## Nov -37.54250 32.78800
## Dic -45.32700 32.82525

global(prec,  range, na.rm = TRUE)

##     X1   X2
## Ene  0  908
## Feb  0  793
## Mar  0  720
## Abr  0 1004
## May  0 2068
## Jun  0 2210
## Jul  0 2381
## Ago  0 1674
## Sep  0 1955
## Oct  0 2328
## Nov  0  718
## Dic  0  806

La resolución de la grilla es de 0.1666667° (unos 18 km en el ecuador), así que cada celda resume el clima promedio de un área grande. Esto es importante cuando se extraiga el clima de puntos puntuales del Valle del Cauca.

2.2 Shape global de países

# Si el caso trae un shape global se usa ese; si no, el mapa de países de spData
shp <- list.files(ruta_datos, pattern = "\\.shp$", recursive = TRUE, full.names = TRUE)
if (length(shp) > 0) {
  paises <- st_read(shp[1], quiet = TRUE)
} else {
  paises <- spData::world
}
paises <- st_transform(paises, 4326)
# Columna con el nombre del país (cambia según la fuente del shape)
col_nombre <- intersect(c("name_long", "NAME", "NAME_LONG", "name", "COUNTRY", "ADMIN"),
                        names(paises))[1]
paises$pais <- paises[[col_nombre]]
paises <- paises[paises$pais != "Antarctica", ]   # misma extensión que los raster
nrow(paises)

## [1] 176

2.3 Funciones de apoyo para graficar

ggplot2 no grafica objetos SpatRaster directamente, así que se pasan a tabla con as.data.frame(xy = TRUE). Todas las figuras usan la misma escala de color por variable para poder comparar los meses entre sí.

# Raster -> tabla larga (x, y, capa, valor) para ggplot.
# Solo para graficar, los raster de más de 0.25° se agregan (promedio de 3x3 celdas)
# para que las figuras sean livianas; los cálculos se hacen con la resolución completa.
raster_a_tabla <- function(r, agregar = res(r)[1] < 0.25) {
  if (agregar) r <- aggregate(r, fact = 3, fun = "mean", na.rm = TRUE)
  as.data.frame(r, xy = TRUE, na.rm = TRUE) %>%
    pivot_longer(-c(x, y), names_to = "capa", values_to = "valor") %>%
    mutate(capa = factor(capa, levels = names(r)))
}

# Mapa de un raster continuo con el contorno de los países
mapa_continuo <- function(r, titulo, leyenda, paleta = "viridis", direccion = 1,
                          limites = NULL, contorno = paises,
                          agregar = res(r)[1] < 0.25) {
  ggplot(raster_a_tabla(r, agregar)) +
    geom_raster(aes(x, y, fill = valor)) +
    {if (!is.null(contorno)) geom_sf(data = contorno, fill = NA, colour = "grey30",
                                     linewidth = 0.1, inherit.aes = FALSE)} +
    scale_fill_viridis_c(option = paleta, direction = direccion, limits = limites,
                         name = leyenda) +
    coord_sf(expand = FALSE) +
    labs(title = titulo, x = NULL, y = NULL) +
    theme_minimal(base_size = 10)
}

2.4 Exploración de la temperatura y la precipitación mensual

ggplot(raster_a_tabla(tmean)) +
  geom_raster(aes(x, y, fill = valor)) +
  facet_wrap(~ capa, ncol = 4) +
  scale_fill_distiller(palette = "RdBu", name = "°C") +   # misma escala para los 12 meses
  coord_equal(expand = FALSE) +
  labs(title = "Temperatura media mensual (línea base)", x = NULL, y = NULL) +
  theme_void(base_size = 10)

ggplot(raster_a_tabla(prec)) +
  geom_raster(aes(x, y, fill = pmin(valor, 500))) +
  facet_wrap(~ capa, ncol = 4) +
  scale_fill_distiller(palette = "YlGnBu", direction = 1,
                       name = "mm\n(≥500\nagrupado)") +
  coord_equal(expand = FALSE) +
  labs(title = "Precipitación mensual (línea base)", x = NULL, y = NULL) +
  theme_void(base_size = 10)

Para el análisis de aptitud se resumen los 12 meses en variables anuales: la temperatura media anual, la precipitación total anual y la temperatura del mes más frío.

t_anual   <- mean(tmean)   # temperatura media anual (°C)
p_anual   <- sum(prec)     # precipitación total anual (mm)
t_mes_frio <- min(tmean)   # temperatura del mes más frío (°C)
names(t_anual) <- "t_anual"; names(p_anual) <- "p_anual"; names(t_mes_frio) <- "t_frio"

mapa_continuo(t_anual, "Temperatura media anual", "°C", paleta = "inferno") +
  mapa_continuo(clamp(p_anual, upper = 4000), "Precipitación anual", "mm\n(≥4000\nagrupado)",
                paleta = "mako", direccion = -1) +
  plot_layout(ncol = 1)

3 Punto 1. Mapas de aptitud climática para la caña de azúcar

3.1 Rangos óptimos del cultivo

Se usan los requerimientos de la base de datos ECOCROP de la FAO para Saccharum officinarum. Para cada variable hay un rango absoluto (fuera de él el cultivo no prospera) y un rango óptimo (donde el cultivo tiene su mejor desempeño):

rangos <- tibble::tribble(
  ~variable,                        ~abs_min, ~opt_min, ~opt_max, ~abs_max,
  "Temperatura media anual (°C)",         15,       24,       37,       41,
  "Precipitación anual (mm)",           1000,     1500,     2500,     5000,
  "Temperatura del mes más frío (°C)",    15,       18,       99,       99   # sin heladas ni frío prolongado
)
kable(rangos, caption = "Rangos de temperatura y precipitación para la caña de azúcar (FAO ECOCROP)")

Rangos de temperatura y precipitación para la caña de azúcar (FAO ECOCROP)

variable	abs_min	opt_min	opt_max	abs_max
Temperatura media anual (°C)	15	24	37	41
Precipitación anual (mm)	1000	1500	2500	5000
Temperatura del mes más frío (°C)	15	18	99	99

La caña es una planta C4 de clima tropical: crece mejor con temperaturas altas, necesita mucha agua durante el crecimiento y no tolera heladas. Por eso se agrega la temperatura del mes más frío: con menos de 15 °C la caña prácticamente deja de crecer.

3.2 Índice de aptitud (función trapezoidal)

Como en el modelo ECOCROP, la aptitud de cada variable es una función trapezoidal: vale 0 fuera del rango absoluto, 1 dentro del rango óptimo y crece o decrece linealmente entre ambos. La aptitud climática total es el mínimo de las aptitudes parciales (ley del mínimo de Liebig: el factor más limitante define la aptitud).

# Aptitud trapezoidal: 0 fuera de [abs_min, abs_max], 1 en [opt_min, opt_max]
aptitud_trapecio <- function(r, abs_min, opt_min, opt_max, abs_max) {
  subida <- (r - abs_min) / (opt_min - abs_min)   # tramo creciente
  bajada <- (abs_max - r) / (abs_max - opt_max)   # tramo decreciente
  apt <- min(c(subida, bajada))                   # el menor de los dos tramos
  clamp(apt, lower = 0, upper = 1)                # se acota entre 0 y 1
}

apt_t    <- aptitud_trapecio(t_anual,    15,   24,   37,   41)
apt_p    <- aptitud_trapecio(p_anual,  1000, 1500, 2500, 5000)
apt_frio <- aptitud_trapecio(t_mes_frio, 15,   18,   99,   99)

aptitud <- min(c(apt_t, apt_p, apt_frio))   # ley del mínimo
names(apt_t) <- "Temperatura"; names(apt_p) <- "Precipitación"
names(apt_frio) <- "Mes más frío"; names(aptitud) <- "Aptitud climática"

parciales <- c(apt_t, apt_p, apt_frio)
ggplot(raster_a_tabla(parciales)) +
  geom_raster(aes(x, y, fill = valor)) +
  facet_wrap(~ capa, ncol = 1) +
  scale_fill_distiller(palette = "YlGn", direction = 1, limits = c(0, 1),
                       name = "Aptitud") +
  coord_equal(expand = FALSE) +
  labs(title = "Aptitud parcial por variable climática", x = NULL, y = NULL) +
  theme_void(base_size = 10)

Para leer el mapa final, el índice se agrupa en cuatro clases:

clases <- c("No apta", "Marginal", "Moderada", "Alta")
aptitud_cl <- classify(aptitud, rcl = matrix(c(-Inf, 0,    0,
                                                0,   0.5,  1,
                                                0.5, 0.8,  2,
                                                0.8, Inf,  3), ncol = 3, byrow = TRUE),
                       include.lowest = TRUE, right = TRUE)
levels(aptitud_cl) <- data.frame(id = 0:3, clase = clases)
colores_clase <- c("No apta" = "grey85", "Marginal" = "#fee08b",
                   "Moderada" = "#a6d96a", "Alta" = "#1a9641")

# Para graficar se agrega con la moda (clase más frecuente en cada bloque de 3x3)
tabla_aptitud <- as.data.frame(aggregate(aptitud_cl, fact = 3, fun = "modal", na.rm = TRUE),
                               xy = TRUE, na.rm = TRUE)

ggplot(tabla_aptitud) +
  geom_raster(aes(x, y, fill = clase)) +
  geom_sf(data = paises, fill = NA, colour = "grey40", linewidth = 0.1,
          inherit.aes = FALSE) +
  scale_fill_manual(values = colores_clase, name = "Aptitud") +
  coord_sf(expand = FALSE) +
  labs(title = "Aptitud climática para la caña de azúcar (rangos óptimos FAO ECOCROP)",
       x = NULL, y = NULL) +
  theme_minimal(base_size = 10)

Se calcula el área de cada clase con cellSize(), que tiene en cuenta que las celdas son más pequeñas cerca de los polos.

area_celda <- cellSize(aptitud_cl, unit = "km")
area_clase <- zonal(area_celda, aptitud_cl, fun = "sum", na.rm = TRUE)
names(area_clase) <- c("clase", "area_km2")
area_clase %>%
  mutate(area_millones_km2 = round(area_km2 / 1e6, 2),
         porcentaje = percent(area_km2 / sum(area_km2), accuracy = 0.1)) %>%
  select(-area_km2) %>%
  kable(caption = "Superficie continental por clase de aptitud")

Superficie continental por clase de aptitud

clase	area_millones_km2	porcentaje
No apta	111.16	78.5%
Marginal	7.37	5.2%
Moderada	6.13	4.3%
Alta	16.89	11.9%

Para saber qué limita la aptitud en cada lugar, se identifica la variable con el menor índice parcial en las celdas no aptas o marginales:

limitante <- which.min(parciales)                       # capa con la aptitud más baja
limitante <- mask(limitante, aptitud > 0.5, maskvalues = 1)  # solo donde la aptitud es baja
levels(limitante) <- data.frame(id = 1:3, factor = names(parciales))

ggplot(as.data.frame(aggregate(limitante, fact = 3, fun = "modal", na.rm = TRUE),
                     xy = TRUE, na.rm = TRUE)) +
  geom_raster(aes(x, y, fill = factor)) +
  geom_sf(data = paises, fill = NA, colour = "grey40", linewidth = 0.1,
          inherit.aes = FALSE) +
  scale_fill_manual(values = c("Temperatura" = "#d7301f", "Precipitación" = "#2b8cbe",
                               "Mes más frío" = "#756bb1"), name = "Factor limitante") +
  coord_sf(expand = FALSE) +
  labs(title = "Factor que más limita la aptitud (celdas con aptitud ≤ 0.5)",
       x = NULL, y = NULL) +
  theme_minimal(base_size = 10)

4 Punto 2. Países con alto potencial y corte con el shape global

Para escoger los países de forma objetiva se cruza el mapa de aptitud con el shape de países: se rasteriza el shape sobre la misma grilla del clima y se suma el área con aptitud alta en cada país.

paises_r <- rasterize(vect(paises), aptitud, field = "pais")      # país de cada celda
area_alta <- area_celda * (aptitud > 0.8)                          # km² con aptitud alta
ranking <- zonal(c(area_alta, area_celda), paises_r, fun = "sum", na.rm = TRUE)
names(ranking) <- c("pais", "area_alta_km2", "area_total_km2")
ranking <- ranking %>%
  mutate(pct_alta = area_alta_km2 / area_total_km2) %>%
  arrange(desc(area_alta_km2))

ranking %>%
  head(10) %>%
  mutate(area_alta_miles_km2 = round(area_alta_km2 / 1000),
         pct_alta = percent(pct_alta, accuracy = 0.1)) %>%
  select(pais, area_alta_miles_km2, pct_alta) %>%
  kable(caption = "Diez países con más superficie de aptitud climática alta")

Diez países con más superficie de aptitud climática alta

pais	area_alta_miles_km2	pct_alta
Brazil	5152	60.6%
Democratic Republic of the Congo	1665	71.7%
Indonesia	1059	58.3%
Colombia	645	55.9%
Peru	538	41.1%
Venezuela	535	58.9%
Bolivia	401	37.0%
Cameroon	331	72.0%
Central African Republic	315	50.5%
Republic of the Congo	297	87.5%

Se toman los dos países con más área de aptitud alta y se agrega Colombia, que es el país de referencia del caso.

paises_sel <- unique(c(head(ranking$pais, 2), "Colombia"))
paises_sel

## [1] "Brazil"                           "Democratic Republic of the Congo"
## [3] "Colombia"

El corte se hace en dos pasos, como con crop y mask del paquete raster: crop() recorta el raster a la caja que contiene el país y mask() deja en NA las celdas que caen fuera de su polígono.

mapa_pais <- function(nombre) {
  pol <- paises[paises$pais == nombre, ]
  r <- mask(crop(aptitud, vect(pol)), vect(pol))     # corte con el shape del país
  ggplot(as.data.frame(r, xy = TRUE, na.rm = TRUE)) +
    geom_raster(aes(x, y, fill = `Aptitud climática`)) +
    geom_sf(data = pol, fill = NA, colour = "black", linewidth = 0.3, inherit.aes = FALSE) +
    scale_fill_distiller(palette = "YlGn", direction = 1, limits = c(0, 1),
                         name = "Aptitud") +
    coord_sf(expand = FALSE) +
    labs(title = nombre, x = NULL, y = NULL) +
    theme_minimal(base_size = 9)
}
wrap_plots(lapply(paises_sel, mapa_pais), nrow = 1, guides = "collect")

5 Punto 3. Clima de sitios en el Valle del Cauca

5.1 Geocodificación de los sitios

Se escogieron tres sitios del Valle del Cauca con ayuda de Google Maps, buscando que caigan en celdas distintas de la grilla y que representen ambientes diferentes:

• Palmira: centro de la zona azucarera del valle geográfico del río Cauca.
• Cartago: extremo norte del valle, también con cultivos de caña.
• Bajo Calima (Buenaventura): llanura del Pacífico, un clima muy lluvioso donde no se cultiva caña comercial.

Las coordenadas se tomaron de Google Maps. Como en el ejemplo de la unidad, también se pueden obtener geocodificando el nombre con geocode_OSM() de tmaptools (el chunk solo se ejecuta si la librería está instalada y hay conexión).

sitios <- tibble::tribble(
  ~sitio,                        ~lon,     ~lat,
  "Palmira",                  -76.303,    3.539,
  "Cartago",                  -75.912,    4.747,
  "Bajo Calima (Buenaventura)", -76.970,  3.990
)

library(tmaptools)
geo <- tryCatch(geocode_OSM(c("Palmira, Valle del Cauca, Colombia",
                              "Cartago, Valle del Cauca, Colombia")),
                error = function(e) NULL)
geo

##                                query      lat       lon  lat_min  lat_max
## 1 Palmira, Valle del Cauca, Colombia 3.530837 -76.29880 3.451623 3.718068
## 2 Cartago, Valle del Cauca, Colombia 4.710659 -75.93239 4.612901 4.808081
##     lon_min   lon_max
## 1 -76.49174 -75.94967
## 2 -76.01666 -75.82212

sitios_sf <- st_as_sf(sitios, coords = c("lon", "lat"), crs = 4326)
sitios$celda <- cellFromXY(tmean, as.matrix(sitios[, c("lon", "lat")]))
sitios   # se verifica que cada sitio caiga en una celda diferente

## # A tibble: 3 × 4
##   sitio                        lon   lat   celda
##   <chr>                      <dbl> <dbl>   <dbl>
## 1 Palmira                    -76.3  3.54 1119503
## 2 Cartago                    -75.9  4.75 1104385
## 3 Bajo Calima (Buenaventura) -77.0  3.99 1115179

valle <- st_as_sfc(st_bbox(c(xmin = -78, xmax = -75, ymin = 2.5, ymax = 5.5), crs = 4326))
mapa_continuo(crop(t_anual, vect(valle)), "Sitios sobre la temperatura media anual",
              "°C", paleta = "inferno", contorno = NULL, agregar = FALSE) +
  geom_sf(data = sitios_sf, colour = "cyan", size = 3, inherit.aes = FALSE) +
  geom_sf_label(data = sitios_sf, aes(label = sitio), nudge_y = 0.25, size = 3,
                inherit.aes = FALSE)

5.2 Extracción del clima y series mensuales

Con terra::extract() se obtiene el valor de las 12 capas de temperatura y de precipitación en cada sitio.

coords <- as.matrix(sitios[, c("lon", "lat")])
clima_sitios <- bind_rows(
  data.frame(sitio = sitios$sitio, variable = "tmean", terra::extract(tmean, coords)),
  data.frame(sitio = sitios$sitio, variable = "prec",  terra::extract(prec,  coords))
) %>%
  pivot_longer(all_of(meses), names_to = "mes", values_to = "valor") %>%
  mutate(mes = factor(mes, levels = meses)) %>%
  pivot_wider(names_from = variable, values_from = valor)

clima_sitios %>%
  group_by(sitio) %>%
  summarise(`T media anual (°C)` = round(mean(tmean), 1),
            `T mes más frío (°C)` = round(min(tmean), 1),
            `Precipitación anual (mm)` = round(sum(prec))) %>%
  kable(caption = "Resumen climático de los sitios")

Resumen climático de los sitios

sitio	T media anual (°C)	T mes más frío (°C)	Precipitación anual (mm)
Bajo Calima (Buenaventura)	26.1	25.8	5631
Cartago	23.1	22.5	1896
Palmira	22.4	21.9	1484

g_t <- ggplot(clima_sitios, aes(mes, tmean, colour = sitio, group = sitio)) +
  annotate("rect", xmin = -Inf, xmax = Inf, ymin = 24, ymax = 37,
           fill = "darkgreen", alpha = 0.08) +           # rango óptimo de temperatura
  geom_line(linewidth = 1) + geom_point() +
  labs(title = "Temperatura media mensual", x = NULL, y = "°C", colour = NULL,
       caption = "Franja verde: rango óptimo de temperatura (24 a 37 °C)") +
  theme_minimal(base_size = 10)

g_p <- ggplot(clima_sitios, aes(mes, prec, fill = sitio)) +
  geom_col(position = "dodge") +
  labs(title = "Precipitación mensual", x = NULL, y = "mm", fill = NULL) +
  theme_minimal(base_size = 10)

g_t / g_p + plot_layout(guides = "collect")

6 Punto 4. Mapas de similaridad climática

6.1 Métrica: distancia euclidiana estandarizada

El clima de cada celda se describe con un vector de 24 variables (12 temperaturas y 12 precipitaciones mensuales). Antes de medir la distancia se hacen dos ajustes:

• La precipitación mensual es muy asimétrica (muchas celdas secas y pocas con más de 1.000 mm), así que se transforma con \(\log(1 + P)\). Así, una diferencia de 100 mm pesa más entre un mes seco y uno húmedo que entre dos meses muy lluviosos.
• Como la temperatura (°C) y la precipitación (mm) tienen escalas muy distintas, cada capa se estandariza con su media y desviación estándar global (puntaje z). Así ambas variables pesan lo mismo en la distancia:

\[d(c, s) = \sqrt{\sum_{k=1}^{24} \left(z_{k}(c) - z_{k}(s)\right)^2}\]

donde \(z_k(c)\) es el valor estandarizado de la variable \(k\) en la celda \(c\) y \(z_k(s)\) el del sitio de referencia \(s\). Una distancia pequeña significa que la celda tiene un clima análogo al del sitio, con el mismo ciclo estacional. Para facilitar la lectura, la distancia se convierte en una similaridad entre 0 y 1: \(sim = 1 - d / d_{max}\).

clima <- c(tmean, log1p(prec))     # precipitación en escala log(1 + P)
names(clima) <- c(paste0("t_", meses), paste0("p_", meses))
medias <- global(clima, "mean", na.rm = TRUE)$mean
desv   <- global(clima, "sd",   na.rm = TRUE)$sd
clima_z <- (clima - medias) / desv            # puntaje z de cada capa
z_sitios <- terra::extract(clima_z, coords)          # vector z de cada sitio

similaridad_sitio <- function(i) {
  d <- sqrt(sum((clima_z - unlist(z_sitios[i, ]))^2))   # distancia euclidiana
  1 - d / global(d, "max", na.rm = TRUE)$max             # similaridad 0-1
}
similaridad <- rast(lapply(seq_len(nrow(sitios)), similaridad_sitio))

## |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          

names(similaridad) <- sitios$sitio

ggplot(raster_a_tabla(similaridad)) +
  geom_raster(aes(x, y, fill = valor)) +
  geom_sf(data = paises, fill = NA, colour = "grey40", linewidth = 0.1,
          inherit.aes = FALSE) +
  geom_sf(data = sitios_sf, colour = "red", size = 1.5, inherit.aes = FALSE) +
  facet_wrap(~ capa, ncol = 1) +
  scale_fill_viridis_c(option = "magma", name = "Similaridad") +
  coord_sf(expand = FALSE) +
  labs(title = "Similaridad climática con cada sitio del Valle del Cauca",
       x = NULL, y = NULL) +
  theme_minimal(base_size = 10)

La similaridad decae de forma continua, así que para hablar de zonas análogas se toma el 2 % de las celdas continentales más parecidas a cada sitio (unos 2,5 millones de km², un área comparable a la de los grandes países productores):

umbral <- global(similaridad, quantile, probs = 0.98, na.rm = TRUE)[, 1]
analogos <- similaridad >= umbral        # TRUE en el 2 % de celdas más similares
analogos <- mask(analogos, analogos, maskvalues = 0)

# Para graficar: 1 = análogo, 0 = resto del continente (máximo en bloques de 3x3)
analogos_graf <- mask(classify(analogos, cbind(NA, 0)), t_anual)
analogos_graf <- aggregate(analogos_graf, fact = 3, fun = "max", na.rm = TRUE)
names(analogos_graf) <- sitios$sitio

ggplot(raster_a_tabla(analogos_graf, agregar = FALSE)) +
  geom_raster(aes(x, y, fill = factor(valor, levels = 0:1, labels = c("No", "Sí")))) +
  geom_sf(data = paises, fill = NA, colour = "grey50", linewidth = 0.1,
          inherit.aes = FALSE) +
  geom_sf(data = sitios_sf, colour = "black", size = 1.5, inherit.aes = FALSE) +
  scale_fill_manual(values = c("No" = "grey90", "Sí" = "#d95f02"), name = "Análogo") +
  facet_wrap(~ capa, ncol = 1) +
  coord_sf(expand = FALSE) +
  labs(title = "Zonas análogas: 2 % de celdas con clima más parecido a cada sitio",
       x = NULL, y = NULL) +
  theme_minimal(base_size = 10)

area_analogos <- zonal(area_celda * classify(analogos, cbind(NA, 0)), paises_r,
                       fun = "sum", na.rm = TRUE)
names(area_analogos) <- c("pais", sitios$sitio)
area_analogos %>%
  pivot_longer(-pais, names_to = "sitio", values_to = "km2") %>%
  group_by(sitio) %>%
  slice_max(km2, n = 5) %>%
  mutate(miles_km2 = round(km2 / 1000)) %>%
  select(sitio, pais, miles_km2) %>%
  kable(caption = "Cinco países con más área análoga a cada sitio")

Cinco países con más área análoga a cada sitio

sitio	pais	miles_km2
Bajo Calima (Buenaventura)	Indonesia	983
Bajo Calima (Buenaventura)	Brazil	510
Bajo Calima (Buenaventura)	Colombia	491
Bajo Calima (Buenaventura)	Peru	297
Bajo Calima (Buenaventura)	Malaysia	252
Cartago	Democratic Republic of the Congo	981
Cartago	Brazil	495
Cartago	Indonesia	440
Cartago	Peru	339
Cartago	Colombia	314
Palmira	Democratic Republic of the Congo	1095
Palmira	Brazil	606
Palmira	Peru	334
Palmira	Indonesia	212
Palmira	Republic of the Congo	212

7 Punto 5. Comparación de las dos aproximaciones

Se compara la aptitud alta o moderada (índice > 0.5) con las zonas análogas a Palmira, que representa la zona azucarera del Valle del Cauca. Cada celda queda en una de cuatro categorías:

apto <- aptitud > 0.5
analogo_palmira <- classify(analogos[["Palmira"]], cbind(NA, 0)) == 1
analogo_palmira <- mask(analogo_palmira, aptitud)

comparacion <- apto * 1 + analogo_palmira * 2   # 0 ninguno, 1 solo aptitud, 2 solo análogo, 3 ambos
levels(comparacion) <- data.frame(id = 0:3,
  categoria = c("Ninguno", "Solo aptitud (ECOCROP)", "Solo análogo a Palmira", "Ambos"))

ggplot(as.data.frame(aggregate(comparacion, fact = 3, fun = "modal", na.rm = TRUE),
                     xy = TRUE, na.rm = TRUE)) +
  geom_raster(aes(x, y, fill = categoria)) +
  geom_sf(data = paises, fill = NA, colour = "grey40", linewidth = 0.1,
          inherit.aes = FALSE) +
  scale_fill_manual(values = c("Ninguno" = "grey90", "Solo aptitud (ECOCROP)" = "#a6d96a",
                               "Solo análogo a Palmira" = "#fdae61", "Ambos" = "#1a9641"),
                    name = NULL) +
  coord_sf(expand = FALSE) +
  labs(title = "Coincidencia entre aptitud por rangos óptimos y similaridad con Palmira",
       x = NULL, y = NULL) +
  theme_minimal(base_size = 10) +
  theme(legend.position = "bottom")

area_comp <- zonal(area_celda, comparacion, fun = "sum", na.rm = TRUE)
names(area_comp) <- c("categoria", "area_km2")
area_comp %>%
  mutate(miles_km2 = round(area_km2 / 1000),
         porcentaje = percent(area_km2 / sum(area_km2), accuracy = 0.1)) %>%
  select(-area_km2) %>%
  kable(caption = "Área en cada categoría de coincidencia")

Área en cada categoría de coincidencia

categoria	miles_km2	porcentaje
Ninguno	118021	83.4%
Solo aptitud (ECOCROP)	19606	13.9%
Solo análogo a Palmira	506	0.4%
Ambos	3413	2.4%

a_ambos  <- area_comp$area_km2[area_comp$categoria == "Ambos"]
a_apto   <- sum(area_comp$area_km2[area_comp$categoria %in% c("Ambos", "Solo aptitud (ECOCROP)")])
a_analog <- sum(area_comp$area_km2[area_comp$categoria %in% c("Ambos", "Solo análogo a Palmira")])
coincidencia <- tibble::tibble(
  indicador = c("% de las zonas análogas a Palmira que son aptas",
                "% de las zonas aptas que son análogas a Palmira",
                "Índice de Jaccard (área común / área unida)"),
  valor = c(percent(a_ambos / a_analog, 0.1), percent(a_ambos / a_apto, 0.1),
            round(a_ambos / (a_apto + a_analog - a_ambos), 2)))
kable(coincidencia)

indicador	valor
% de las zonas análogas a Palmira que son aptas	87.1%
% de las zonas aptas que son análogas a Palmira	14.8%
Índice de Jaccard (área común / área unida)	0.15

Las dos aproximaciones también se pueden comparar celda a celda con la correlación de Spearman entre el índice de aptitud y la similaridad con cada sitio:

valores <- as.data.frame(c(aptitud, similaridad), na.rm = TRUE)
sapply(sitios$sitio, function(s) round(cor(valores[[1]], valores[[s]],
                                           method = "spearman"), 2))

##                    Palmira                    Cartago 
##                       0.52                       0.53 
## Bajo Calima (Buenaventura) 
##                       0.55

8 Interpretación de resultados

Antes de interpretar, hay que tener en cuenta qué representan los datos. WorldClim 2.1 es el clima promedio de 1970-2000 en celdas de unos 18 km. El modelo de aptitud solo considera el clima (no el suelo, la pendiente, el riego ni la infraestructura) y supone un cultivo de secano. Por eso los mapas muestran el potencial climático y no dónde se cultiva caña hoy.

Aptitud global (punto 1). Solo el 11,9 % de la superficie continental (16,9 millones de km²) tiene aptitud climática alta para la caña. Otro 4,3 % es moderada y el 5,2 % marginal. El 78,5 % restante no es apto. Las zonas de aptitud alta forman un cinturón en el trópico húmedo: la Amazonía y el centro de Brasil, la cuenca del Congo y el golfo de Guinea, el sudeste asiático e Indonesia, Centroamérica y buena parte de Colombia. El mapa de factor limitante explica por qué el resto del mundo queda por fuera:

• La temperatura limita en 64,7 millones de km²: latitudes altas y zonas de montaña (Andes, Himalaya, Tíbet, el altiplano de África oriental).
• La precipitación limita en 47,9 millones de km²: los desiertos, pero también las sabanas y regiones monzónicas con estación seca larga, como el noreste de Brasil, la mayor parte de la India, el norte de Australia y el sur de África.
• El mes más frío limita en 5,9 millones de km², en los bordes subtropicales: sureste de Estados Unidos, sur de China, sur de Brasil y Paraguay.

Un resultado llamativo es que varios de los grandes productores mundiales de caña salen con poca aptitud alta. La India tiene solo el 9,3 % de su territorio en esa clase, China el 0,3 %, Australia el 1,5 %, y Pakistán y Sudáfrica el 0 %. En casi todos los casos el factor limitante es la lluvia. Esto no es un error del mapa: esos países cultivan la caña con riego, y el modelo de secano no lo tiene en cuenta.

Países con alto potencial (punto 2). Según el área con aptitud alta, los tres primeros son:

• Brasil: 5,15 millones de km², el 60,6 % de su territorio. La aptitud alta cubre la Amazonía y el centro-oeste. El noreste semiárido no es apto por falta de lluvia y el sur no lo es por el frío. Es coherente con que Brasil sea el mayor productor de caña del mundo.
• República Democrática del Congo: 1,67 millones de km², el 71,7 % de su territorio. Tiene un potencial climático enorme y casi sin aprovechar. Esto muestra que el clima es una condición necesaria, pero no suficiente, para tener una industria azucarera.
• Colombia: 645 mil km², el 55,9 % del país, cuarto lugar en el ranking mundial. El corte muestra aptitud alta en la Orinoquía, la Amazonía, el Magdalena Medio y el Caribe. Los Andes no son aptos por temperatura y el Pacífico no lo es por exceso de lluvia. El valle del río Cauca aparece como una franja angosta de aptitud moderada a alta entre las dos cordilleras.

Sitios del Valle del Cauca (punto 3). Las series mensuales muestran lo que se espera en un clima ecuatorial:

• La temperatura casi no cambia en el año: varía menos de 1 °C entre meses en los tres sitios.
• La estacionalidad la marca la lluvia. En Palmira y Cartago el régimen es bimodal, con lluvias en marzo-mayo y en septiembre-noviembre y dos temporadas secas: diciembre-febrero y junio-agosto, cuando Palmira baja a 60-70 mm al mes.

Palmira (22,4 °C; 1.484 mm al año) y Cartago (23,1 °C; 1.896 mm) tienen aptitud alta (0,82 y 0,90). Aun así, Palmira queda un poco por debajo del óptimo de temperatura (24 °C) y de lluvia (1.500 mm). Esto concuerda con que la zona azucarera del valle complementa con riego en las temporadas secas. Bajo Calima es más cálido (26,1 °C), pero recibe 5.631 mm al año, más que el máximo absoluto de 5.000 mm, por lo que su aptitud es 0. Hay que recordar que una celda de 18 km mezcla el fondo del valle con el piedemonte, por eso la temperatura de Palmira sale un poco más baja que la medida en las estaciones del valle.

Similaridad climática (punto 4). Los mapas continuos muestran el gran contraste entre los climas tropicales húmedos (similaridad alta) y los desiertos y zonas polares (baja). Las diferencias entre sitios se ven en las zonas análogas (el 2 % más similar):

• Palmira: sus análogos están en el centro de la cuenca del Congo (1,09 millones de km²), el sur de Brasil y Paraguay, el piedemonte amazónico de Perú y Bolivia, Uganda e islas de Indonesia. Todos son climas cálidos sin exceso de lluvia y con dos temporadas secas moderadas.
• Cartago: tiene un patrón parecido, pero con más área en Colombia, Perú e Indonesia.
• Bajo Calima: sus análogos son las selvas superhúmedas de Indonesia (Borneo, Sumatra), Malasia, Papúa y la Amazonía noroccidental. Son zonas muy distintas a las de los otros dos sitios.

Esto confirma que la métrica distingue bien un clima de valle interandino de uno de selva pluvial, aunque estén a menos de 80 km de distancia.

Comparación de las dos aproximaciones (punto 5).

• Las zonas análogas a Palmira están casi siempre dentro de la zona apta: el 87,1 % de ellas tiene aptitud mayor que 0,5 según ECOCROP. Las dos aproximaciones son consistentes.
• Lo contrario no se cumple: solo el 14,8 % del área apta es análoga a Palmira (índice de Jaccard de 0,15), y la correlación de Spearman entre los dos índices es moderada (0,52 a 0,55). ECOCROP identifica todo el clima que la caña tolera. La similaridad, en cambio, identifica solo los lugares que se parecen a un sitio en particular, incluido su ciclo de lluvias.
• Las zonas que coinciden (centro del Congo, sur de Brasil y Paraguay, piedemonte amazónico y algunas islas de Indonesia) son las de mayor confianza.
• Las zonas solo análogas son pocas (0,4 % del área). Están en bordes subtropicales y de montaña, donde el frío o la lluvia salen un poco de los rangos óptimos.
• El caso de Bajo Calima muestra el riesgo de usar solo la similaridad: sus análogos se parecen al sitio, pero no son aptos para la caña, porque el sitio de referencia tampoco lo es.

9 Conclusiones

1. El potencial climático de la caña de azúcar se concentra en el trópico húmedo. Solo el 11,9 % de la superficie continental tiene aptitud alta, y Brasil, la República Democrática del Congo, Indonesia y Colombia tienen las mayores áreas. Fuera del trópico la caña está limitada por la temperatura y el frío. Dentro del trópico, el factor que más limita es la lluvia.
2. El clima no explica por sí solo dónde se cultiva la caña. Varios grandes productores (India, China, Pakistán, Australia, Sudáfrica) tienen poca aptitud de secano y producen gracias al riego. En cambio, el Congo tiene un potencial enorme que no aprovecha. Un mapa de aptitud para tomar decisiones debe agregar riego, suelos, pendiente, acceso a mercados e infraestructura agroindustrial.
3. Colombia tiene un gran potencial climático (56 % de su territorio en aptitud alta), pero no en todas partes. El valle del río Cauca es una franja angosta con aptitud alta que está cerca del límite inferior de temperatura y lluvia. Su productividad depende del manejo (riego en las temporadas secas, variedades adaptadas). Las nuevas áreas con clima apto están sobre todo en la Orinoquía y el Caribe.
4. Las dos aproximaciones son complementarias. La aptitud por rangos óptimos responde dónde puede crecer la caña, pero depende de los rangos que se elijan. La similaridad responde dónde el clima se parece a una zona productiva conocida. No necesita parámetros del cultivo y tiene en cuenta el ciclo estacional, por lo que sirve para transferir variedades y prácticas del Valle del Cauca a otras regiones. Su debilidad es que depende del sitio de referencia y del umbral, y puede señalar análogos de un sitio no apto, como Bajo Calima. Donde ambas coinciden está la recomendación más sólida.
5. Hay varias limitaciones y pasos siguientes.
   • La resolución de 10 minutos mezcla el fondo del valle con las cordilleras. Para el análisis regional convendría usar la versión de 30 segundos (1 km) de WorldClim.
   • El clima de 1970-2000 debería compararse con escenarios de cambio climático (CMIP6) para saber si la aptitud del Valle del Cauca y de sus análogos se mantendrá.
   • El análisis se puede extender con variables de suelo y con la distribución real de cultivos de caña para validar los mapas.

10 Referencias

• FAO. ECOCROP: Saccharum officinarum L. Base de datos de requerimientos ecológicos de cultivos. https://gaez.fao.org/pages/ecocrop
• Hijmans, R. J., Guarino, L., Cruz, M., y Rojas, E. (2001). Computer tools for spatial analysis of plant genetic resources data: 1. DIVA-GIS. Plant Genetic Resources Newsletter, 127, 15-19 (modelo ECOCROP).
• Fick, S. E., y Hijmans, R. J. (2017). WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology, 37(12), 4302-4315.
• Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G., y Jarvis, A. (2005). Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology, 25, 1965-1978 (WorldClim).
• Arango-Londoño, D. (2023). Módulo 1 - Unidad 2. RPubs. https://rpubs.com/darango/M1U2