Caso Cultivo Caña de Azucar

Introducción

Una organización desea indentificar los rangos optimos para el cultivo de caña de azucar buscando las condiciones optimas que generen la mayor productividad es decir la mayor cantidad kg/ Ha de caña de azucar, se sabe que el clima y la precipitación son variabes muy importantes que impactan la productividad de estos cultivos.

Reto

El reto principal consiste en determinar los lugares optimos para cultivar caña de azucar, se sabe que la caña de azucar obtiene un mejor desempeño si se cultiva en zonas con las siguientes caractteristicas:

• Temperatua media entre 22.5 y 28 Grados Centigrados.
• Precipitación anual entre 1500 y 3000 mm.
• Precipitación mensual entre 125 y 290 mm.

1. Determinación de zonas óptimas para cultivo.

Para obtención de los datos nos apoyaremos en datos publicados en la web :

https://www.worldclim.org/ , en esta página tendremos datos de temperatura y precipitación a nivel global mes a mes.

Carga de librerias necesarias

require (terra)
require(geodata)
library(mapview)
library(sf)
library(rnaturalearth)
library(tmap)
library(tmaptools)      
library(dplyr)          
library(tidyr)          
library(leaflet)        
library(ggplot2) 
library(raster)

Carga del archivo de la base de datos

Descargamos los datos de la web para generamos los mapas temáticos identificando los lugares que cumplen las condiones optimas para el cultivo de caña de azucar:

# 1. DESCARGA de datos WorldClim (12 meses)

prec <- worldclim_global(var = "prec", res = 10, path = "data/")   # precipitación mensual (mm)
temp <- worldclim_global(var = "tavg", res = 10, path = "data/")  # temperatura media mensual (°C)

# 2. RENOMBRE de capas con abreviaturas de mes
meses <- month.abb   # "Jan","Feb",…,"Dec"
names(prec) <- meses
names(temp) <- meses

# 3. MÁSCARAS LÓGICAS

#  3.1. Temperatura entre 22.5 y 28 °C, mes a mes
temp_mask <- (temp >= 22.5) & (temp <= 28)
names(temp_mask) <- meses

#  3.2. Precipitación entre 125 y 290 mm, mes a mes
prec_mask_monthly <- (prec  >= 125) & (prec <= 290)
names(prec_mask_monthly) <- meses

#  3.3. Precipitación anual acumulada entre 1500 y 3500 mm
#      – suma las 12 capas de precipitación
annual_prec   <- app(prec, sum)            
names(annual_prec) <- "Annual_prec_mm"

#      – máscara sobre la suma anual
annual_prec_mask <- (annual_prec >= 1500) & (annual_prec <= 3500)
names(annual_prec_mask) <- "1500–3500 mm"

# 4. FUNCIÓN AUXILIAR para graficar máscaras mes a mes en 4 filas × 3 columnas
plot_monthly_mask <- function(mask_rast, titulo, paleta = c("white","steelblue")) {
  op <- par(
    mfrow = c(4, 3),     # 4 filas × 3 columnas
    mar   = c(2, 2, 2, 1),# márgenes interiores
    oma   = c(0, 0, 2, 0) # margen exterior para título global
  )
  on.exit(par(op), add = TRUE)
  
  for (i in seq_len(nlyr(mask_rast))) {
    plot(
      mask_rast[[i]], 
      col    = paleta, 
      legend = TRUE, 
      axes   = TRUE, 
      box    = TRUE,
      main   = names(mask_rast)[i]  # muestra "Jan","Feb",…
    )
  }
  mtext(titulo, outer = TRUE, cex = 1.2)
}

# 5. GRAFICAR
#  5.1. Temperatura mensual en rango
plot_monthly_mask(
  temp_mask,
  "Meses con temp. promedio entre 22.5–28 °C",
  paleta = c("white","tomato")
)

#  5.2. Precipitación mensual en rango
plot_monthly_mask(
  prec_mask_monthly,
  "Meses con precipitación entre 125–290 mm",
  paleta = c("white","navy")
)

#  5.3. Precipitación anual acumulada en rango (mapa único)
par(mfrow = c(1,1), mar = c(4,4,3,2))
plot(
  annual_prec_mask,
  col    = c("white","darkgreen"),
  legend = TRUE,
  main   = "Áreas con precipitación anual 1500–3500 mm"
)
legend(
  "bottomleft",
  legend = c("Fuera de rango","Dentro de rango"),
  fill   = c("white","darkgreen"),
  bg     = "white"
)

Generamos un mapa tematico identificando las zonas que cumplen con las tres condiciones optimas para los cultivos de caña de azucar:

