Nota: Documento académico para la Maestría en ciencia de datos de la Universidad Javeriana de Cali

Introducción

Las instrucciones para realizar el trabajo se encuentran en el vídeo de YouTube titulado “[M1U2] Caso: situacion 5”:

En resumen, se solicita:

Generar un código en R utilizando datos de clima de línea base a nivel global para construir mapas de aptitud climática de la caña de azúcar, basándose en los rangos óptimos. Representar los mapas con una escala de colores adecuada.
Identificar 2 o 3 países con áreas de alto potencial para la caña de azúcar y hacer un recorte de estas zonas usando el shape global. Presentar los mapas con una escala de colores adecuada.
Escoger algunos puntos (2 o 3) al azar en la región del Valle del Cauca (usando Google Maps) y extraer la información climática. Graficar las series de tiempo de temperatura y precipitación.
Usar una métrica de similaridad, como la distancia euclidiana, para desarrollar un código en R que genere mapas de similaridad a nivel global para los sitios identificados en el punto anterior. Visualizar estos mapas con una escala de colores adecuada.
Comparar los mapas generados por ambas aproximaciones y presentar conclusiones.

Según la información del video, los rangos óptimos de temperatura son entre 22.5 °C y 28 °C, con una pluviosidad anual de entre 1500 y 3500 mm (mensual entre 125 y 290 mm).

Los datos climáticos históricos se pueden obtener de WorldClim, versión 2.1, para el período 1970-2000 (versión lanzada en enero de 2020).

1 Utilizando los datos de clima de línea base a nivel global, genere un código en R que permita construir mapas de aptitud en términos climáticos para la caña de azúcar (con base en los rangos óptimos). Grafique los mapas con una escala de colores adecuada.

##Importar Datos Espaciales en R
#install.packages("rgdal") #no voy a usar rgdal porque por mi versión de rtools ya no sirve.
#install.packages("sf")
#install.packages("terra")
#install.packages("rasterVis")
#install.packages("raster")

library(terra)

## terra 1.7.78

library(sf)

## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

library(rasterVis)

## Cargando paquete requerido: lattice

library(raster)

## Cargando paquete requerido: sp

## 
## Adjuntando el paquete: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

# Listar y cargar los archivos de precipitación en un solo objeto de SpatRaster
archivos_precip <- list.files(pattern = "wc2.1_10m_prec_.*\\.tif$", full.names = TRUE)
precipitacion <- rast(archivos_precip)

# Asignar nombres de los meses
names(precipitacion) <- month.name

# Convertir a RasterStack para usar con rasterVis
precipitacion_stack <- stack(precipitacion)

# Crear un raster binario: verde si está entre 125 y 290, gris si no
precipitacion_binaria <- calc(precipitacion_stack, fun = function(x) {
  ifelse(x >= 125 & x <= 290, 1, 0)
})

# Definir los colores: verde para 1, gris para 0
colores <- c("gray", "green")
plot(precipitacion)

# Graficar usando levelplot con la paleta de colores binaria y nombres de los meses en los paneles
levelplot(precipitacion_binaria, col.regions = colores, at = c(0, 0.5, 1),
          main = "Precipitación Mensual (125-290 en Verde)", names.attr = month.name)

library(terra)
library(raster)
library(rasterVis)

# Define el directorio de trabajo
setwd("C:/Users/doar9/OneDrive/Javeriana/analisis espacial/unidad 2/worldata")

# Listar y cargar los archivos de temperatura promedio en un solo objeto de SpatRaster
archivos_tavg <- list.files(pattern = "wc2.1_10m_tavg_.*\\.tif$", full.names = TRUE)
temperatura_promedio <- rast(archivos_tavg)

# Asignar nombres de los meses
names(temperatura_promedio) <- month.name

# Convertir a RasterStack para usar con rasterVis
temperatura_promedio_stack <- stack(temperatura_promedio)

