Planteamiento del Problema

Se desea conocer en qué zonas del mundo se puede establecer la producción comercial apícola para obtener los mejores rendimientos teniendo en cuenta los rangos de aptitud climática

Revisando información de la UPRA se conocen los rangos optimos de condición climática:

Se desea conocer qué zonas en colombia presentan estas condiciones y conocer otros sitios en el mundo que podrían tener condiciones similares.

Ver https://sipra.upra.gov.co/

Preparación de los datos

Obtención de los datos

Se hace uso de los datos de clima de línea base a nivel global del siguiente enlace, los cuales corresponden al periodo 1970 - 2000, cuya última actualización fue en enero 2020.

https://www.worldclim.org/data/worldclim21.html

Se realiza la carga de datos de precipitación y temperatura media anual en una resolución espacial de 5 minutos.

#Carga de librerías
library(leaflet)           # Para la creación de mapas interactivos.
library(sf)                # Para trabajar con datos geoespaciales.
library(mapview)           # Para la visualización interactiva de mapas.
library(terra)             # Para el manejo de datos raster y análisis espacial.
library(leafsync)   # Para la sincronización de múltiples mapas interactivos.
library(raster)     # Para trabajar con datos espaciales como imagenes
library(rasterVis)
library(sp)         # Para trabajar con datos espaciales en formatos vectoriales
library(dplyr)             #Para la manipulación y transformación de datos.
library(rnaturalearth)     #Para obtener mapas a distintas escalas. 
library(rnaturalearthdata)
library(exactextractr)
library(geodata)
#Cargue de datos de precipitación
Prec=list.files("D:/MaestriaCienciaDatos/S2_AnalisisInfoGeografica/M1U2_Exploracion/CasoAccion/DatosClima/Prec/", pattern = ".tif", full.names = TRUE)
precipitacion=raster::stack(Prec)

Visualización mensual de la precipitación media (mm)

require(rasterVis)
names(precipitacion)=month.name
levelplot(precipitacion,
          par.settings = BuRdTheme,
          at = seq(0, 1500, by = 50))

# Calcular la media mensual por capa
monthly_mean <- cellStats(precipitacion, stat = 'mean', na.rm = TRUE)

# Graficar
barplot(monthly_mean,
        names.arg = month.name,
        las = 2,
        col = "skyblue",
        main = "Precipitación promedio mensual",
        xlab = "Meses", ylab = "Precipitación promedio (mm)")

Visualización mensual de la temperatura media (C°)

#Cargue de datos de temperatura
require(rasterVis)
Temp=list.files("D:/MaestriaCienciaDatos/S2_AnalisisInfoGeografica/M1U2_Exploracion/CasoAccion/DatosClima/Temp/", full.names = TRUE)
temperatura=stack(Temp)
names(temperatura)=month.name
levelplot(temperatura,
          par.settings = BuRdTheme,
          at = seq(-40, 50, by = 5))

# Calcular la media mensual por capa
monthly_meant <- cellStats(temperatura, stat = 'mean', na.rm = TRUE)

# Graficar
barplot(monthly_meant,
        names.arg = month.name,
        las = 2,
        col = "skyblue",
        main = "Temperatura promedio mensual",
        xlab = "Meses", ylab = "Temperatura promedio (°C)")

Se destaca que la tendencia en la mayor cantidad de países del mundo es mantener una precipítación media baja, cuyos mayores valores se observan en los países ubicados a lo largo de la línea del Ecuador. Los meses con mayor precipitación son Julio y Agosto influenciados particularmente por regiones tropicales como Suramérica, el Sureste de Asia y algunos sectores de África.

En el caso de la temperatura, se destaca una gran variabilidad climática en dónde predominan los valores extremos altos sobre los bajos. Lo cual impacta en gran medida la temperatura promedio mensual, generando valores por debajo de 0°C, se destacan los meses junio, julio y agosto como aquellos con mayor temperatura.

Generación del Mapa de Aptitud Climática

Selección de rangos óptimos

Se filtraron las áreas teniendo en cuenta los rangos óptimos establecidos tanto para la precipitación como para la temperatura en la producción comercial apícola.

Zonas con precipitación óptima

#Precipitación media anual: entre 500 y 1800 mm/año
#Precipitación mensual optima: entre 42 y 150 mm
prec_apicola <- precipitacion>42 & precipitacion<150
plot(prec_apicola)

Zonas con temperatura óptima

#Temperatura media anual: mayor a 12 grados centígrados
temp_apicola <- temperatura>12 
plot(temp_apicola)