# 6. CONDICIÓN COMPUESTA
# 6.1. Que cumpla temp Y precip **cada** mes
monthly_both <- temp_mask & prec_mask_monthly

# 6.2. Restricción: que cumpla eso en **todos** los meses
all_months <- app(monthly_both, fun=function(v) all(v == 1))
names(all_months) <- "AllMonthsOK"

# 6.3. Y además la precipitación anual
final_mask <- all_months & annual_prec_mask
names(final_mask) <- "Cumple_todas"

# 7. PLOT INTERACTIVO



# Asignar valores categóricos para mejor visualización
final_mask_factor <- classify(final_mask, matrix(c(0, NA, 1, 1), ncol = 2, byrow = TRUE))
levels(final_mask_factor) <- data.frame(ID = 1, Categoria = "Cumple_todas")

# Establecer paleta personalizada
mapviewOptions(fgb = FALSE)  # evitar fallos con visualización en algunos navegadores

# Visualización interactiva con zoom
mapview(
  final_mask_factor,
  col.regions = c("darkgreen"),
  na.color = "lightgrey",
  legend = TRUE,
  layer.name = "Zonas que cumplen las 3 condiciones"
)

2. Escogemos 3 paises (Colombia, Perú, Indonesia) y verificamos el potencial

# 1. Modo estático: incrusta mapas en el HTML
tmap_mode("plot")

# 2. Descargar rásteres globales WorldClim (12 capas = meses)

#prec <- worldclim_global(var = "prec", res = 10, path = "data/")
#temp <- worldclim_global(var = "tavg", res = 10, path = "data/")

# 3. Crear máscaras con terra::app()
# 3a. Temperatura mensual 22.5–28 °C
temp_mask  <- app(temp, fun = function(x) ifelse(x >= 22.5 & x <= 28, 1, NA))
# 4a. Precipitación mensual 125–290 mm
prec_mask  <- app(prec, fun = function(x) ifelse(x >= 125  & x <= 290, 1, NA))
# 4b. Precipitación anual 1500–3500 mm
anual_mask <- app(prec, fun = function(x) {
  tot <- sum(x, na.rm = TRUE)
  ifelse(tot >= 1500 & tot <= 3500, 1, NA)
})

Obtenemos y filtramos los poligonos de Colombia, Perú e Indonesia

# 5. Obtener y filtrar polígonos de Colombia, Perú e Indonesia
paises <- ne_countries(scale = "medium", returnclass = "sf")
sel3   <- subset(paises, name %in% c("Colombia", "Peru", "Indonesia"))

# 6. Generar los cuatro tipos de gráficos para cada país

for (i in seq_len(nrow(sel3))) {
  poly <- st_transform(sel3[i, ], crs = crs(temp))
  pais <- sel3$name[i]
  
  # 6.1 Panel 4×3 de temperatura mensual
  temp_plots <- vector("list", 12)
  for (m in 1:12) {
    t1 <- mask(crop(temp_mask[[m]], poly), poly)
    temp_plots[[m]] <- tm_shape(t1) +
      tm_raster(palette="Oranges", title = month.name[m]) +
      tm_layout(main.title = paste(pais, "Temp 22,5-28 °C", month.name[m]))
  }
  print(tmap_arrange(grobs = temp_plots, ncol = 4, nrow = 3))
  
  # 6.2 Panel 4×3 de precipitación mensual
  prec_plots <- vector("list", 12)
  for (m in 1:12) {
    p1 <- mask(crop(prec_mask[[m]], poly), poly)
    prec_plots[[m]] <- tm_shape(p1) +
      tm_raster(palette="Blues", title = month.name[m]) +
      tm_layout(main.title = paste(pais, "Prec 125-290 mm", month.name[m]))
  }
  print(tmap_arrange(grobs = prec_plots, ncol = 4, nrow = 3))
  
  # 6.3 Mapa único de precipitación anual
  a1 <- mask(crop(anual_mask, poly), poly)
  map_anual <- tm_shape(poly) + tm_borders() +
    tm_shape(a1) + tm_raster(palette="Greens", title="Precip Anual 1500-3500 mm") +
    tm_layout(main.title = paste(pais, "Precip Anual"))
  print(map_anual)
  
  # 6.4 Mapa de intersección: condiciones mensuales EN TODOS los meses + condición anual
  # 6.4a. Pila mensual donde 1 = cumple temp & prec ese mes
  monthly_intersect <- temp_mask & prec_mask
  # 6.4b. Reducir a un ráster único que vale 1 sólo si todas las capas == 1
  inter_all <- app(monthly_intersect, fun = function(x) {
    ifelse(all(x == 1, na.rm = FALSE), 1, NA)
  })
  # 6.4c. Recortar y enmascarar por el país
  inter_all <- mask(crop(inter_all, poly), poly)
  # 6.4d. Aplicar también la condición anual
  inter_all <- mask(inter_all, mask(crop(anual_mask, poly), poly))
  # 6.4e. Dibujar el mapa
  map_all <- tm_shape(poly) + tm_borders() +
    tm_shape(inter_all) + tm_raster(palette="Purples", title="Cumple todas condiciones") +
    tm_layout(main.title = paste(pais, 
      "— Temp & Prec en TODOS los meses + Precip Anual 1500-3500 mm"))
  print(map_all)
}