# Crear un raster binario: verde si está en el rango deseado (ajusta el rango si es necesario)
temperatura_binaria <- calc(temperatura_promedio_stack, fun = function(x) {
  ifelse(x >= 22.5 & x <= 28, 1, 0)
})

# Definir los colores: verde para 1 (en rango), gris para 0 (fuera de rango)
colores <- c("gray", "green")
plot(temperatura_promedio)

# Graficar usando levelplot con la paleta de colores binaria y nombres de los meses en los paneles
levelplot(temperatura_binaria, col.regions = colores, at = c(0, 0.5, 1),
          main = "Temperatura Promedio Mensual (22.5 - 28 en Verde)", names.attr = month.name)

Aunque el mapa muestra lugares con temperatura promedio adecuada más al norte en meses de verano boreal y y más al sur en meses de verano austral, en los siguientes análisis se corrige para que tome el promedio apto anual y no envíe lugares más al norte y al sur.

# Convertir el SpatRaster a un data frame
precip_df <- as.data.frame(precipitacion, xy = TRUE, na.rm = TRUE)
# Convertir el SpatRaster de temperatura promedio a un data frame
temp_avg_df <- as.data.frame(temperatura_promedio, xy = TRUE, na.rm = TRUE)

# Crear la columna 'suma_precipitacion' como la suma por filas
precip_df$suma_precipitacion <- rowSums(precip_df[, 3:14])

# Crear la columna 'precip:adecuada' con la nueva condición
precip_df$`precip:adecuada` <- ifelse(precip_df$suma_precipitacion >= 1500 & precip_df$suma_precipitacion <= 3500, 1, 0)

# Crear la columna 'prom_fila' como el promedio por filas
temp_avg_df$prom_fila <- rowMeans(temp_avg_df[, 3:14])

# Crear la columna 'temp_prom_adecuada' con la condición de estar entre 22.5 y 28
temp_avg_df$temp_prom_adecuada <- ifelse(temp_avg_df$prom_fila >= 22.5 & temp_avg_df$prom_fila <= 28, 1, 0)

# Seleccionar las columnas necesarias de cada data frame
precip_filtrado <- precip_df[, c("x", "y", "precip:adecuada")]
temp_avg_filtrado <- temp_avg_df[, c("x", "y", "temp_prom_adecuada")]


# Combinar los data frames por las columnas 'x' e 'y'
mapa_cana_aptitud <- merge(precip_filtrado, temp_avg_filtrado, by = c("x", "y"))

# Mostrar el data frame combinado
head(mapa_cana_aptitud)

# Crear la columna 'apto_para_cana' con la condición de que las tres variables sean 1
mapa_cana_aptitud$apto_para_cana <- ifelse(
  mapa_cana_aptitud$`precip:adecuada` == 1 & 
  mapa_cana_aptitud$temp_prom_adecuada == 1 , 
  1, 0
)

table(mapa_cana_aptitud$`precip:adecuada`)

## 
##      0      1 
## 739671  68382

table(mapa_cana_aptitud$temp_prom_adecuada)

## 
##      0      1 
## 668780 139273

table(mapa_cana_aptitud$apto_para_cana)

## 
##      0      1 
## 758158  49895

library(ggplot2)
library(sf)

# Convertir el data frame completo en un objeto sf
mapa_cana_sf <- st_as_sf(mapa_cana_aptitud, coords = c("x", "y"), crs = 4326)

# Crear el mapa con ggplot2, coloreando según 'apto_para_cana'
ggplot() +
  geom_sf(data = mapa_cana_sf, aes(color = factor(apto_para_cana)), size = 0.5) +
  scale_color_manual(values = c("0" = "red", "1" = "green"), 
                     labels = c("No apto", "Apto"),
                     name = "Aptitud para caña") +
  labs(title = "Mapa de aptitud para caña de azúcar",
       x = "Longitud", y = "Latitud") +
  theme_minimal()

2 Identifique 2 o 3 países con áreas de alto potencial para la caña de azúcar y realice un corte para estas zonas con el shape global. Grafique los mapas con una escala de colores adecuada.