Identificación de países con alto potencial para establecer la producción comercial apícola

prec_bin <- prec_apicola*1
# -------------------------
# Se entiende la celda apta si algún mes del año está en el rango (binario 0/)
prec_agg <- calc(prec_bin, fun = function(x) { as.integer(any(x == 1, na.rm = TRUE)) })
 
# -------------------------
# Reproyectar a proyección métrica / equal-area para cálculo de área real (si no está ya en metros)
# Verificar si las unidades del raster son metros (crs tiene +units=m) o EPSG proyectado
crs_r <- crs(prec_agg)

# Si el CRS no es métrico, reproyecta; usamos Mollweide (equal area) como ejemplo
is_metric_crs <- grepl("units=m|+units=m|+proj=aea|+proj=moll|EPSG:\\d+", as.character(crs_r), fixed = FALSE)
# (simplista: si no estás seguro, reproyecta de todas formas)
prec_m <- projectRaster(prec_agg, crs = "+proj=moll +units=m", method = "ngb")

# -------------------------
# Preparar países en mismo CRS
paises <- ne_countries(scale = "medium", returnclass = "sf")
paises_m <- st_transform(paises, crs = st_crs(prec_m))

# -------------------------
# Extraer solo países que intersectan la extensión del raster (evita mismatches)
ext_poly <- as(extent(prec_m), "SpatialPolygons")
crs(ext_poly) <- crs(prec_m)
ext_sf <- st_as_sf(ext_poly)
intersect_idx <- st_intersects(paises_m, ext_sf, sparse = FALSE)[,1]
paises_to_extract <- paises_m[intersect_idx, ]

# -------------------------
# Usar exact_extract sobre el raster AGREGADO (ya single-layer)
# Esto devuelve para cada país la SUMA de valores del raster dentro del polígono
ae <- exact_extract(prec_m, paises_to_extract, 'sum', progress = FALSE)

# Convertir a vector (un número por país)
area_count_subset <- sapply(ae, function(x) {
  # 'x' puede ser un número, vector, df ... aseguramos sumar
  if (is.data.frame(x)) return(sum(as.matrix(x), na.rm = TRUE))
  if (is.numeric(x)) return(sum(x, na.rm = TRUE))
  if (is.list(x)) { y <- unlist(x); if (is.numeric(y)) return(sum(y, na.rm = TRUE)) }
  return(0)
})

# Rellenar vector para todos los países (0 si no intersectaron)
area_count_all <- numeric(nrow(paises_m))
area_count_all[which(intersect_idx)] <- area_count_subset

# -------------------------
# Convertir conteo a área en km² 
#      * area de celdas aptas en total de meses = area_count_all * area_celda (esto da "meses·area")
area_celda_m2 <- prod(res(prec_m))   # m^2 por celda (res en metros)

area_apta_km2 <- (area_count_all * area_celda_m2) / 1e6

# -------------------------
# Calcular área total país en km2 y proporción
area_total_km2 <- as.numeric(st_area(paises_m)) / 1e6
# Prevenir div/0
area_total_km2[area_total_km2 == 0] <- NA

prop_apta <- area_apta_km2 / area_total_km2

# -------------------------
# Adjuntar y mostrar resultados
paises_m$area_total_km2 <- area_total_km2
paises_m$area_apta_km2  <- area_apta_km2
paises_m$prop_apta      <- prop_apta

resultados <- paises_m %>%
  st_drop_geometry() %>%
  select(admin, area_total_km2, area_apta_km2, prop_apta) %>%
  arrange(desc(prop_apta))

knitr::kable(head(resultados, 20), caption = "Top 20 países con precipitación optima")
Top 20 países con precipitación optima
admin area_total_km2 area_apta_km2 prop_apta
Saint Barthelemy 24.57829 24.57829 1
Andorra 442.79947 442.79949 1
Liechtenstein 138.27521 138.27522 1
Malawi 120062.21902 120062.22410 1
North Macedonia 25420.81768 25420.81869 1
Rwanda 25442.41080 25442.41180 1
Macao S.A.R 22.79087 22.79088 1
Republic of Serbia 77663.53604 77663.53862 1
Paraguay 401803.25237 401803.26555 1
Lesotho 30094.79370 30094.79464 1
eSwatini 17030.08532 17030.08580 1
Monaco 12.44278 12.44278 1
Slovakia 48447.50660 48447.50784 1
Ivory Coast 323018.43686 323018.44409 1
Saint Vincent and the Grenadines 311.78799 311.78799 1
South Sudan 625360.44540 625360.45469 1
Montserrat 79.96314 79.96314 1
Zambia 757286.69822 757286.70738 1
Moldova 33191.15272 33191.15312 1
Bosnia and Herzegovina 51800.29019 51800.29077 1 
temp_bin <- temp_apicola*1
# -------------------------
# Se entiende la celda apta si algún mes del año está en el rango (binario 0/)
temp_agg <- calc(temp_bin, fun = function(x) { as.integer(any(x == 1, na.rm = TRUE)) })
 