3. Se escogen 2 lugares del Valle y graficamos las series de tiempo

Para la ciudad de Palmira, Valle:

# 1. Crear un RasterBrick con raster::brick()
prec_stack <- brick(prec)   # 12 capas de precipitación
temp_stack <- brick(temp)   # 12 capas de temperatura

# 2. Geocodificar Palmira
gc <- geocode_OSM("Palmira, Valle")
coords_df <- data.frame(x = gc$coords[1], y = gc$coords[2])

# 3. Extraer los 12 valores mensuales
prec_vals <- raster::extract(prec_stack, coords_df)
temp_vals <- raster::extract(temp_stack, coords_df)

# 4. Verificar que tienes un vector de longitud 12
print(prec_vals)  # debería dar algo como:      [,1]  [,2] … [,12]

     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
[1,]  92 102 132 172 146 112  61  73 105 195 173 121

print(temp_vals)

          Jan     Feb      Mar      Apr     May    Jun     Jul      Aug
[1,] 22.50275 22.6155 22.73675 22.44175 22.3585 22.328 22.6225 22.74675
          Sep      Oct      Nov      Dec
[1,] 22.57575 22.02075 21.88275 22.14775

# 5. Convertir a serie y graficar
#library(dplyr); library(tidyr); library(ggplot2)

ts_df <- tibble(
  month = 1:12,
  prec  = as.numeric(prec_vals),
  temp  = as.numeric(temp_vals)
) %>%
  pivot_longer(c(prec, temp), names_to = "variable", values_to = "value")

ggplot(ts_df, aes(month, value)) +
  geom_line() +
  facet_wrap(~variable, scales = "free_y",
             labeller = labeller(variable = c(
               prec = "Precipitación (mm)",
               temp = "Temp. media (°C)"
             ))) +
  scale_x_continuous(breaks = 1:12, labels = month.abb) +
  theme_minimal() +
  labs(x = "Mes", y = NULL,
       title = "Serie climática mensual en Palmira, Valle")

Para la ciudad de Cartago, Valle:

# 1. Crear un RasterBrick con raster::brick()
prec_stack <- brick(prec)   # 12 capas de precipitación
temp_stack <- brick(temp)   # 12 capas de temperatura

# 2. Geocodificar Cartago
gc <- geocode_OSM("Cartago, Valle")
coords_df <- data.frame(x = gc$coords[1], y = gc$coords[2])

# 3. Extraer los 12 valores mensuales
prec_vals <- raster::extract(prec_stack, coords_df)
temp_vals <- raster::extract(temp_stack, coords_df)

# 4. Verificar que tienes un vector de longitud 12
print(prec_vals)  # debería dar algo como:      [,1]  [,2] … [,12]

     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
[1,] 106 109 150 217 202 159 110 129 169 241 197 107

print(temp_vals)

         Jan      Feb      Mar      Apr      May    Jun    Jul      Aug
[1,] 23.2375 23.45825 23.57375 23.26525 23.08175 23.006 23.382 23.42725
          Sep     Oct    Nov      Dec
[1,] 23.08925 22.5165 22.501 22.84375

# 5. Convertir a serie y graficar
#library(dplyr); library(tidyr); library(ggplot2)

ts_df <- tibble(
  month = 1:12,
  prec  = as.numeric(prec_vals),
  temp  = as.numeric(temp_vals)
) %>%
  pivot_longer(c(prec, temp), names_to = "variable", values_to = "value")

ggplot(ts_df, aes(month, value)) +
  geom_line() +
  facet_wrap(~variable, scales = "free_y",
             labeller = labeller(variable = c(
               prec = "Precipitación (mm)",
               temp = "Temp. media (°C)"
             ))) +
  scale_x_continuous(breaks = 1:12, labels = month.abb) +
  theme_minimal() +
  labs(x = "Mes", y = NULL,
       title = "Serie climática mensual en Cartago, Valle")