Se identifican países intertropicales cercanos al ecuador pero con alta pluviosidad, por ejemplo Brasil, Colombia o países en el sudeste asiático como Indonesia.

#install.packages("rnaturalearth")
#install.packages("rnaturalearthdata")
#install.packages("rnaturalearthhires")

library(ggplot2)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)

## 
## Adjuntando el paquete: 'rnaturalearthdata'

## The following object is masked from 'package:rnaturalearth':
## 
##     countries110

# Convertir el data frame completo en un objeto sf
mapa_cana_sf <- st_as_sf(mapa_cana_aptitud, coords = c("x", "y"), crs = 4326)

# Cargar los datos de países
paises <- ne_countries(scale = "medium", returnclass = "sf")

# Crear el mapa con ggplot2, agregando los límites y siglas de los países
ggplot() +
  geom_sf(data = paises, fill = NA, color = "black") +  # Límites de los países
  geom_sf(data = mapa_cana_sf, aes(color = factor(apto_para_cana)), size = 0.5) +
  scale_color_manual(values = c("0" = "red", "1" = "green"), 
                     labels = c("No apto", "Apto"),
                     name = "Aptitud para caña") +
  geom_sf_text(data = paises, aes(label = iso_a3), size = 2.5, color = "black") +  # Siglas de los países
  labs(title = "Mapa de aptitud para caña de azúcar",
       x = "Longitud", y = "Latitud") +
  theme_minimal()

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

3 Identificar algunos puntos (2 o 3) al azar en la región del valle del cauca (use google maps) y extraer la información de clima. Grafique las series de tiempo de temperatura y precipitación.

Se encuentran puntos hacia el centro y hacia el norte del Valle del Cauca como se muestra en el mapa

library(leaflet)
library(sf)
library(dplyr)
library(rnaturalearth)

# Obtener los límites de Colombia a nivel de departamentos
colombia <- ne_states(country = "Colombia", returnclass = "sf")

# Filtrar solo el Valle del Cauca
valle_del_cauca <- colombia %>% filter(name == "Valle del Cauca")

# Convertir el data frame de puntos aptos en un objeto sf
mapa_cana_sf <- st_as_sf(mapa_cana_aptitud, coords = c("x", "y"), crs = 4326)

# Filtrar solo los puntos aptos que están dentro del Valle del Cauca
mapa_cana_valle <- mapa_cana_sf %>%
  filter(apto_para_cana == 1) %>%
  filter(st_within(., valle_del_cauca, sparse = FALSE))

## Warning: Using one column matrices in `filter()` was deprecated in dplyr 1.1.0.
## ℹ Please use one dimensional logical vectors instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Crear el mapa con leaflet y agregar los puntos aptos en el Valle del Cauca
leaflet() %>%
  addTiles() %>%
  addPolygons(data = valle_del_cauca, color = "black", weight = 2, fill = NA) %>% # Mostrar el límite de Valle del Cauca
  addCircleMarkers(data = mapa_cana_valle,
                   color = "green", radius = 3,
                   label = ~paste("Punto Apto"),
                   popup = ~paste("Longitud:", st_coordinates(mapa_cana_valle)[,1], 
                                  "<br>Latitud:", st_coordinates(mapa_cana_valle)[,2])) %>%
  addControl("Puntos Aptos para Caña en el Valle del Cauca", position = "topright")

4 Por medio de alguna métrica de similaridad (ejemplo: distancia euclidiana) genere un código en R que permita identificar mapas de similaridad a nivel global para los sitios identificados en 3. Grafique los mapas con una escala de colores adecuada.

library(dplyr)
library(sf)
library(rnaturalearth)

# Obtener los límites de Colombia a nivel de departamentos
colombia <- ne_states(country = "Colombia", returnclass = "sf")

# Filtrar solo el Valle del Cauca
valle_del_cauca <- colombia %>% filter(name == "Valle del Cauca")