# -------------------------
# Reproyectar a proyección métrica / equal-area para cálculo de área real (si no está ya en metros)
# Verificar si las unidades del raster son metros (crs tiene +units=m) o EPSG proyectado
crs_r <- crs(temp_agg)

# Si el CRS no es métrico, reproyecta; usamos Mollweide (equal area) como ejemplo
is_metric_crs <- grepl("units=m|+units=m|+proj=aea|+proj=moll|EPSG:\\d+", as.character(crs_r), fixed = FALSE)
# (simplista: si no estás seguro, reproyecta de todas formas)
temp_m <- projectRaster(temp_agg, crs = "+proj=moll +units=m", method = "ngb")

# -------------------------
# Preparar países en mismo CRS
paises <- ne_countries(scale = "medium", returnclass = "sf")
paises_m <- st_transform(paises, crs = st_crs(temp_m))

# -------------------------
# Extraer solo países que intersectan la extensión del raster (evita mismatches)
ext_poly <- as(extent(temp_m), "SpatialPolygons")
crs(ext_poly) <- crs(temp_m)
ext_sf <- st_as_sf(ext_poly)
intersect_idx <- st_intersects(paises_m, ext_sf, sparse = FALSE)[,1]
paises_to_extract <- paises_m[intersect_idx, ]

# -------------------------
# Usar exact_extract sobre el raster AGREGADO (ya single-layer)
# Esto devuelve para cada país la SUMA de valores del raster dentro del polígono
ae <- exact_extract(temp_m, paises_to_extract, 'sum', progress = FALSE)

# Convertir a vector (un número por país)
area_count_subset <- sapply(ae, function(x) {
  # 'x' puede ser un número, vector, df ... aseguramos sumar
  if (is.data.frame(x)) return(sum(as.matrix(x), na.rm = TRUE))
  if (is.numeric(x)) return(sum(x, na.rm = TRUE))
  if (is.list(x)) { y <- unlist(x); if (is.numeric(y)) return(sum(y, na.rm = TRUE)) }
  return(0)
})

# Rellenar vector para todos los países (0 si no intersectaron)
area_count_all <- numeric(nrow(paises_m))
area_count_all[which(intersect_idx)] <- area_count_subset

# -------------------------
# Convertir conteo a área en km² 
#      * area de celdas aptas en total de meses = area_count_all * area_celda (esto da "meses·area")
area_celda_m2 <- prod(res(temp_m))   # m^2 por celda (res en metros)

area_apta_km2 <- (area_count_all * area_celda_m2) / 1e6

# -------------------------
# Calcular área total país en km2 y proporción
area_total_km2 <- as.numeric(st_area(paises_m)) / 1e6
# Prevenir div/0
area_total_km2[area_total_km2 == 0] <- NA

prop_apta <- area_apta_km2 / area_total_km2

# -------------------------
# Adjuntar y mostrar resultados
paises_m$area_total_km2 <- area_total_km2
paises_m$area_apta_km2  <- area_apta_km2
paises_m$prop_apta      <- prop_apta

resultados <- paises_m %>%
  st_drop_geometry() %>%
  select(admin, area_total_km2, area_apta_km2, prop_apta) %>%
  arrange(desc(prop_apta))

