Diagnóstico espacio–temporal 2024 de variables meteorológicas de la REDMET (México)

Author

Jairo Andrés Angel, Sara Nathalia Reina

1. Contexto y fuente de datos

Este trabajo utiliza datos horarios del año 2024 provenientes de la BASE DE DATOS DE LA RED DE METEOROLOGÍA Y RADIACIÓN SOLAR (REDMET) de la Ciudad de México.
La consulta y descarga de información se realiza desde el portal oficial de SEDEMA CDMX:
https://www.aire.cdmx.gob.mx/default.php?opc=‘aKBi’

La REDMET integra estaciones automáticas con registros puntuales en el espacio y horarios en el tiempo, adecuados para análisis geoestadísticos y diagnósticos operativos.

2. Objetivo general

Entregar a la empresa un diagnóstico espacio–temporal con resolución horaria para 2024 de las variables meteorológicas clave medidas por la REDMET, con el fin de soportar decisiones operativas (planeación de actividades en campo, evaluación de riesgos, control de calidad de datos y preparación de modelos predictivos).

3. Variables y unidades

  • Temperatura del aire (TMP) — grados Celsius (°C)
  • Humedad relativa (RH) — porcentaje (%)
  • Dirección del viento (WDR) — grados azimut (°)
  • Velocidad del viento (WSP) — metros por segundo (m/s)

Resolución temporal: series horarias con sellos de tiempo tipo dd/mm/aaaa hh:mm (p. ej., 01/01/2024 01:00, 01/01/2024 02:00, …).
Cobertura espacial: estaciones puntuales georreferenciadas en el territorio de la CDMX y áreas aledañas.

4. Planteamiento del problema

Para una operación eficiente, la empresa necesita cuantificar cómo varían en el espacio y a lo largo del tiempo las condiciones atmosféricas que inciden directamente en la seguridad, productividad y desempeño de sus procesos. En particular, se requiere:

  1. Integrar y depurar las series horarias 2024 (detección de faltantes/atípicos y estandarización de unidades).
  2. Localizar y mapear las estaciones, superponiéndolas al área de interés para identificar gradientes y zonas con condiciones atípicas.
  3. Caracterizar patrones temporales (ciclo diurno/semanal) y patrones espaciales básicos (mapas de valores y semivariograma exploratorio).
  4. Evaluar estacionariedad en media para habilitar modelos geoestadísticos sobre un proceso residual de media cero en etapas posteriores.

El enfoque parte de cortes horarios representativos para validar la metodología y luego escala a todo 2024 con rutinas reproducibles (agregaciones diarias/mensuales y tableros/figuras de soporte).

Mapa de estaciones

knitr::include_graphics("C:/Users/sarar/Documents/Maestría Estadística/Espacial/Estadística Espacial/Rplot01.png")
Figure 1: Mapa de estaciones 

library(dplyr)

Adjuntando el paquete: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
library(readr)
library(lubridate)

Adjuntando el paquete: 'lubridate'
The following objects are masked from 'package:base':

    date, intersect, setdiff, union
library(sf)
Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library(janitor)

Adjuntando el paquete: 'janitor'
The following objects are masked from 'package:stats':

    chisq.test, fisher.test
# --- Parámetros ---
TZ <- "America/Mexico_City"
dia_obj    <- as.Date("2024-11-25")
hora_texto <- "12:00"

# --- Rutas ---
path_datos      <- "C:/Users/sarar/Documents/Maestría Estadística/Datos/meteorologia_2024.csv"
path_estaciones <- "C:/Users/sarar/Documents/Maestría Estadística/Espacial/Estadística Espacial/estaciones_planas.csv"

# --- Lectura ---
dat_raw <- read_delim(path_datos, delim = ";", locale = locale(encoding = "ISO-8859-1"),show_col_types = FALSE) |> clean_names()

est_raw <- read_csv(path_estaciones, locale = locale(encoding = "ISO-8859-1"), show_col_types = FALSE) |> clean_names()