# Convertir el data frame de puntos aptos en un objeto sf
mapa_cana_sf <- st_as_sf(mapa_cana_aptitud, coords = c("x", "y"), crs = 4326)

# Filtrar solo los puntos aptos que están dentro del Valle del Cauca
mapa_cana_valle <- mapa_cana_sf %>%
  filter(apto_para_cana == 1) %>%
  filter(st_within(., valle_del_cauca, sparse = FALSE))

# Extraer las coordenadas de los puntos aptos en el Valle del Cauca y crear una tabla
coords_valle <- st_coordinates(mapa_cana_valle)
coords_valle_df <- as.data.frame(coords_valle) %>%
  rename(longitud = X, latitud = Y) %>%
  mutate(punto = paste0("punto ", row_number()))

# Mostrar la tabla de coordenadas
coords_valle_df

library(dplyr)
library(ggplot2)
library(tidyr)

## 
## Adjuntando el paquete: 'tidyr'

## The following object is masked from 'package:raster':
## 
##     extract

## The following object is masked from 'package:terra':
## 
##     extract

# Filtrar precip_df y temp_avg_df para seleccionar solo los puntos en coords_valle_df
precip_valle <- precip_df %>%
  semi_join(coords_valle_df, by = c("x" = "longitud", "y" = "latitud"))

temp_valle <- temp_avg_df %>%
  semi_join(coords_valle_df, by = c("x" = "longitud", "y" = "latitud"))

# Convertir precip_valle a formato largo
precip_long <- precip_valle %>%
  pivot_longer(cols = January:December, names_to = "Mes", values_to = "Precipitación") %>%
  left_join(coords_valle_df, by = c("x" = "longitud", "y" = "latitud"))

# Convertir temp_valle a formato largo
temp_long <- temp_valle %>%
  pivot_longer(cols = January:December, names_to = "Mes", values_to = "Temperatura") %>%
  left_join(coords_valle_df, by = c("x" = "longitud", "y" = "latitud"))