knitr::kable(head(resultados, 20), caption = "Top 20 países con temperatura optima")
Top 20 países con temperatura optima
admin area_total_km2 area_apta_km2 prop_apta
Saint Barthelemy 2.457829e+01 2.457829e+01 1
Malawi 1.200622e+05 1.200622e+05 1
Macao S.A.R 2.279087e+01 2.279088e+01 1
Republic of Serbia 7.766354e+04 7.766354e+04 1
Paraguay 4.018033e+05 4.018033e+05 1
Republic of the Congo 3.476507e+05 3.476507e+05 1
eSwatini 1.703009e+04 1.703009e+04 1
Chad 1.273841e+06 1.273841e+06 1
Monaco 1.244278e+01 1.244278e+01 1
Ivory Coast 3.230184e+05 3.230184e+05 1
Iraq 4.384353e+05 4.384353e+05 1
Saint Vincent and the Grenadines 3.117880e+02 3.117880e+02 1
South Sudan 6.253604e+05 6.253605e+05 1
Montserrat 7.996314e+01 7.996314e+01 1
Western Sahara 9.099750e+04 9.099750e+04 1
Zambia 7.572867e+05 7.572867e+05 1
Moldova 3.319115e+04 3.319115e+04 1
Bosnia and Herzegovina 5.180029e+04 5.180029e+04 1
San Marino 6.841198e+01 6.841198e+01 1
Turkmenistan 4.716821e+05 4.716821e+05 1

Revisando el top 20 de países con precipitación y temperatura optima se encuentra que el Sahara Occidental y Chad, territorios ubicados en África, cuentan con alto potencial para la la producción comercial apícola.

# --- 1. Obtener el polígono de Sahara Occidental ---
map_units <- ne_download(
  scale = "medium",
  type = "admin_0_map_units",
  category = "cultural",
  returnclass = "sf"
)
ws_sf <- map_units %>%
  filter(GEOUNIT == "Western Sahara")
ws_sf <- st_transform(ws_sf, crs(precipitacion))
ws_sp <- as(ws_sf, "Spatial")  # convertir a formato Spatial (para raster)

# --- 2. Filtrar zonas con precipitación óptima ---
# Suponiendo que 'precipitacion' es tu RasterStack o RasterBrick mensual
prec_apicola <- clamp(precipitacion, lower = 42, upper = 150, useValues = TRUE)



# --- 3. Recortar y enmascarar el raster a Western Sahara ---
prec_apicola_ws <- crop(prec_apicola, ws_sp)
prec_apicola_ws <- mask(prec_apicola_ws, ws_sp)

# --- 4. Graficar el resultado ---
levelplot(prec_apicola_ws,
          par.settings = BuRdTheme,
          at = seq(42, 150, by = 2),
          main = "Zonas con precipitacion optima para apicultura (Sahara Occidental)",
          margin = FALSE,
          sp.layout = list(list("sp.polygons", ws_sp, col = "black", lwd = 1.5)))

#ws_sf <- st_transform(ws_sf, crs(temperatura))
ws_sp <- as(ws_sf, "Spatial")  # convertir a formato Spatial (para raster)

# --- 2. Filtrar zonas con precipitación óptima ---
# Suponiendo que 'precipitacion' es tu RasterStack o RasterBrick mensual
temp_apicola <- clamp(temperatura, lower = 12, useValues = TRUE)

# --- 3. Recortar y enmascarar el raster a Western Sahara ---
temp_apicola_ws <- crop(temp_apicola, ws_sp)
temp_apicola_ws <- mask(temp_apicola_ws, ws_sp)

# --- 4. Graficar el resultado ---
levelplot(temp_apicola_ws,
          par.settings = BuRdTheme,
          at = seq(12, 50, by = 2),
          main = "Zonas con temperatura optima para apicultura (Sahara Occidental)",
          margin = FALSE,
          sp.layout = list(list("sp.polygons", ws_sp, col = "black", lwd = 1.5)))

# --- 1. Obtener el polígono de Chad ---
Chad_sf <- ne_countries(scale = "medium", country = "Chad", returnclass = "sf")
#Chad_sf <- st_transform(Chad_sf, crs = st_crs(prec_apicola_Chad))
Chad_sp <- as(Chad_sf, "Spatial")

# --- 2. Filtrar zonas con precipitación óptima ---
# Suponiendo que 'precipitacion' es tu RasterStack o RasterBrick mensual
prec_apicola <- mask(precipitacion, precipitacion > 42 & precipitacion < 150)


# --- 3. Recortar y enmascarar el raster a Chad ---
prec_apicola_Chad <- crop(prec_apicola, Chad_sp)
prec_apicola_Chad <- mask(prec_apicola_Chad, Chad_sp)

# --- 4. Graficar el resultado ---
levelplot(prec_apicola_Chad,
          par.settings = BuRdTheme,
          at = seq(42, 150, by = 2),
          main = "Zonas con precipitacion optima para apicultura (Chad)",
          margin = FALSE,
          sp.layout = list(
            list("sp.polygons", Chad_sp, col = "black", lwd = 1.5))
          )