# --- Tipos + fecha/hora robusta ---
dat <- dat_raw |>
  mutate(
    id_station   = as.character(id_station),
    id_parameter = as.character(id_parameter),
    datetime     = parse_date_time(as.character(date), orders = c("dmy HM","dmy"), tz = TZ),
    date         = as.Date(datetime)
  )

est <- est_raw |>
  mutate(
    cve_estac = as.character(cve_estac),
    epsg      = as.integer(epsg)
  )

# --- Merge ---
datos_merged <- inner_join(dat, est, by = c("id_station" = "cve_estac"))

# --- EPSG y sf ---
epsg_crs <- datos_merged$epsg |> na.omit() |> unique()
if (length(epsg_crs) == 0) stop("No hay EPSG en estaciones.")
if (length(epsg_crs) > 1)  warning("Múltiples EPSG: ", paste(epsg_crs, collapse = ", "), " (uso el primero).")
epsg_crs <- epsg_crs[1]

stopifnot(all(c("easting","northing") %in% names(datos_merged)))

datos_sf <- st_as_sf(datos_merged, coords = c("easting","northing"), crs = epsg_crs, remove = FALSE)

# --- Instante objetivo (misma hora para todos) ---
instante_obj <- ymd_hm(paste(dia_obj, hora_texto), tz = TZ)

datos_instante <- datos_sf |>
  filter(date == dia_obj, floor_date(datetime, "hour") == instante_obj)

# --- Vista rápida (top 10) ---
knitr::kable(
  datos_instante |>
    st_drop_geometry() |>
    select(datetime, id_station, id_parameter, value, unit, nom_estac, easting, northing) |>
    slice_head(n = 10),
  caption = paste("REDMET —", format(instante_obj, "%Y-%m-%d %H:%M"), TZ)
)
REDMET — 2024-11-25 12:00 America/Mexico_City
datetime id_station id_parameter value unit nom_estac easting northing
2024-11-25 12:00:00 ACO RH 21.0 6 Acolman 509226.0 2171149
2024-11-25 12:00:00 ACO TMP 18.3 5 Acolman 509226.0 2171149
2024-11-25 12:00:00 ACO WDR 272.0 4 Acolman 509226.0 2171149
2024-11-25 12:00:00 ACO WSP 3.1 3 Acolman 509226.0 2171149
2024-11-25 12:00:00 AJM RH 14.0 6 Ajusco Medio 478170.7 2130955
2024-11-25 12:00:00 AJM TMP 19.3 5 Ajusco Medio 478170.7 2130955
2024-11-25 12:00:00 AJM WDR NA 4 Ajusco Medio 478170.7 2130955
2024-11-25 12:00:00 AJM WSP NA 3 Ajusco Medio 478170.7 2130955
2024-11-25 12:00:00 AJU RH 17.0 6 Ajusco 482901.0 2117907
2024-11-25 12:00:00 AJU TMP 15.9 5 Ajusco 482901.0 2117907

Gráficos Exploratorios

# --- Pares (x, y, variable) para un parámetro ---
pairs_param <- function(df, param, etiqueta) {
  d <- df |>
    sf::st_drop_geometry() |>
    dplyr::filter(id_parameter == param, !is.na(value)) |>
    dplyr::select(easting, northing, value)

  if (nrow(d) == 0) {
    message("Sin datos para ", param); return(invisible(NULL))
  }

  # Etiquetas en los paneles (x, y, nombre de la variable)
  pairs(~ easting + northing + value, data = d,
        labels = c("x", "y", etiqueta),
        pch = 1)
}

# --- Llamadas (una figura por variable) ---
pairs_param(datos_instante, "TMP", "Temperatura (°C)")

pairs_param(datos_instante, "RH",  "Humedad relativa (%)")

pairs_param(datos_instante, "WSP", "Velocidad viento (m/s)")

Temperatura (TMP)

En la primera figura no se aprecia una tendencia espacial marcada para la temperatura. No obstante, se observan valores claramente periféricos que podrían influir de forma desproporcionada en los ajustes. Por prudencia, se repite la exploración excluyéndolas y se comparan los resultados.