4. Determinación de mapas de similaridad para Palmira, Valle.

# 1. Descarga datos climáticos globales

#prec_stack <- worldclim_global(var = "prec", res = 10, path = "data/")
#temp_stack <- worldclim_global(var = "tavg", res = 10, path = "data/")
clim_global <- c(prec, temp)
#message("Datos climáticos cargados: ", nlyr(clim_global), " capas")

# 2. Coordenadas de Palmira
#message("Geocodificando Palmira...")
gc <- geocode_OSM("Palmira, Valle")
#if (is.null(gc$coords)) stop("No se encontraron coordenadas para Palmira.")
coords <- vect(data.frame(x = gc$coords[1], y = gc$coords[2]), 
               geom = c("x", "y"), crs = crs(clim_global))

# 3. Extracción de clima en Palmira
ref_vals <- terra::extract(clim_global, coords)[,-1]
ref_vec <- as.numeric(ref_vals)
#message("Valores extraídos en Palmira:")
print(ref_vec)

 [1]  92.00000 102.00000 132.00000 172.00000 146.00000 112.00000  61.00000
 [8]  73.00000 105.00000 195.00000 173.00000 121.00000  22.50275  22.61550
[15]  22.73675  22.44175  22.35850  22.32800  22.62250  22.74675  22.57575
[22]  22.02075  21.88275  22.14775

# 4. Estadísticas globales por capa
means_df <- global(clim_global, fun = "mean", na.rm = TRUE)
sds_df   <- global(clim_global, fun = "sd", na.rm = TRUE)

means <- as.numeric(means_df[,1])
sds   <- as.numeric(sds_df[,1])

# 5. Validar y limpiar
valid_idx <- which(!is.na(sds) & sds > 0 & !is.na(ref_vec))
#message("Capas válidas después de limpiar SDs cero o NA: ", length(valid_idx))
#if (length(valid_idx) < 6) stop("Muy pocas capas válidas para análisis")

# Aplicar filtros
clim_global <- clim_global[[valid_idx]]
means <- means[valid_idx]
sds <- sds[valid_idx]
ref_vec <- ref_vec[valid_idx]

# 6. Escalamiento
clim_std <- (clim_global - means) / sds
ref_std <- (ref_vec - means) / sds
#message("Escalamiento realizado")

# 7. Calcular distancias
message("Calculando distancias euclidianas...")
dist_map <- terra::app(clim_std, fun = function(cell_vals) {
  if (all(is.na(cell_vals))) return(NA)
  sqrt(sum((cell_vals - ref_std)^2, na.rm = TRUE))
})

# 8. Diagnóstico: visualización previa
vals <- values(dist_map)
vals <- vals[!is.na(vals)]
#message("Número de celdas válidas: ", length(vals))
#if (length(vals) == 0) stop("El mapa de distancias está vacío.")


# 9. Aplicar umbral (percentil 5%)
threshold <- quantile(vals, probs = 0.05, na.rm = TRUE)
message("Umbral (5%): ", threshold)

similar_zones <- dist_map < threshold

# 10. Verificar zonas encontradas
zone_vals <- values(similar_zones)
zone_vals <- zone_vals[!is.na(zone_vals)]
#message("Zonas similares encontradas: ", sum(zone_vals), " celdas")

# Filtrar solo celdas con zonas similares (valor TRUE)
pts_similares <- as.points(similar_zones, na.rm = TRUE)

# Verificar si hay puntos válidos
if (!is.null(pts_similares) && nrow(pts_similares) > 0) {
  ext_zoom <- ext(pts_similares)
} else {
  warning("No se encontraron zonas similares, se usará la extensión completa.")
  ext_zoom <- ext(dist_map)  # uso total como fallback
}

# Mapa de distancias con zoom
plot(dist_map, main = "Zonas climáticas similares a Palmira", col = hcl.colors(100, "Viridis"),
     xlim = c(ext_zoom[1], ext_zoom[2]), ylim = c(ext_zoom[3], ext_zoom[4]))

# Zonas similares sobrepuestas
plot(similar_zones, col = c(NA, "#00640088"), legend = FALSE, add = TRUE)

# Marcar Palmira
points(coords, col = "red", pch = 19, cex = 1.2)
text(coords, labels = "Palmira", pos = 4, col = "red")

5. Determinación de mapas de similaridad para Cartago, Valle.

# 1. Descarga datos climáticos globales