# --- Filtrar zonas con precipitación óptima ---
# Suponiendo que 'precipitacion' es tu RasterStack o RasterBrick mensual
temp_apicola <- mask(temperatura, temperatura > 12)


# --- Recortar y enmascarar el raster a Chad ---
temp_apicola_Chad <- crop(temp_apicola, Chad_sp)
temp_apicola_Chad <- mask(temp_apicola_Chad, Chad_sp)

# --- Graficar el resultado ---
levelplot(temp_apicola_Chad,
          par.settings = BuRdTheme,
          at = seq(12, 50, by = 2),
          main = "Zonas con temperatura optima para apicultura (Chad)",
          margin = FALSE,
          sp.layout = list(
            list("sp.polygons", Chad_sp, col = "black", lwd = 1.5))
          )

# --- 1. Obtener el polígono de Colombia ---
colombia_sf <- ne_countries(scale = "medium", country = "Colombia", returnclass = "sf")
colombia_sp <- as(colombia_sf, "Spatial")  # convertir a formato Spatial (para raster)

# --- 2. Filtrar zonas con precipitación óptima ---
# Suponiendo que 'precipitacion' es tu RasterStack o RasterBrick mensual
prec_apicola <- mask(precipitacion, precipitacion > 42 & precipitacion < 150)


# --- 3. Recortar y enmascarar el raster a Colombia ---
prec_apicola_col <- crop(prec_apicola, colombia_sp)
prec_apicola_col <- mask(prec_apicola_col, colombia_sp)

# --- 4. Graficar el resultado ---
levelplot(prec_apicola_col,
          par.settings = BuRdTheme,
          at = seq(42, 150, by = 2),
          main = "Zonas con precipitacion optima para apicultura (Colombia)")

# --- Filtrar zonas con temperatura óptima ---
temp_apicola <- mask(temperatura, temperatura > 12)


# --- Recortar y enmascarar el raster a Colombia ---
temp_apicola_col <- crop(temp_apicola, colombia_sp)
temp_apicola_col <- mask(temp_apicola_col, colombia_sp)

# --- Graficar el resultado ---
levelplot(temp_apicola_col,
          par.settings = BuRdTheme,
          at = seq(12, 40, by = 2),
          main = "Zonas con temperatura optima para apicultura (Colombia)")

Identificación de puntos al azar en Colombia

Priorización de áreas

Para unificar el ejercicio se promedian los resultados de temperatura y precipitación en condiciones optimas para la producción comercial apícola. A partir de estos resultados se seleccionan aquellas zonas por encima del 80% con el fin de seleccionar las mejores áreas.

# Promedio anual
prec_optima_apicultura <- sum(prec_apicola) / 12
temp_optima_apicultura <- sum(temp_apicola) / 12

# Recortar a Colombia
optimo_colombia_prec <- crop(prec_optima_apicultura, colombia_sp)
optimo_colombia_temp <- crop(temp_optima_apicultura, colombia_sp)

# Calcular umbrales del 80% del promedio
umbral_prec <- mean(values(optimo_colombia_prec), na.rm = TRUE) * 0.8
umbral_temp <- mean(values(optimo_colombia_temp), na.rm = TRUE) * 0.8

# Asignar NA a celdas no aptas (menores al 80% del promedio)
optimo_colombia_prec[optimo_colombia_prec < umbral_prec] <- NA
optimo_colombia_temp[optimo_colombia_temp < umbral_temp] <- NA

# Enmascarar con el shape de Colombia
optimo_col_prec <- mask(optimo_colombia_prec, colombia_sp)
optimo_col_temp <- mask(optimo_colombia_temp, colombia_sp)

# Crear mapas interactivos
m1 <- mapview(optimo_col_prec, maxpixels = 2332800)
m2 <- mapview(optimo_col_temp, maxpixels = 2332800)

# Sincronizar
m <- leafsync::sync(m1, m2)
m

Series temporales de puntos al azar

Inicialmente se establecen las zonas en las que se cumplen las condiciones de precipitación y temperatura, posteriormente se eligen 3 puntos en esas áreas y se extraen las respectivas series temporales.

# --- Combinar zonas óptimas
# Convertir a la misma resolución y extensión
optimo_colombia_prec <- resample(optimo_colombia_prec, optimo_colombia_temp)

# Intersección (solo donde ambas son válidas)
zona_optima <- !is.na(optimo_colombia_prec) & !is.na(optimo_colombia_temp)