Humedad relativa (RH)

En el caso de la humedad relativa se sugiere la presencia de una tendencia positiva: los valores tienden a aumentar a medida que cambia la posición espacial. Esto amerita modelar y/o remover la tendencia antes de ajustar el variograma.

quitar_estaciones_extremas_por_param <- function(df_instante, param = "TMP", n_extremos = 2) {
  df_instante <- df_instante |> dplyr::mutate(value = as.numeric(value))

  tmp_ordenado <- df_instante |>
    dplyr::filter(id_parameter == param, !is.na(value)) |>
    dplyr::arrange(value) |>
    sf::st_drop_geometry()

  if (nrow(tmp_ordenado) == 0L) {
    warning("No hay datos para el parámetro: ", param)
    return(list(df_filtrado = df_instante, ids_quitados = character(0)))
  }

  n_small <- min(n_extremos, nrow(tmp_ordenado))
  n_large <- min(n_extremos, nrow(tmp_ordenado))

  ids_bajas <- dplyr::slice_head(tmp_ordenado, n = n_small)$id_station
  ids_altas <- dplyr::slice_tail(tmp_ordenado, n = n_large)$id_station
  ids_quitar <- unique(c(ids_bajas, ids_altas))

  message("Estaciones removidas por ", param, " extrema: ",
          paste(ids_quitar, collapse = ", "))

  df_filtrado <- df_instante |>
    dplyr::filter(!id_station %in% ids_quitar)

  list(df_filtrado = df_filtrado, ids_quitados = ids_quitar)
}

# Ejecuta:
res <- quitar_estaciones_extremas_por_param(datos_instante, param = "TMP", n_extremos = 2)
Estaciones removidas por TMP extrema: AJU, INN, MER, FAC
datos_instante_sin_tmp_extrema <- res$df_filtrado

# Gráficos exploratorios nuevamente:
pairs_param(datos_instante_sin_tmp_extrema, "TMP", "Temperatura (°C)")

pairs_param(datos_instante_sin_tmp_extrema, "RH",  "Humedad relativa (%)")

pairs_param(datos_instante_sin_tmp_extrema, "WSP", "Velocidad viento (m/s)")

Temperatura (TMP)

Los puntos parecen bastante dispersos y sin un patrón evidente. Las diferencias entre valores son pequeñas (rango estrecho, ~18.8 – 20 °C), lo cual sugiere que quitar las estaciones extremas eliminó los posibles outliers que distorsionaban la distribución. Esto indica que el supuesto de estacionariedad (media constante en el espacio) puede ser razonable para continuar con el análisis geoestadístico.

Humedad relativa (RH)

Aunque sigue habiendo cierta dispersión, no se observa ya una tendencia fuerte. La dispersión es mayor que en TMP, pero no concentrada en un solo sector del área, lo cual también sugiere ausencia de tendencia clara.

Esto permite asumir que la variabilidad de la humedad es más local que global, por lo que el variograma debería capturar estructuras de dependencia a pequeña escala.