# Crear un gráfico de líneas para la precipitación
ggplot(precip_long, aes(x = Mes, y = Precipitación, group = punto, color = punto)) +
  geom_line() +
  labs(title = "Serie de Tiempo de Precipitación - Puntos Aptos en el Valle del Cauca",
       x = "Mes", y = "Precipitación (mm)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Crear un gráfico de líneas para la temperatura
ggplot(temp_long, aes(x = Mes, y = Temperatura, group = punto, color = punto)) +
  geom_line() +
  labs(title = "Serie de Tiempo de Temperatura - Puntos Aptos en el Valle del Cauca",
       x = "Mes", y = "Temperatura (°C)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

library(leaflet)
library(dplyr)
library(sf)

# Añadir la columna `punto` a `mapa_cana_valle` directamente
mapa_cana_valle$punto <- coords_valle_df$punto

# Crear el mapa con leaflet y agregar los puntos aptos en el Valle del Cauca con etiquetas
leaflet() %>%
  addTiles() %>%
  addPolygons(data = valle_del_cauca, color = "black", weight = 2, fill = NA) %>% # Mostrar el límite de Valle del Cauca
  addCircleMarkers(data = mapa_cana_valle,
                   color = "green", radius = 3,
                   label = ~punto,  # Etiqueta con el número del punto
                   popup = ~paste("Longitud:", st_coordinates(mapa_cana_valle)[,1], 
                                  "<br>Latitud:", st_coordinates(mapa_cana_valle)[,2])) %>%
  addControl("Puntos Aptos para Caña en el Valle del Cauca", position = "topright")

library(dplyr)
library(sf)
library(ggplot2)

# Paso 1: Calcular la suma de precipitación y el promedio de temperatura para los puntos de referencia en el Valle del Cauca

# Filtrar precip_valle y temp_valle para calcular la suma y el promedio de los 8 puntos aptos en el Valle del Cauca
precip_valle <- precip_df %>%
  semi_join(coords_valle_df, by = c("x" = "longitud", "y" = "latitud")) %>%
  mutate(suma_precipitacion = rowSums(.[, 3:14]))  # Usando las posiciones de las columnas de los meses

temp_valle <- temp_avg_df %>%
  semi_join(coords_valle_df, by = c("x" = "longitud", "y" = "latitud")) %>%
  mutate(prom_fila = rowMeans(.[, 3:14]))  # Usando las posiciones de las columnas de los meses

# Crear un dataframe con las medidas de referencia para comparación
valle_referencia <- data.frame(
  x = precip_valle$x,
  y = precip_valle$y,
  suma_precipitacion = precip_valle$suma_precipitacion,
  prom_fila = temp_valle$prom_fila
)

# Paso 2: Calcular la distancia euclidiana de todos los puntos a los puntos de referencia y marcar los similares

# Unir precip_df y temp_avg_df con sus columnas de suma y promedio
precip_df <- precip_df %>% mutate(suma_precipitacion = rowSums(.[, 3:14]))
temp_avg_df <- temp_avg_df %>% mutate(prom_fila = rowMeans(.[, 3:14]))

# Unir precip_df y temp_avg_df en un solo dataframe para calcular distancias
puntos <- precip_df %>%
  dplyr::select(x, y, suma_precipitacion) %>%
  left_join(temp_avg_df %>% dplyr::select(x, y, prom_fila), by = c("x", "y"))

# Calcular la distancia euclidiana para cada punto en comparación con los puntos de referencia
puntos <- puntos %>%
  rowwise() %>%
  mutate(
    distancia_minima = min(sqrt((suma_precipitacion - valle_referencia$suma_precipitacion)^2 +
                                (prom_fila - valle_referencia$prom_fila)^2))
  ) %>%
  ungroup()

# Definir un umbral de similitud (puedes ajustar este valor según los datos)
umbral_distancia <- quantile(puntos$distancia_minima, 0.01) # Seleccionar el 1% más cercano

# Marcar los puntos similares
puntos <- puntos %>%
  mutate(parecidos_a_valle_del_cauca = ifelse(distancia_minima <= umbral_distancia, 1, 0))

# Paso 3: Convertir a un objeto sf y graficar en ggplot2

# Convertir puntos a sf
puntos_sf <- st_as_sf(puntos, coords = c("x", "y"), crs = 4326)

# Crear el mapa en ggplot2, coloreando los puntos según su similitud con los puntos del Valle del Cauca
ggplot() +
  geom_sf(data = puntos_sf, aes(color = factor(parecidos_a_valle_del_cauca)), size = 0.5) +
  scale_color_manual(values = c("0" = "gray", "1" = "blue"), 
                     labels = c("No similar", "Similar"),
                     name = "Similitud con Valle del Cauca") +
  labs(title = "Mapa de Similitud de Puntos con Valle del Cauca",
       x = "Longitud", y = "Latitud") +
  theme_minimal()

5 Compare los mapas generados por ambas aproximaciones y concluya.

Se realizó un análisis para identificar puntos climáticamente similares a aquellos en el Valle del Cauca, utilizando la distancia euclidiana basada en la precipitación acumulada y la temperatura promedio anual. Para ello, se estableció un umbral del 1% (percentil 0.01) para seleccionar puntos con características de precipitación y temperatura similares a los puntos aptos de caña de azúcar en el Valle del Cauca. Los mapas y gráficos obtenidos muestran tanto la distribución geográfica de los puntos similares en el mundo como las series de tiempo de precipitación y temperatura para los puntos aptos en el Valle del Cauca. Esta metodología permite encontrar áreas con condiciones adecuadas para el cultivo de caña de azúcar en otras regiones. Las técnicas de visualización en los mapas proporcionan una representación geográfica precisa, facilitando la identificación de ubicaciones con características similares en distintas partes del mundo.

Fin del documento

Análisis de información Geográfica y Espacial

David Aguilar.