# Aplicar máscara a Colombia
zona_optima <- mask(zona_optima, colombia_sp)
# --- Selección de 3 puntos aleatorios dentro de la zona

# Convertir zona_optima (o cualquier RasterLayer) a SpatRaster
zona_optima_spat <- terra::rast(zona_optima)

set.seed(456)
puntos_random <- terra::spatSample(
  zona_optima_spat,
  size = 3,
  method = "random",
  as.points = TRUE,
  na.rm = TRUE  # <- ignora los NA
)

plot(zona_optima_spat)
points(puntos_random, pch = 19, col = "red")

# --- Extraer las series temporales para cada punto

# Convertir puntos a Spatial (clase de raster)
puntos_random_sp <- as(puntos_random, "Spatial")

# Extraer valores
valores_prec <- raster::extract(precipitacion, puntos_random_sp)
valores_temp <- raster::extract(temperatura, puntos_random_sp)

# Crear vector con los nombres de los meses (suponiendo que son 12 capas mensuales)
meses <- 1:nlayers(precipitacion)

# Convertir a data.frame para graficar fácilmente
df_prec <- as.data.frame(t(valores_prec))
df_temp <- as.data.frame(t(valores_temp))

colnames(df_prec) <- paste0("Punto_", 1:ncol(df_prec))
colnames(df_temp) <- paste0("Punto_", 1:ncol(df_temp))

df_prec$Mes <- 1:nrow(df_prec)
df_temp$Mes <- 1:nrow(df_temp)

# ─── Graficar precipitación ───────────────────────────────────────────────
matplot(df_prec$Mes, df_prec[, 1:3], type = "l", lwd = 2, lty = 1,
        col = c("blue", "red", "green"),
        xlab = "Mes", ylab = "Precipitacion (mm)",
        main = "Serie de tiempo de precipitacion")

legend("topright", legend = paste("Punto", 1:3),
       col = c("blue", "red", "green"), lwd = 2, bty = "n")

# ─── Graficar temperatura ────────────────────────────────────────────────
matplot(df_temp$Mes, df_temp[, 1:3], type = "l", lwd = 2, lty = 1,
        col = c("orange", "purple", "brown"),
        xlab = "Mes", ylab = "Temperatura (C)",
        main = "Serie de tiempo de temperatura")

legend("topright", legend = paste("Punto", 1:3),
       col = c("orange", "purple", "brown"), lwd = 2, bty = "n")

Mapas de similaridad a nivel global

A partir de la distancia euclidiana se generan mapas de similaridad a nivel global para los sitios anteriormente identificados. Valores bajos (azules) representan aquellos con climas más similares al punto. Valores altos (rojos) representan aquellos con climas más diferentes del punto.

library(RColorBrewer)

# ─── 1. Convertir tus capas a RasterBrick si aún no lo están ─────────────
prec_brick <- brick(precipitacion)
temp_brick <- brick(temperatura)

# ─── 2. Extraer series temporales promedio en los 3 puntos ──────────────
valores_prec <- raster::extract(prec_brick, puntos_random_sp)
valores_temp <- raster::extract(temp_brick, puntos_random_sp)

# Calcular el promedio mensual por punto (manteniendo estructura mensual)
valores_clima <- lapply(1:3, function(i) {
  c(valores_prec[i, ], valores_temp[i, ])
})
names(valores_clima) <- paste0("Punto_", 1:3)

# ─── 3. Combinar los raster globales de precipitación y temperatura ─────
clima_stack <- stack(prec_brick, temp_brick)

# ─── 4. Calcular distancia euclidiana de cada píxel a cada punto ────────
distancias <- lapply(1:3, function(i) {
  punto_ref <- valores_clima[[i]]
  raster::calc(clima_stack, fun = function(x) {
    if (any(is.na(x))) return(NA)
    sqrt(sum((x - punto_ref)^2, na.rm = TRUE))
  })
})

# ─── 5. Asignar nombres y combinar en un stack ──────────────────────────
names(distancias) <- paste0("Distancia_Punto_", 1:3)
dist_stack <- stack(distancias)

# ─── 6. Visualizar mapas con escala de colores adecuada ─────────────────
paleta <- rev(brewer.pal(11, "RdYlBu"))

par(mfrow = c(1, 3))
for (i in 1:3) {
  plot(dist_stack[[i]],
       main = paste(names(dist_stack)[i]),
       col = paleta,
       zlim = c(0, quantile(dist_stack[[i]][], 0.95, na.rm = TRUE))) # recorte a 95%
}

Conclusiones