# --- helper: prueba tendencia y entrega residuales con media≈0 ---
fit_trend <- function(sf_obj, param) {
  d <- sf::st_drop_geometry(sf_obj) |>
    dplyr::filter(id_parameter == param) |>
    dplyr::mutate(
      z = as.numeric(value),                 # (1) asegurar numérico
      x = as.numeric(easting),
      y = as.numeric(northing)
    ) |>
    dplyr::filter(!is.na(z), !is.na(x), !is.na(y))

  stopifnot(nrow(d) > 3)

  # estandarizar x,y para estabilidad numérica
  d$x_sc <- as.numeric(scale(d$x))
  d$y_sc <- as.numeric(scale(d$y))

  # (3) si x o y quedan (casi) constantes, evitamos términos que exploten
  var_x <- stats::var(d$x_sc, na.rm = TRUE)
  var_y <- stats::var(d$y_sc, na.rm = TRUE)
  usar_x <- !is.na(var_x) && var_x > .Machine$double.eps
  usar_y <- !is.na(var_y) && var_y > .Machine$double.eps

  f1 <- as.formula(
    paste0("z ~ ",
           paste(c(if (usar_x) "x_sc", if (usar_y) "y_sc"), collapse = " + "))
  )
  f2 <- as.formula(
    paste0("z ~ ",
           paste(c(if (usar_x) "x_sc", if (usar_y) "y_sc",
                   if (usar_x) "I(x_sc^2)", if (usar_y) "I(y_sc^2)",
                   if (usar_x && usar_y) "I(x_sc*y_sc)"), collapse = " + "))
  )

  # modelos candidatos (2) misma data en todos
  m0 <- lm(z ~ 1, data = d) #sin tendencia
  m1 <- if (usar_x || usar_y) lm(f1, data = d) else m0  #Lineal
  m2 <- if (usar_x || usar_y) lm(f2, data = d) else m1  # Cuadrático

  # comparación (ANOVA parcial + AIC)
  a01 <- tryCatch(anova(m0, m1), error = function(e) NA)
  a12 <- tryCatch(anova(m1, m2), error = function(e) NA)
  aic <- c(m0 = AIC(m0), m1 = AIC(m1), m2 = AIC(m2))

  # regla de decisión
  best <- m0; best_name <- "m0: sin tendencia"
  p01 <- if (is.list(a01)) a01$`Pr(>F)`[2] else 1
  p12 <- if (is.list(a12)) a12$`Pr(>F)`[2] else 1

  if ((usar_x || usar_y) && p01 < 0.05 && (AIC(m1) <= AIC(m0) - 2)) {
    best <- m1 ; best_name <- "m1: lineal (x,y)"
    if (p12 < 0.05 && (AIC(m2) <= AIC(m1) - 2)) {
      best <- m2 ; best_name <- "m2: cuadrático (x,y,x²,y²,xy)"
    }
  }

  # residuales y chequeo de media
  res <- resid(best)
  res <- res - mean(res, na.rm = TRUE)  # centrado exacto a 0

  list(
    data     = d,
    model    = best,
    which    = best_name,
    anova01  = a01,
    anova12  = a12,
    AIC      = aic,
    residual = res
  )
}

# --- correr para cada variable ---
fit_TMP <- fit_trend(datos_instante, "TMP")
fit_RH  <- fit_trend(datos_instante, "RH")
fit_WSP <- fit_trend(datos_instante, "WSP")

# --- resumen corto en consola ---
cat("\nTMP →", fit_TMP$which, " | AIC:", paste(names(fit_TMP$AIC), round(fit_TMP$AIC,1), collapse="  "), "\n")

TMP → m0: sin tendencia  | AIC: m0 60.1  m1 63.1  m2 49.4 
cat("RH  →", fit_RH$which,  " | AIC:", paste(names(fit_RH$AIC),  round(fit_RH$AIC,1),  collapse="  "), "\n")
RH  → m1: lineal (x,y)  | AIC: m0 125  m1 113.2  m2 115 
cat("WSP →", fit_WSP$which, " | AIC:", paste(names(fit_WSP$AIC), round(fit_WSP$AIC,1), collapse="  "), "\n")
WSP → m1: lineal (x,y)  | AIC: m0 80  m1 54.3  m2 50.5 
# --- diagnóstico visual rápido: residuales vs x,y (deberían no mostrar tendencia) ---
diag_plot <- function(fit, titulo) {
  par(mfrow = c(1,2), mar = c(4,4,2,1))
  plot(fit$data$x, fit$residual, pch=1, xlab="x (Easting)",  ylab="Residual", main=paste(titulo,"resid vs x"))
  abline(h=0, lty=2, col="grey50")
  plot(fit$data$y, fit$residual, pch=1, xlab="y (Northing)", ylab="Residual", main=paste(titulo,"resid vs y"))
  abline(h=0, lty=2, col="grey50")
  par(mfrow = c(1,1))
}
diag_plot(fit_TMP, "TMP")