#prec_stack <- worldclim_global(var = "prec", res = 10, path = "data/")
#temp_stack <- worldclim_global(var = "tavg", res = 10, path = "data/")
clim_global <- c(prec, temp)
#message("Datos climáticos cargados: ", nlyr(clim_global), " capas")

# 2. Coordenadas de Cartago
#message("Geocodificando Cartago...")
gc <- geocode_OSM("Cartago, Valle")
#if (is.null(gc$coords)) stop("No se encontraron coordenadas para Palmira.")
coords <- vect(data.frame(x = gc$coords[1], y = gc$coords[2]), 
               geom = c("x", "y"), crs = crs(clim_global))

# 3. Extracción de clima en Cartago
ref_vals <- terra::extract(clim_global, coords)[,-1]
ref_vec <- as.numeric(ref_vals)
#message("Valores extraídos en Cartago:")
print(ref_vec)

 [1] 106.00000 109.00000 150.00000 217.00000 202.00000 159.00000 110.00000
 [8] 129.00000 169.00000 241.00000 197.00000 107.00000  23.23750  23.45825
[15]  23.57375  23.26525  23.08175  23.00600  23.38200  23.42725  23.08925
[22]  22.51650  22.50100  22.84375

# 4. Estadísticas globales por capa
means_df <- global(clim_global, fun = "mean", na.rm = TRUE)
sds_df   <- global(clim_global, fun = "sd", na.rm = TRUE)

means <- as.numeric(means_df[,1])
sds   <- as.numeric(sds_df[,1])

# 5. Validar y limpiar
valid_idx <- which(!is.na(sds) & sds > 0 & !is.na(ref_vec))
#message("Capas válidas después de limpiar SDs cero o NA: ", length(valid_idx))
#if (length(valid_idx) < 6) stop("Muy pocas capas válidas para análisis")

# Aplicar filtros
clim_global <- clim_global[[valid_idx]]
means <- means[valid_idx]
sds <- sds[valid_idx]
ref_vec <- ref_vec[valid_idx]

# 6. Escalamiento
clim_std <- (clim_global - means) / sds
ref_std <- (ref_vec - means) / sds
#message("Escalamiento realizado")

# 7. Calcular distancias
message("Calculando distancias euclidianas...")
dist_map <- terra::app(clim_std, fun = function(cell_vals) {
  if (all(is.na(cell_vals))) return(NA)
  sqrt(sum((cell_vals - ref_std)^2, na.rm = TRUE))
})

# 8. Diagnóstico: visualización previa
vals <- values(dist_map)
vals <- vals[!is.na(vals)]
#message("Número de celdas válidas: ", length(vals))
#if (length(vals) == 0) stop("El mapa de distancias está vacío.")


# 9. Aplicar umbral (percentil 5%)
threshold <- quantile(vals, probs = 0.05, na.rm = TRUE)
message("Umbral (5%): ", threshold)

similar_zones <- dist_map < threshold

# 10. Verificar zonas encontradas
zone_vals <- values(similar_zones)
zone_vals <- zone_vals[!is.na(zone_vals)]
#message("Zonas similares encontradas: ", sum(zone_vals), " celdas")

# Filtrar solo celdas con zonas similares (valor TRUE)
pts_similares <- as.points(similar_zones, na.rm = TRUE)

# Verificar si hay puntos válidos
if (!is.null(pts_similares) && nrow(pts_similares) > 0) {
  ext_zoom <- ext(pts_similares)
} else {
  warning("No se encontraron zonas similares, se usará la extensión completa.")
  ext_zoom <- ext(dist_map)  # uso total como fallback
}

# Mapa de distancias con zoom
plot(dist_map, main = "Zonas climáticas similares a Cartago", col = hcl.colors(100, "Viridis"),
     xlim = c(ext_zoom[1], ext_zoom[2]), ylim = c(ext_zoom[3], ext_zoom[4]))

# Zonas similares sobrepuestas
plot(similar_zones, col = c(NA, "#00640088"), legend = FALSE, add = TRUE)

# Marcar Cartago
points(coords, col = "red", pch = 19, cex = 1.2)
text(coords, labels = "Cartago", pos = 4, col = "red")

6. Conclusiones

La consulta de la base de datos de la web WordClimb resulta muy útil para poder tomar decisiones de busqueda de las zonas que cumplen determinadas condiciones, mediante el análisis de los parámetros de Temperatura y Precipitación se logra determinar las zonas del mundo que cumple con las condiciones ideales, de esta manera se pueden tomar decisiones releventes, por otro lado conocer las condiones de determinado lugar nos permite buscar zonas con condiciones similares según nuestros intereses.