diag_plot(fit_RH,  "RH")

diag_plot(fit_WSP, "WSP")

# --- residuales listos para variograma (media≈0) ---
res_TMP <- dplyr::mutate(fit_TMP$data, resid = fit_TMP$residual)
res_RH  <- dplyr::mutate(fit_RH$data,  resid = fit_RH$residual)
res_WSP <- dplyr::mutate(fit_WSP$data, resid = fit_WSP$residual)

Temperatura (TMP)

Se evaluó la presencia de tendencia espacial ajustando tres modelos, m0 sin tendencia; m1 lineal y m2 cuadrático. La selección consideró ANOVA secuencial y AIC.

Aunque m2 mejora moderadamente frente a m0 (()AIC ()), las pruebas de significancia no avalaron incluir términos de tendencia. Se elige m0 (sin tendencia). Para el instante analizado, no se evidencia un gradiente espacial relevante de temperatura tras retirar las estaciones con valores extremos. El supuesto de estacionariedad en media es razonable.

Humedad relativa (RH)

La comparación m1 vs. m0 muestra una mejora clara (()AIC ()), mientras que m2 vs. m1 aporta ganancia marginal (()AIC ()). Se selecciona m1 (tendencia lineal en (x,y)).

RH presenta un gradiente espacial lineal para el instante analizado; no se justifica añadir curvatura. La variabilidad residual es predominantemente local.

Con esto, entonces preparamos los datos para el variograma:

# --- TMP: usar valores (sin tendencia) ---
df_tmp <- fit_TMP$data |>
  dplyr::transmute(x = x, y = y, z = z)  # z = TMP

# --- RH y WSP: usar residuales (con tendencia lineal) ---
df_rh  <- fit_RH$data  |>
  dplyr::transmute(x = x, y = y, z = fit_RH$residual)

df_wsp <- fit_WSP$data |>
  dplyr::transmute(x = x, y = y, z = fit_WSP$residual)

# Función auxiliar para parámetros de binning (en metros)
variogram_binning <- function(x, y, prop_cutoff = 0.6, nbins = 15) {
  # distancia máxima observable
  Dmax <- as.numeric(max(dist(cbind(x, y))))
  cutoff <- prop_cutoff * Dmax
  width  <- cutoff / nbins
  list(cutoff = cutoff, width = width)
}

Semivariograma Empírico

library(gstat)

# --- Definir variogram_binning  ---
variogram_binning <- function(x, y, prop_cutoff = 0.6, nbins = 15) {
  Dmax <- as.numeric(max(stats::dist(cbind(x, y))))
  cutoff <- prop_cutoff * Dmax
  width  <- cutoff / nbins
  list(cutoff = cutoff, width = width)
}

# --- Parámetros por variable ---
p_tmp <- variogram_binning(df_tmp$x, df_tmp$y, 0.6, 15)
p_rh  <- variogram_binning(df_rh$x,  df_rh$y,  0.6, 15)
p_wsp <- variogram_binning(df_wsp$x, df_wsp$y, 0.6, 15)

# =========================
#   Semivariograma clásico
# =========================
vg_TMP <- variogram(z ~ 1, ~ x + y, data = df_tmp,
                    cutoff = p_tmp$cutoff, width = p_tmp$width)
vg_RH  <- variogram(z ~ 1, ~ x + y, data = df_rh,
                    cutoff = p_rh$cutoff,  width = p_rh$width)
vg_WSP <- variogram(z ~ 1, ~ x + y, data = df_wsp,
                    cutoff = p_wsp$cutoff, width = p_wsp$width)

par(mfrow = c(1,3))
plot(vg_TMP, main = "TMP — semivariograma (clásico)")

plot(vg_RH,  main = "RH  — semivariograma (clásico, residuales)")

par(mfrow = c(1,1))

Semivariograma resistente a datos atípicos

# =========================
#   Semivariograma robusto
# =========================
vg_TMP_rob <- variogram(z ~ 1, ~ x + y, data = df_tmp,
                        cutoff = p_tmp$cutoff, width = p_tmp$width,
                        cressie = TRUE)
vg_RH_rob  <- variogram(z ~ 1, ~ x + y, data = df_rh,
                        cutoff = p_rh$cutoff,  width = p_rh$width,
                        cressie = TRUE)
vg_WSP_rob <- variogram(z ~ 1, ~ x + y, data = df_wsp,
                        cutoff = p_wsp$cutoff, width = p_wsp$width,
                        cressie = TRUE)

par(mfrow = c(1,3))
plot(vg_TMP_rob, main = "TMP — semivariograma (robusto Cressie)")

plot(vg_RH_rob,  main = "RH  — semivariograma (robusto, residuales)")

par(mfrow = c(1,1))

Anexos

# ----- Conversión de estaciones MX a coordenadas planas (UTM) -----

#Paquetes
libs <- c("readr","dplyr","stringr","sf")
inst <- setdiff(libs, rownames(installed.packages()))
if (length(inst)) install.packages(inst, quiet = TRUE)
invisible(lapply(libs, library, character.only = TRUE))

#Ruta del CSV 
ruta_csv <- "C:/Users/sarar/Documents/Maestría Estadística/Datos/cat_estacion.csv"

stopifnot(file.exists(ruta_csv))

# Detección simple del separador (coma ; tab)
primera <- readr::read_lines(ruta_csv, n_max = 1)
delim <- if (grepl(";", primera)) ";" else if (grepl("\t", primera)) "\t" else ","

#Leer CSV
df <- readr::read_delim(ruta_csv, delim = delim, show_col_types = FALSE, trim_ws = TRUE)

# Verifica columnas mínimas
req <- c("cve_estac","nom_estac","longitud","latitud")
faltan <- setdiff(req, names(df))
if (length(faltan)) stop("Faltan columnas: ", paste(faltan, collapse = ", "))

# Limpiar/normalizar coordenadas
rm_dots <- function(x) as.numeric(gsub("\\.", "", as.character(x)))

scale_into_range <- function(x, low, high) {
  x <- as.numeric(x)
  for (i in 1:8) {
    ok <- is.na(x) | (x >= low & x <= high)
    if (all(ok)) break
    x[!ok] <- x[!ok] / 10
  }
  x
}

df <- df %>%
  mutate(
    lon_raw = rm_dots(longitud),
    lat_raw = rm_dots(latitud),
    # corrige signo de longitud (debe ser negativa en MX)
    lon = ifelse(!is.na(lon_raw) & lon_raw > 0, -lon_raw, lon_raw),
    lat = lat_raw,
    # escalar a rangos válidos
    lon = scale_into_range(lon, -180, -86),   # primero deja <= -86
    lon = pmin(lon, -86),                     # recorte suave por si quedó -85.x
    lon = pmax(lon, -180),
    lat = scale_into_range(lat, 14, 90),      # sube hasta >= 14
    lat = pmin(lat, 33),                      # recorte suave por si quedó >33
    lat = pmax(lat, 14)
  )

# A sf WGS84
g <- sf::st_as_sf(df, coords = c("lon","lat"), crs = 4326, remove = FALSE)

# Proyección a UTM: zona automática por longitud (11..16 Norte)
zonas  <- floor((df$lon + 180)/6) + 1
zonas[!(zonas %in% 11:16)] <- NA_integer_
epsg   <- ifelse(is.na(zonas), NA_integer_, 32600 + zonas)

g$zona_utm <- zonas
g$epsg_utm <- epsg

# Separar por EPSG válido
parts <- split(g, g$epsg_utm, drop = TRUE)

mk_xy <- function(sfx, crs) {
  tr <- sf::st_transform(sfx, crs)
  xy <- sf::st_coordinates(tr)
  cbind(sf::st_drop_geometry(tr),
        Easting = xy[,1], Northing = xy[,2],
        sistema = paste0("UTM_", crs - 32600, "N"), epsg = crs)
}

if (length(parts)) {
  out_list <- mapply(mk_xy, parts, as.integer(names(parts)), SIMPLIFY = FALSE)
  out <- dplyr::bind_rows(out_list)
} else {
  # Respaldo: todo a UTM 14N (CDMX) si algo falló
  crs_fallback <- 32614
  tr <- sf::st_transform(g, crs_fallback)
  xy <- sf::st_coordinates(tr)
  out <- cbind(sf::st_drop_geometry(tr),
               Easting = xy[,1], Northing = xy[,2],
               sistema = "UTM_14N", epsg = crs_fallback)
}

# Salida ordenada
salida <- out %>%
  dplyr::select(
    cve_estac, nom_estac, longitud, latitud, alt,
    Easting, Northing, sistema, epsg
  )

readr::write_csv(salida, "estaciones_planas.csv")
message("Guardado: estaciones_planas.csv (filas: ", nrow(salida), ")")
Guardado: estaciones_planas.csv (filas: 69)
plot(salida$Easting, salida$Northing, asp = 1,
     xlab = "Easting (m)", ylab = "Northing (m)",
     main = "Estaciones UTM 14N — CDMX")

library(dplyr)
library(leaflet)
library(htmltools)
library(sf)
library(purrr)

# --- Partimos de 'salida' ya creado arriba (Easting, Northing, epsg, sistema, etc.) ---

# 1) Construir un sf por EPSG UTM y transformar a WGS84 (4326)
split_epsg <- split(salida, salida$epsg, drop = TRUE)

to_ll <- function(df_epsg) {
  crs_utm <- unique(na.omit(df_epsg$epsg))
  if (length(crs_utm) != 1L) stop("EPSG inconsistente en el grupo.")
  sf_utm <- st_as_sf(df_epsg, coords = c("Easting","Northing"), crs = crs_utm, remove = FALSE)
  sf_ll  <- st_transform(sf_utm, 4326)
  coords <- st_coordinates(sf_ll)
  sf_ll |>
    st_drop_geometry() |>
    mutate(lon = coords[,1], lat = coords[,2])
}

pts <- map_dfr(split_epsg, to_ll)

# 2) Etiquetas: usa un nombre si existe; si no, “Punto i”
safe_utf8 <- function(x) enc2utf8(iconv(ifelse(is.na(x), "", x), from = "", to = "UTF-8"))
col_bonita <- intersect(c("nom_estac","NOMBRE","nombre","ESTACION","estacion"), names(pts))
pts <- pts |>
  mutate(
    label_txt = if (length(col_bonita)) safe_utf8(.data[[col_bonita[1]]]) else paste0("Punto ", row_number()),
    popup_txt = paste0(
      "<b>", htmlEscape(label_txt), "</b><br/>",
      "EPSG: ", epsg, " (", sistema, ")<br/>",
      "Easting: ", format(Easting, big.mark = " ", scientific = FALSE), " m<br/>",
      "Northing: ", format(Northing, big.mark = " ", scientific = FALSE), " m<br/>",
      "Lon/Lat: ", sprintf('%.6f', lon), ", ", sprintf('%.6f', lat)
    )
  )

# 3) Mapa Leaflet (tiles en Web Mercator; puntos en WGS84)
leaflet(pts, options = leafletOptions(minZoom = 4)) |>
  addProviderTiles("CartoDB.Positron") |>
  addCircleMarkers(
    lng = ~lon, lat = ~lat,
    radius = 6, stroke = TRUE, weight = 1, fillOpacity = 0.9,
    label = ~htmlEscape(label_txt),
    popup = ~HTML(popup_txt)
  ) |>
  fitBounds(
    lng1 = min(pts$lon, na.rm = TRUE), lat1 = min(pts$lat, na.rm = TRUE),
    lng2 = max(pts$lon, na.rm = TRUE), lat2 = max(pts$lat, na.rm = TRUE)
  )