knitr::include_graphics("Img/Rplot01.png")
PARCIAL 2: Diagnóstico espacio–temporal 2024 de variables meteorológicas de la REDMET (México)
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.
Objetivo general
Entregar 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).
Variables y unidades
- Temperatura del aire (TMP) — grados Celsius (°C)
- Humedad relativa (RH) — porcentaje (%)
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.
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:
- Integrar y depurar las series horarias 2024 (detección de faltantes/atípicos y estandarización de unidades).
- Localizar y mapear las estaciones, superponiéndolas al área de interés para identificar gradientes y zonas con condiciones atípicas.
- Caracterizar patrones temporales (ciclo diurno/semanal) y patrones espaciales básicos (mapas de valores y semivariograma exploratorio).
- 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
# --- 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),
# --- Normalizar formato "DD/MM/YYYY H:MM" → "DD/MM/YYYY HH:MM"
date_clean = date |>
stringr::str_replace("^([0-9]{2}/[0-9]{2}/[0-9]{4}) ([0-9]):",
"\\1 0\\2:"),
# --- Ahora sí se puede parsear correctamente
datetime = lubridate::dmy_hm(date_clean, tz = TZ),
date = as.Date(datetime)
)Warning: There was 1 warning in `mutate()`.
ℹ In argument: `datetime = lubridate::dmy_hm(date_clean, tz = TZ)`.
Caused by warning:
! 40880 failed to parse.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_crs <- unique(na.omit(datos_merged$epsg))[1]
datos_sf <- st_as_sf(datos_merged,
coords=c("easting","northing"),
crs=epsg_crs, remove=FALSE)
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)
datos_instante |>
dplyr::select(datetime, id_station, id_parameter, value, unit, nom_estac, easting, northing)Simple feature collection with 112 features and 8 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 460020.4 ymin: 2117907 xmax: 511970 ymax: 2180751
Projected CRS: WGS 84 / UTM zone 14N
# A tibble: 112 × 9
datetime id_station id_parameter value unit nom_estac easting
<dttm> <chr> <chr> <dbl> <dbl> <chr> <dbl>
1 2024-11-25 12:00:00 ACO RH 21 6 Acolman 509226.
2 2024-11-25 12:00:00 ACO TMP 18.3 5 Acolman 509226.
3 2024-11-25 12:00:00 ACO WDR 272 4 Acolman 509226.
4 2024-11-25 12:00:00 ACO WSP 3.1 3 Acolman 509226.
5 2024-11-25 12:00:00 AJM RH 14 6 Ajusco Medio 478171.
6 2024-11-25 12:00:00 AJM TMP 19.3 5 Ajusco Medio 478171.
7 2024-11-25 12:00:00 AJM WDR NA 4 Ajusco Medio 478171.
8 2024-11-25 12:00:00 AJM WSP NA 3 Ajusco Medio 478171.
9 2024-11-25 12:00:00 AJU RH 17 6 Ajusco 482901.
10 2024-11-25 12:00:00 AJU TMP 15.9 5 Ajusco 482901.
# ℹ 102 more rows
# ℹ 2 more variables: northing <dbl>, geometry <POINT [m]># --- Vista rápida (top 10) ---
knitr::kable(
datos_instante |>
st_drop_geometry() |>
dplyr::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)
)| 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 |
1. Análisis geoestadístico univariado y bivariado de al menos dos variables espaciales:
1.1 Análisis y modelamiento de cada una de las medias.
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 (%)")
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.
Detección y remoción de outliers
calc_outliers <- function(df_var) {
vals <- df_var$value
Q1 <- quantile(vals, 0.25, na.rm = TRUE)
Q3 <- quantile(vals, 0.75, na.rm = TRUE)
IQR <- Q3 - Q1
lower <- Q1 - 1.5 * IQR
upper <- Q3 + 1.5 * IQR
df_var |> filter(value < lower | value > upper) |> pull(id_station)
}
# --- TMP ---
df_tmp <- datos_instante |>
filter(id_parameter == "TMP") |>
st_drop_geometry() |>
mutate(value = as.numeric(value))
ids_tmp_out <- calc_outliers(df_tmp)
cat("Outliers TMP:", ifelse(length(ids_tmp_out)==0, "ninguno", paste(ids_tmp_out, collapse=", ")), "\n")Outliers TMP: AJU, FAC, INN datos_sin_tmp <- datos_instante |> filter(!(id_station %in% ids_tmp_out))
# --- RH ---
df_rh <- datos_sin_tmp |>
filter(id_parameter == "RH") |>
st_drop_geometry() |>
mutate(value = as.numeric(value))
ids_rh_out <- calc_outliers(df_rh)
cat("Outliers RH:", ifelse(length(ids_rh_out)==0, "ninguno", paste(ids_rh_out, collapse=", ")), "\n")Outliers RH: ninguno datos_limpios <- datos_sin_tmp |> filter(!(id_station %in% ids_rh_out))
# --- Gráficos finales ---
pairs_param(datos_limpios, "TMP", "Temperatura (°C)")
pairs_param(datos_limpios, "RH", "Humedad relativa (%)")
A partir de la evaluación gráfica inicial y del criterio intercuartílico (IQR), se quitaron las estaciones AJU, FAC e INN las cuales presentaban mediciones anomalas para la temperatura. En el caso de laTemperatura (TMP) los puntos parecen bastante dispersos y sin un patrón evidente, 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.
En la 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.
Modelamiento de la media
## --- función para preparar un df por variable ---
prep_df_utm <- function(datos,param){
datos |> filter(id_parameter==param) |> st_drop_geometry() |>
transmute(
x = easting,
y = northing,
z = as.numeric(value)
) |> filter(!is.na(z))
}
df_tmp <- prep_df_utm(datos_limpios,"TMP")
df_rh <- prep_df_utm(datos_limpios,"RH")1. Temperatura (TMP)
# Modelos para la media
m0_tmp <- lm(z ~ 1, data = df_tmp)
m1_tmp <- lm(z ~ x + y, data = df_tmp)
m2_tmp <- lm(z ~ x + y + I(x^2) + I(y^2) + I(x*y), data = df_tmp)
summary(m0_tmp)
Call:
lm(formula = z ~ 1, data = df_tmp)
Residuals:
Min 1Q Median 3Q Max
-0.95333 -0.15333 -0.05333 0.34667 0.84667
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.253 0.139 138.5 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5383 on 14 degrees of freedomsummary(m1_tmp)
Call:
lm(formula = z ~ x + y, data = df_tmp)
Residuals:
Min 1Q Median 3Q Max
-0.5471 -0.3917 -0.1247 0.3596 0.7757
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.543e+01 1.685e+01 3.883 0.00218 **
x -6.306e-06 1.186e-05 -0.532 0.60449
y -2.004e-05 7.907e-06 -2.534 0.02623 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4558 on 12 degrees of freedom
Multiple R-squared: 0.3854, Adjusted R-squared: 0.283
F-statistic: 3.763 on 2 and 12 DF, p-value: 0.05388summary(m2_tmp)
Call:
lm(formula = z ~ x + y + I(x^2) + I(y^2) + I(x * y), data = df_tmp)
Residuals:
Min 1Q Median 3Q Max
-0.5760 -0.2657 -0.1323 0.4042 0.5621
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.346e+03 2.985e+03 -1.456 0.179
x 4.870e-04 1.734e-03 0.281 0.785
y 3.965e-03 2.633e-03 1.506 0.166
I(x^2) -8.188e-11 1.302e-09 -0.063 0.951
I(y^2) -9.028e-10 5.651e-10 -1.598 0.145
I(x * y) -1.944e-10 8.953e-10 -0.217 0.833
Residual standard error: 0.4563 on 9 degrees of freedom
Multiple R-squared: 0.5381, Adjusted R-squared: 0.2815
F-statistic: 2.097 on 5 and 9 DF, p-value: 0.1578anova(m0_tmp, m1_tmp, m2_tmp)Analysis of Variance Table
Model 1: z ~ 1
Model 2: z ~ x + y
Model 3: z ~ x + y + I(x^2) + I(y^2) + I(x * y)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 14 4.0573
2 12 2.4935 2 1.56386 3.7549 0.0652 .
3 9 1.8742 3 0.61929 0.9913 0.4398
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Se evaluaron tres modelos para la media espacial:
- m0: modelo sin tendencia
- m1: tendencia lineal en ((x, y))
- m2: tendencia cuadrática en ((x, y))
La comparación mediante ANOVA mostró que el modelo lineal no mejora significativamente al modelo sin tendencia ((p = 0.065)).
Asimismo, los términos cuadráticos incluidos en (m2) no aportaron una mejora estadística respecto a (m1) ((p = 0.44)).
Aunque dentro del modelo lineal el coeficiente asociado a la coordenada (y) aparece marginalmente significativo, el efecto global del modelo completo no supera al de (m0). Esto indica que, tras la eliminación de las estaciones con valores atípicos (AJU, FAC, INN), no se observa un gradiente espacial relevante en la temperatura para el instante analizado.
En consecuencia, se selecciona el modelo sin tendencia (m0) y se asume estacionariedad en la media para TMP.
RESIDUOS PARA EL VARIOGRAMA
df_tmp$resid <- residuals(m0_tmp)
df_tmp$resid <- df_tmp$resid - mean(df_tmp$resid, na.rm = TRUE) # centrar exacto en 02. Humedad relativa (RH)
m0_rh <- lm(z ~ 1, data = df_rh)
m1_rh <- lm(z ~ x + y, data = df_rh)
m2_rh <- lm(z ~ x + y + I(x^2) + I(y^2) + I(x*y), data = df_rh)
summary(m0_rh)
Call:
lm(formula = z ~ 1, data = df_rh)
Residuals:
Min 1Q Median 3Q Max
-6.8889 -2.6389 -0.8889 3.6111 7.1111
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 20.8889 0.9249 22.58 4.07e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.924 on 17 degrees of freedomsummary(m1_rh)
Call:
lm(formula = z ~ x + y, data = df_rh)
Residuals:
Min 1Q Median 3Q Max
-4.7402 -2.5088 0.1192 2.2854 7.0244
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.376e+02 1.260e+02 -1.886 0.0788 .
x 1.788e-04 7.552e-05 2.367 0.0318 *
y 7.937e-05 5.662e-05 1.402 0.1813
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.388 on 15 degrees of freedom
Multiple R-squared: 0.3424, Adjusted R-squared: 0.2547
F-statistic: 3.905 on 2 and 15 DF, p-value: 0.04314summary(m2_rh)
Call:
lm(formula = z ~ x + y + I(x^2) + I(y^2) + I(x * y), data = df_rh)
Residuals:
Min 1Q Median 3Q Max
-4.8526 -1.3401 0.2186 1.5106 5.2337
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.335e+04 2.163e+04 -1.079 0.302
x 2.225e-02 1.128e-02 1.973 0.072 .
y 1.651e-02 1.976e-02 0.835 0.420
I(x^2) -1.050e-08 7.785e-09 -1.349 0.202
I(y^2) -3.194e-09 4.356e-09 -0.733 0.477
I(x * y) -5.447e-09 4.575e-09 -1.191 0.257
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.13 on 12 degrees of freedom
Multiple R-squared: 0.5509, Adjusted R-squared: 0.3637
F-statistic: 2.944 on 5 and 12 DF, p-value: 0.05822anova(m0_rh, m1_rh, m2_rh)Analysis of Variance Table
Model 1: z ~ 1
Model 2: z ~ x + y
Model 3: z ~ x + y + I(x^2) + I(y^2) + I(x * y)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 17 261.78
2 15 172.15 2 89.625 4.5738 0.03338 *
3 12 117.57 3 54.582 1.8570 0.19069
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La comparación mediante ANOVA mostró que el modelo lineal mejora significativamente al modelo sin tendencia ((p = 0.033)), indicando la presencia de un gradiente espacial suave en la humedad relativa.
Por otro lado, la inclusión de términos cuadráticos en (m2) no aportó una mejora estadística respecto al modelo lineal ((p = 0.19)), y ninguno de los términos no lineales resultó significativo.
En consecuencia, se selecciona el modelo lineal (m1) para describir la media espacial de RH.
RESIDUOS PARA EL VARIOGRAMA
df_rh$resid <- residuals(m1_rh)
df_rh$resid <- df_rh$resid - mean(df_rh$resid, na.rm = TRUE)1.2 Selección de los semivariogramas empíricos individuales, más adecuados para la estimación del modelo teórico
Variograma cloud
# variograma cloud de residuales
vg_TMP_cloud <- variogram(resid ~ 1, ~ x + y, data = df_tmp, cloud = TRUE)
vg_RH_cloud <- variogram(resid ~ 1, ~ x + y, data = df_rh, cloud = TRUE)
par(mfrow = c(1,2))
plot(vg_TMP_cloud, main = "TMP — Cloud (residuales)")
plot(vg_RH_cloud, main = "RH — Cloud (residuales)")
par(mfrow = c(1,1))El variograma cloud de TMP muestra una nube compacta, sin valores atípicos ni semivarianzas extremas. La variación entre pares es baja y estable, indicando que el cálculo del variograma clásico es adecuado.
En contraste, el variograma cloud de RH presenta una dispersión mucho mayor, con semivarianzas que alcanzan valores altos incluso a distancias intermedias. Esto refleja variabilidad local alta y posibles pares influyentes. Por ello, resulta pertinente complementar el análisis con un estimador robusto.
Semivariograma Clásico (binned)
variogram_binning <- function(df, prop_cutoff = 0.6, nbins = 15) {
Dmax <- max(dist(df[, c("x", "y")]))
cutoff <- prop_cutoff * Dmax
width <- cutoff / nbins
list(cutoff = cutoff, width = width)
}
p_tmp <- variogram_binning(df_tmp)
p_rh <- variogram_binning(df_rh)
# Binned variogram (clásico)
vg_TMP <- variogram(resid ~ 1, ~ x + y, data = df_tmp,
cutoff = p_tmp$cutoff, width = p_tmp$width)
vg_RH <- variogram(resid ~ 1, ~ x + y, data = df_rh,
cutoff = p_rh$cutoff, width = p_rh$width)
par(mfrow = c(1,2))
plot(vg_TMP, main = "TMP — Semivariograma clásico")
plot(vg_RH, main = "RH — Semivariograma clásico")
par(mfrow = c(1,1))El semivariograma clásico de TMP muestra una estructura espacial suave, con semivarianzas muy bajas. Esto confirma que la temperatura es un proceso relativamente estable, sin valores extremos.
En contraste, el semivariograma clásico de RH (calculado sobre los residuales del modelo lineal) revela una dependencia espacial más marcada y una variabilidad considerablemente mayor. La forma general sigue siendo creciente, aunque con mayor dispersión, lo que refleja la naturaleza más heterogénea de la humedad relativa.
Semivariograma resistente a datos atípicos
vg_TMP_rob <- variogram(resid ~ 1, ~ x + y, data = df_tmp,
cutoff = p_tmp$cutoff, width = p_tmp$width,
cressie = TRUE)
vg_RH_rob <- variogram(resid ~ 1, ~ x + y, data = df_rh,
cutoff = p_rh$cutoff, width = p_rh$width,
cressie = TRUE)
par(mfrow = c(1,2))
plot(vg_TMP_rob, main = "TMP — Semivariograma robusto")
plot(vg_RH_rob, main = "RH — Semivariograma robusto")
par(mfrow = c(1,1))En el caso de TMP, el semivariograma robusto (Cressie–Hawkins) presenta una forma casi idéntica al clásico. Esto confirma la ausencia de valores atípicos influyentes y la estabilidad de la estructura espacial de segundo orden. Por tanto, el estimador clásico es suficiente para ajustar un modelo teórico.
Para RH, el variograma robusto suaviza la variabilidad inicial y atenúa la influencia de pares de observaciones con diferencias extremas. Esto produce una estructura espacial más clara, con un aumento más regular de la semivarianza con la distancia. Debido a la mayor heterogeneidad local de la humedad relativa, el estimador robusto resulta más adecuado para la selección del modelo teórico.
1.3 Estimación de los semivariogramas individuales teóricos a sentimiento
El ajuste visual del semivariograma para la temperatura (TMP) se realizó mediante la herramienta eyefit() del paquete geoR.
El modelo que mejor representó la forma del variograma empírico fue esférico.
Parámetros seleccionados (aprox.):
- Nugget ((^2)) ≈ 0.10
- Sill parcial ((^2)) ≈ 0.215
- Rango (()) ≈ 18 km
knitr::include_graphics("Img/Semivariograma_tem.png")
El modelo esférico captura razonablemente el crecimiento inicial del variograma y alcanza un nivel cercano a la meseta alrededor de 18 km, consistente con la estructura espacial suave observada para TMP. La presencia de un nugget pequeño concuerda con una medición relativamente estable y sin discontinuidades fuertes entre estaciones cercanas.
Para la humedad relativa (RH), el semivariograma presenta una mayor dispersión y valores altos de semivarianza. El modelo que mejor representó la forma creciente y cóncava del variograma empírico fue el exponencial.
Parámetros seleccionados (aprox.):
- Nugget ((^2)) ≈ 1.85
- Sill parcial ((^2)) ≈ 14.7
- Rango (()) ≈ 11.5 km
knitr::include_graphics("img/Semivariograma_RH.png")
El modelo exponencial refleja de manera adecuada la rápida subida inicial del variograma y la aproximación asintótica hacia la meseta. El nugget relativamente grande es coherente con la alta variabilidad local detectada en RH, incluso entre estaciones cercanas.
##1.4. Estimación de los semivariogramas individuales y del modelo teórico usando MCO, MCP, MV, CL.
Vamos a usar geor:
# --- Geodata usando los residuales ---
geo_TMP <- as.geodata(df_tmp,
coords.col = c("x", "y"),
data.col = "resid")
geo_RH <- as.geodata(df_rh,
coords.col = c("x", "y"),
data.col = "resid"
)
var_TMP <- variog(geo_TMP, max.dist = p_tmp$cutoff)variog: computing omnidirectional variogramvar_RH <- variog(geo_RH, max.dist = p_rh$cutoff)variog: computing omnidirectional variogramplot(var_TMP, main = "TMP — Variograma empírico (geoR)")
plot(var_RH, main = "RH — Variograma empírico (geoR)")
Parámetros iniciales según lo eyefit:
# ---- TMP (modelo esférico) ----
ini_TMP <- c(sigmasq = 0.215, phi = 18000)
nug_TMP <- 0.10
# ---- RH (modelo exponencial) ----
ini_RH <- c(sigmasq = 14.7, phi = 11500)
nug_RH <- 1.85Temperatura:
# ---- TMP: MCO ----
fit_TMP_MCO <- variofit(
vario = var_TMP,
ini.cov.pars = ini_TMP,
cov.model = "spherical",
fix.nugget = TRUE,
nugget = nug_TMP,
weights = "equal"
)variofit: covariance model used is spherical
variofit: weights used: equal
variofit: minimisation function used: optim # ---- TMP: MCP (npairs) ----
fit_TMP_npairs <- variofit(
vario = var_TMP,
ini.cov.pars = ini_TMP,
cov.model = "spherical",
fix.nugget = TRUE,
nugget = nug_TMP,
weights = "npairs"
)variofit: covariance model used is spherical
variofit: weights used: npairs
variofit: minimisation function used: optim # ---- TMP: MCP (Cressie) ----
fit_TMP_cressie <- variofit(
vario = var_TMP,
ini.cov.pars = ini_TMP,
cov.model = "spherical",
fix.nugget = TRUE,
nugget = nug_TMP,
weights = "cressie"
)variofit: covariance model used is spherical
variofit: weights used: cressie
variofit: minimisation function used: optim # ---- TMP: Máxima verosimilitud (ML) ----
fit_TMP_ML <- likfit(
geodata = geo_TMP,
coords = geo_TMP$coords,
data = geo_TMP$data,
ini.cov.pars = ini_TMP,
fix.nugget = TRUE,
nugget = nug_TMP,
cov.model = "spherical",
lik.method = "ML"
)kappa not used for the spherical correlation function
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
arguments for the maximisation function.
For further details see documentation for optim.
likfit: It is highly advisable to run this function several
times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.# ---- TMP: REML ----
fit_TMP_REML <- likfit(
geodata = geo_TMP,
coords = geo_TMP$coords,
data = geo_TMP$data,
ini.cov.pars = ini_TMP,
fix.nugget = TRUE,
nugget = nug_TMP,
cov.model = "spherical",
lik.method = "REML"
)kappa not used for the spherical correlation function
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
arguments for the maximisation function.
For further details see documentation for optim.
likfit: It is highly advisable to run this function several
times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.plot(var_TMP$u, var_TMP$v,
xlab = "Distancia (h)",
ylab = "Semivarianza",
main = "TMP — Ajustes del semivariograma")
lines(fit_TMP_MCO, col = 1, lwd = 2)
lines(fit_TMP_npairs, col = 2, lwd = 2)
lines(fit_TMP_cressie, col = 3, lwd = 2)
lines(fit_TMP_ML, col = 4, lwd = 2)
lines(fit_TMP_REML, col = 5, lwd = 2)
legend(
"bottomright",
legend = c("MCO", "MCP npairs", "MCP Cressie", "ML", "REML"),
col = 1:5,
lwd = 2
)
# ---- RH: MCO ----
fit_RH_MCO <- variofit(
vario = var_RH,
ini.cov.pars = ini_RH,
cov.model = "exponential",
fix.nugget = TRUE,
nugget = nug_RH,
weights = "equal"
)variofit: covariance model used is exponential
variofit: weights used: equal
variofit: minimisation function used: optim # ---- RH: MCP (npairs) ----
fit_RH_npairs <- variofit(
vario = var_RH,
ini.cov.pars = ini_RH,
cov.model = "exponential",
fix.nugget = TRUE,
nugget = nug_RH,
weights = "npairs"
)variofit: covariance model used is exponential
variofit: weights used: npairs
variofit: minimisation function used: optim # ---- RH: MCP (Cressie) ----
fit_RH_cressie <- variofit(
vario = var_RH,
ini.cov.pars = ini_RH,
cov.model = "exponential",
fix.nugget = TRUE,
nugget = nug_RH,
weights = "cressie"
)variofit: covariance model used is exponential
variofit: weights used: cressie
variofit: minimisation function used: optim # ---- RH: Máxima verosimilitud (ML) ----
fit_RH_ML <- likfit(
geodata = geo_RH,
coords = geo_RH$coords,
data = geo_RH$data,
ini.cov.pars = ini_RH,
fix.nugget = TRUE,
nugget = nug_RH,
cov.model = "exponential",
lik.method = "ML"
)kappa not used for the exponential correlation function
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
arguments for the maximisation function.
For further details see documentation for optim.
likfit: It is highly advisable to run this function several
times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.# ---- RH: REML ----
fit_RH_REML <- likfit(
geodata = geo_RH,
coords = geo_RH$coords,
data = geo_RH$data,
ini.cov.pars = ini_RH,
fix.nugget = TRUE,
nugget = nug_RH,
cov.model = "exponential",
lik.method = "REML"
)kappa not used for the exponential correlation function
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
arguments for the maximisation function.
For further details see documentation for optim.
likfit: It is highly advisable to run this function several
times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.plot(var_RH$u, var_RH$v,
xlab = "Distancia (h)",
ylab = "Semivarianza",
main = "RH — Ajustes del semivariograma")
lines(fit_RH_MCO, col = 1, lwd = 2)
lines(fit_RH_npairs, col = 2, lwd = 2)
lines(fit_RH_cressie, col = 3, lwd = 2)
lines(fit_RH_ML, col = 4, lwd = 2)
lines(fit_RH_REML, col = 5, lwd = 2)
legend(
"bottomright",
legend = c("MCO", "MCP npairs", "MCP Cressie", "ML", "REML"),
col = 1:5,
lwd = 2
)
Los semivariogramas empíricos presentan alta dispersión, sin evidenciar una silla clara ni una tendencia monótonamente creciente. Esto sugiere que la estructura espacial a la hora analizada es débil o no estacionaria. Por esta razón, los modelos paramétricos (esférico y exponencial) solo logran ajustes aproximados y no necesariamente representan fielmente la variabilidad real.
Semivariogramas recortados:
# diámetro espacial total
Dmax_tmp <- max(dist(df_tmp[, c("x", "y")]))
Dmax_rh <- max(dist(df_rh[, c("x", "y")]))
Dmax_tmp; Dmax_rh[1] 56218.04[1] 60107.79cut_tmp <- 0.40 * Dmax_tmp
cut_rh <- 0.40 * Dmax_rh
# Recalcular variogramas recortados en geoR
var_TMP <- variog(
geo_TMP,
max.dist = cut_tmp
)variog: computing omnidirectional variogramvar_RH <- variog(
geo_RH,
max.dist = cut_rh
)variog: computing omnidirectional variogrampar(mfrow = c(1,2))
plot(var_TMP, main = "TMP — Variograma recortado")
plot(var_RH, main = "RH — Variograma recortado")
par(mfrow = c(1,1))# --- Parámetros iniciales TMP (tomados del eyefit) ---
ini_TMP <- c(sigmasq = 0.215, phi = 18000)
nug_TMP <- 0.10 # nugget fijo
# -------------------------
# 1. MCO (OLS)
# -------------------------
fit_TMP_MCO_REC <- variofit(
vario = var_TMP,
ini.cov.pars = ini_TMP,
cov.model = "spherical",
weights = "equal",
fix.nugget = TRUE,
nugget = nug_TMP
)variofit: covariance model used is spherical
variofit: weights used: equal
variofit: minimisation function used: optim # -------------------------
# 2. MCP (WLS - npairs)
# -------------------------
fit_TMP_npairs_REC <- variofit(
vario = var_TMP,
ini.cov.pars = ini_TMP,
cov.model = "spherical",
weights = "npairs",
fix.nugget = TRUE,
nugget = nug_TMP
)variofit: covariance model used is spherical
variofit: weights used: npairs
variofit: minimisation function used: optim # -------------------------
# 3. MCP (WLS - Cressie)
# -------------------------
fit_TMP_cressie_REC <- variofit(
vario = var_TMP,
ini.cov.pars = ini_TMP,
cov.model = "spherical",
weights = "cressie",
fix.nugget = TRUE,
nugget = nug_TMP
)variofit: covariance model used is spherical
variofit: weights used: cressie
variofit: minimisation function used: optim # -------------------------
# 4. ML (Maximum Likelihood)
# -------------------------
fit_TMP_ML_REC <- likfit(
geodata = geo_TMP,
coords = geo_TMP$coords,
data = geo_TMP$data,
ini.cov.pars = ini_TMP,
fix.nugget = TRUE,
nugget = nug_TMP,
cov.model = "spherical",
lik.method = "ML"
)kappa not used for the spherical correlation function
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
arguments for the maximisation function.
For further details see documentation for optim.
likfit: It is highly advisable to run this function several
times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.# -------------------------
# 5. REML
# -------------------------
fit_TMP_REML_REC <- likfit(
geodata = geo_TMP,
coords = geo_TMP$coords,
data = geo_TMP$data,
ini.cov.pars = ini_TMP,
fix.nugget = TRUE,
nugget = nug_TMP,
cov.model = "spherical",
lik.method = "REML"
)kappa not used for the spherical correlation function
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
arguments for the maximisation function.
For further details see documentation for optim.
likfit: It is highly advisable to run this function several
times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.# --- Parámetros iniciales RH (del eyefit) ---
ini_RH <- c(sigmasq = 14.7, phi = 11500)
nug_RH <- 1.85
# -------------------------
# 1. MCO (OLS)
# -------------------------
fit_RH_MCO_REC <- variofit(
vario = var_RH,
ini.cov.pars = ini_RH,
cov.model = "exponential",
weights = "equal",
fix.nugget = TRUE,
nugget = nug_RH
)variofit: covariance model used is exponential
variofit: weights used: equal
variofit: minimisation function used: optim # -------------------------
# 2. MCP (npairs)
# -------------------------
fit_RH_npairs_REC <- variofit(
vario = var_RH,
ini.cov.pars = ini_RH,
cov.model = "exponential",
weights = "npairs",
fix.nugget = TRUE,
nugget = nug_RH
)variofit: covariance model used is exponential
variofit: weights used: npairs
variofit: minimisation function used: optim # -------------------------
# 3. MCP Cressie
# -------------------------
fit_RH_cressie_REC <- variofit(
vario = var_RH,
ini.cov.pars = ini_RH,
cov.model = "exponential",
weights = "cressie",
fix.nugget = TRUE,
nugget = nug_RH
)variofit: covariance model used is exponential
variofit: weights used: cressie
variofit: minimisation function used: optim # -------------------------
# 4. ML
# -------------------------
fit_RH_ML_REC <- likfit(
geodata = geo_RH,
coords = geo_RH$coords,
data = geo_RH$data,
ini.cov.pars = ini_RH,
fix.nugget = TRUE,
nugget = nug_RH,
cov.model = "exponential",
lik.method = "ML"
)kappa not used for the exponential correlation function
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
arguments for the maximisation function.
For further details see documentation for optim.
likfit: It is highly advisable to run this function several
times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.# -------------------------
# 5. REML
# -------------------------
fit_RH_REML_REC <- likfit(
geodata = geo_RH,
coords = geo_RH$coords,
data = geo_RH$data,
ini.cov.pars = ini_RH,
fix.nugget = TRUE,
nugget = nug_RH,
cov.model = "exponential",
lik.method = "REML"
)kappa not used for the exponential correlation function
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
arguments for the maximisation function.
For further details see documentation for optim.
likfit: It is highly advisable to run this function several
times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.plot(var_TMP$u, var_TMP$v,
pch = 16,
cex = 1,
xlab = "Distancia (h)",
ylab = "Semivarianza",
main = "TMP — Ajustes del semivariograma (recortado)")
lines(fit_TMP_MCO_REC, col = 1, lwd = 2)
lines(fit_TMP_npairs_REC, col = 2, lwd = 2)
lines(fit_TMP_cressie_REC, col = 3, lwd = 2)
lines(fit_TMP_ML_REC, col = 4, lwd = 2)
lines(fit_TMP_REML_REC, col = 5, lwd = 2)
legend(
"bottomright",
legend = c("MCO", "MCP npairs", "MCP Cressie", "ML", "REML"),
col = 1:5,
lwd = 2
)
plot(var_RH$u, var_RH$v,
pch = 16,
cex = 1,
xlab = "Distancia (h)",
ylab = "Semivarianza",
main = "RH — Ajustes del semivariograma (recortado)")
lines(fit_RH_MCO_REC, col = 1, lwd = 2)
lines(fit_RH_npairs_REC, col = 2, lwd = 2)
lines(fit_RH_cressie_REC, col = 3, lwd = 2)
lines(fit_RH_ML_REC, col = 4, lwd = 2)
lines(fit_RH_REML_REC, col = 5, lwd = 2)
legend(
"bottomright",
legend = c("MCO", "MCP npairs", "MCP Cressie", "ML", "REML"),
col = 1:5,
lwd = 2
)
La estructura espacial de las variables Temperatura (TMP) y Humedad Relativa (RH) fue analizada mediante la estimación de semivariogramas empíricos y sus correspondientes modelos teóricos. En ambos casos, los resultados muestran que la variabilidad espacial presente en los datos no se ajusta de manera clara a un modelo paramétrico estándar, aun después de aplicar procedimientos robustos, centrar los residuales y recortar las distancias máximas consideradas.
Para cada variable se ajustaron modelos de semivariograma mediante MCO, MCP ponderado por número de pares, MCP con pesos de Cressie, máxima verosimilitud (ML) y máxima verosimilitud restringida (REML). En el caso de la temperatura (TMP), los modelos ML y REML producen curvas muy similares, que reproducen adecuadamente el nivel de semivarianza a distancias intermedias y un rango coherente con la separación máxima entre estaciones. Además, REML presenta el mejor compromiso entre verosimilitud y número de parámetros (AIC/BIC ligeramente menores) y es menos sesgado en la estimación de la varianza en muestras pequeñas. Por estas razones, el modelo esférico estimado por REML se adopta como modelo teórico final para TMP y se utiliza en la etapa de kriging. Para la humedad relativa (RH), el variograma empírico es más ruidoso; sin embargo, los modelos exponenciales ajustados por MCP con pesos de Cressie y por REML son los que mejor capturan la tendencia creciente inicial. Dado que REML proporciona también una mejor estimación de la varianza del proceso, se selecciona el modelo exponencial estimado por REML como base para las predicciones de kriging.
1.5 Mapas de predicción de kriging acompañados de cada uno de los mapas respectivos de varianza del error de predicción.
# ===============================================================
# Convertir modelos ML a gstat::vgm()
# ===============================================================
# --- TMP ---
vg_TMP_model <- vgm(
psill = fit_TMP_ML$cov.pars[1],
model = "Sph",
nugget = fit_TMP_ML$nugget,
range = fit_TMP_ML$cov.pars[2]
)
# --- RH ---
# Extraemos parámetros
sig_RH <- fit_RH_ML$cov.pars[1]
phi_RH <- fit_RH_ML$cov.pars[2]
nug_RH <- fit_RH_ML$nugget
# Si ML entrega phi_RH = 0 (caso frecuente)
if (phi_RH <= 0) {
message("ML dio rango=0, se fuerza a un mínimo razonable")
phi_RH <- max(2000, 0.15 * Dmax_rh) # regla razonable
}
vg_RH_model <- vgm(
psill = sig_RH,
model = "Exp",
nugget = nug_RH,
range = phi_RH
)
df_tmp_sp <- df_tmp
coordinates(df_tmp_sp) <- ~ x + y
proj4string(df_tmp_sp) <- CRS(paste0("EPSG:",epsg_crs))
df_rh_sp <- df_rh
coordinates(df_rh_sp) <- ~ x + y
proj4string(df_rh_sp) <- CRS(paste0("EPSG:",epsg_crs))
buffer <- 2000
grd <- expand.grid(
x = seq(min(df_tmp$x)-buffer, max(df_tmp$x)+buffer, length.out=120),
y = seq(min(df_tmp$y)-buffer, max(df_tmp$y)+buffer, length.out=120)
)
coordinates(grd) <- ~ x + y
gridded(grd) <- TRUE
proj4string(grd) <- CRS(paste0("EPSG:",epsg_crs))
# ===============================================================
# Kriging ordinario
# ===============================================================
g_TMP <- gstat(id="resid", formula = resid ~ 1,
data = df_tmp_sp, model = vg_TMP_model)
krig_TMP <- predict(g_TMP, grd)[using ordinary kriging]g_RH <- gstat(id="resid", formula = resid ~ 1,
data = df_rh_sp, model = vg_RH_model)
krig_RH <- predict(g_RH, grd)[using ordinary kriging]# TMP no tiene tendencia → predicción = residual krigeado
TMP_pred_final <- krig_TMP$resid.pred
TMP_var_final <- krig_TMP$resid.var
# RH sí tiene tendencia → reconstrucción media + residuales
beta_rh <- coef(m1_rh)
RH_mean_grid <- beta_rh[1] + beta_rh[2]*grd$x + beta_rh[3]*grd$y
RH_pred_final <- RH_mean_grid + krig_RH$resid.pred
RH_var_final <- krig_RH$resid.var
# Añadir al grid
grd$TMP_pred <- TMP_pred_final
grd$RH_pred <- RH_pred_final
grd$TMP_var <- TMP_var_final
grd$RH_var <- RH_var_final
library(raster)
Adjuntando el paquete: 'raster'The following object is masked from 'package:dplyr':
select# ===============================================================
# Raster en UTM (predicción y varianza)
# ===============================================================
# Raster de predicción en UTM
ras_TMP_utm <- raster(grd["TMP_pred"])
ras_RH_utm <- raster(grd["RH_pred"])
# Raster de varianza en UTM
ras_TMPvar_utm <- raster(grd["TMP_var"])
ras_RHvar_utm <- raster(grd["RH_var"])
# --- Proyección a lat/long ---
crs_ll <- CRS("+proj=longlat +datum=WGS84 +no_defs")
ras_TMP_ll <- projectRaster(ras_TMP_utm, crs = crs_ll)
ras_RH_ll <- projectRaster(ras_RH_utm, crs = crs_ll)
ras_TMPvar_ll <- projectRaster(ras_TMPvar_utm, crs = crs_ll)
ras_RHvar_ll <- projectRaster(ras_RHvar_utm, crs = crs_ll)
# ===============================================================
# Preparar estaciones
# ===============================================================
pts_sf <- st_as_sf(df_tmp, coords = c("x","y"), crs = epsg_crs)
pts_ll <- st_transform(pts_sf, 4326)
# ===============================================================
# Paletas de color para cada capa
# ===============================================================
# Valores de TMP
vals_TMP <- values(ras_TMP_ll)
vals_TMPvar <- values(ras_TMPvar_ll)
# Valores de RH
vals_RH <- values(ras_RH_ll)
vals_RHvar <- values(ras_RHvar_ll)
pal_TMP <- colorNumeric(viridis(200), range(vals_TMP, na.rm = TRUE))
pal_TMPvar <- colorNumeric(viridis(200), range(vals_TMPvar, na.rm = TRUE))
pal_RH <- colorNumeric(viridis(200), range(vals_RH, na.rm = TRUE))
pal_RHvar <- colorNumeric(viridis(200), range(vals_RHvar, na.rm = TRUE))
# Extensión para enfocar el mapa
ext <- extent(ras_TMP_ll)
# ===============================================================
# Leaflet con 4 capas (predicción + varianza)
# ===============================================================
leaflet() %>%
addProviderTiles("CartoDB.Positron") %>%
# --- TMP Predicción ---
addRasterImage(
ras_TMP_ll,
colors = pal_TMP,
opacity = 0.7,
group = "TMP — Predicción"
) %>%
# --- TMP Varianza ---
addRasterImage(
ras_TMPvar_ll,
colors = pal_TMPvar,
opacity = 0.7,
group = "TMP — Varianza"
) %>%
# --- RH Predicción ---
addRasterImage(
ras_RH_ll,
colors = pal_RH,
opacity = 0.7,
group = "RH — Predicción"
) %>%
# --- RH Varianza ---
addRasterImage(
ras_RHvar_ll,
colors = pal_RHvar,
opacity = 0.7,
group = "RH — Varianza"
) %>%
# --- Estaciones ---
addCircleMarkers(
data = pts_ll,
lng = ~st_coordinates(geometry)[,1],
lat = ~st_coordinates(geometry)[,2],
radius = 5,
color = "red",
fillOpacity = 0.7
) %>%
# --- Control de capas ---
addLayersControl(
overlayGroups = c(
"TMP — Predicción", "TMP — Varianza",
"RH — Predicción", "RH — Varianza"
),
options = layersControlOptions(collapsed = FALSE)
) %>%
# --- Ajustar vista ---
fitBounds(
lng1 = xmin(ext), lat1 = ymin(ext),
lng2 = xmax(ext), lat2 = ymax(ext)
)2a. Modelos de semivariograma alternativos mediante regresión lineal y no lineal
En esta sección mostramos ejemplos de ajuste de modelos de semivariograma que no están predefinidos en los paquetes geoR y gstat, usando herramientas generales de regresión (lm, nls) y esquemas de mínimos cuadrados ponderados iterativos (IRLS). Trabajamos sobre datos simulados.
library(dplyr)
library(ggplot2)
library(tidyr)
library(sp)
library(knitr)
set.seed(123)2.1a. Simulación de datos espaciales
Creamos un campo espacial sintético con dependencia espacial suave.
n <- 80
coords <- cbind(runif(n, 0, 10), runif(n, 0, 10))
colnames(coords) <- c("x", "y")
# Proceso con estructura espacial (tendencia + ruido correlacionado)
z_true <- 3 + 2*sin(coords[,1]/2) + 1.5*cos(coords[,2]/3)
z <- z_true + rnorm(n, 0, 0.4)
data <- data.frame(coords, z)2.2a. Cálculo empírico del semivariograma
# Matriz de distancias y semivarianzas
dist_mat <- as.matrix(dist(coords))
gamma_mat <- 0.5 * (outer(z, z, "-"))^2
# Agrupamos distancias en bins
breaks <- seq(0, max(dist_mat), length.out = 15)
hcut <- cut(dist_mat, breaks)
emp_variog <- data.frame(
h = breaks[-length(breaks)],
gamma = tapply(gamma_mat, hcut, mean, na.rm = TRUE),
np = tapply(gamma_mat, hcut, length)
)
emp_variog <- emp_variog %>% filter(!is.na(gamma))
ggplot(emp_variog, aes(h, gamma)) +
geom_point() +
theme_minimal() +
labs(title = "Semivariograma empírico",
x = "Distancia h", y = "Semivarianza γ(h)")
2.3a. Ajuste de modelos personalizados (fuera de geoR/gstat)
Los paquetes geoR y gstat ofrecen modelos clásicos (esférico, exponencial, gaussiano, Matérn, etc.). Aquí usamos modelos no estándar:
- Modelo polinomial cuadrático decreciente
- Modelo logístico modificado
Ambos estimados por MCO y Mínimos Cuadrados No Lineales (NLS).
# Modelo polinomial cuadrático (MCO)
poly_model <- lm(gamma ~ poly(h, 2, raw = TRUE), data = emp_variog)
summary(poly_model)
Call:
lm(formula = gamma ~ poly(h, 2, raw = TRUE), data = emp_variog)
Residuals:
Min 1Q Median 3Q Max
-2.24362 -0.14059 0.03095 0.25494 1.50051
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.01771 0.59689 -0.030 0.9769
poly(h, 2, raw = TRUE)1 0.51382 0.23253 2.210 0.0492 *
poly(h, 2, raw = TRUE)2 0.05062 0.01881 2.691 0.0210 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8533 on 11 degrees of freedom
Multiple R-squared: 0.9682, Adjusted R-squared: 0.9624
F-statistic: 167.6 on 2 and 11 DF, p-value: 5.78e-09# Modelo logístico modificado para el semivariograma
logistic_fun <- function(h, nug, sill, range){
nug + (sill - nug)/(1 + exp(-0.5 * (h - range)))
}
nls_model <- nls(
gamma ~ logistic_fun(h, nug, sill, range),
data = emp_variog,
start = list(nug = 0.1, sill = 1, range = 5)
)
summary(nls_model)
Formula: gamma ~ logistic_fun(h, nug, sill, range)
Parameters:
Estimate Std. Error t value Pr(>|t|)
nug 0.4017 0.6340 0.634 0.539
sill 13.6624 1.1503 11.877 1.29e-07 ***
range 7.5432 0.6718 11.228 2.30e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.108 on 11 degrees of freedom
Number of iterations to convergence: 13
Achieved convergence tolerance: 3.423e-062.4a. Ponderaciones fijas y ponderaciones iterativas (WLS)
# Ponderaciones fijas: inversas al número de pares (np)
w_fixed <- 1 / emp_variog$np
# MCP (WLS con pesos fijos)
wls_model <- lm(
gamma ~ poly(h, 2, raw = TRUE),
data = emp_variog,
weights = w_fixed
)
summary(wls_model)
Call:
lm(formula = gamma ~ poly(h, 2, raw = TRUE), data = emp_variog,
weights = w_fixed)
Weighted Residuals:
Min 1Q Median 3Q Max
-0.41812 0.01302 0.04633 0.05888 0.13184
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.71255 1.42572 0.500 0.6271
poly(h, 2, raw = TRUE)1 -0.27545 0.53071 -0.519 0.6140
poly(h, 2, raw = TRUE)2 0.11963 0.03914 3.057 0.0109 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1418 on 11 degrees of freedom
Multiple R-squared: 0.9327, Adjusted R-squared: 0.9205
F-statistic: 76.27 on 2 and 11 DF, p-value: 3.569e-07# Iterativo: pesos según residuos del modelo anterior
iter_model <- poly_model
for(i in 1:5){
res <- residuals(iter_model)
w_iter <- 1 / (abs(res) + 1e-4) # evitar división por cero
iter_model <- lm(
gamma ~ poly(h, 2, raw = TRUE),
data = emp_variog,
weights = w_iter
)
}
summary(iter_model)
Call:
lm(formula = gamma ~ poly(h, 2, raw = TRUE), data = emp_variog,
weights = w_iter)
Weighted Residuals:
Min 1Q Median 3Q Max
-1.60782 -0.31565 0.00508 0.48586 1.02871
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.046210 0.086546 -0.534 0.604
poly(h, 2, raw = TRUE)1 0.479943 0.043616 11.004 2.82e-07 ***
poly(h, 2, raw = TRUE)2 0.056549 0.004696 12.042 1.12e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7239 on 11 degrees of freedom
Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998
F-statistic: 3.663e+04 on 2 and 11 DF, p-value: < 2.2e-162.5a. ML y Composite Likelihood (CL) con modelo exponencial casero
Definimos un modelo exponencial de covarianza y el semivariograma correspondiente y estimamos sus parámetros por:
- ML (verosimilitud gaussiana completa).
- Composite Likelihood (CL) pareada.
# Centramos los datos para suponer media cero en ML
z_c <- z - mean(z)
# Semivariograma exponencial
gamma_exp <- function(h, sigma2, rho, tau2){
# nugget (tau2) + parte estructural
tau2 + sigma2 * (1 - exp(-h / rho))
}
# Matriz de covarianza exponencial a partir de parámetros
make_cov_exp <- function(theta, dist_mat){
# theta en escala log: (log(sigma2), log(rho), log(tau2))
sigma2 <- exp(theta[1])
rho <- exp(theta[2])
tau2 <- exp(theta[3])
C <- sigma2 * exp(-dist_mat / rho)
diag(C) <- diag(C) + tau2 + 1e-6 # término pequeño para estabilidad
C
}
# Log-verosimilitud negativa (ML completo)
negloglik_full <- function(theta, z, dist_mat){
C <- make_cov_exp(theta, dist_mat)
cholC <- chol(C)
logdet <- 2 * sum(log(diag(cholC)))
sol <- backsolve(cholC, z, transpose = TRUE)
quad <- crossprod(sol)
0.5 * (length(z) * log(2 * pi) + logdet + quad)
}
# Valores iniciales
start_theta <- log(c(var(z_c) * 0.7, max(dist_mat) / 3, var(z_c) * 0.3))
fit_ml <- optim(
start_theta, negloglik_full,
z = z_c, dist_mat = dist_mat,
method = "BFGS"
)
theta_ml <- exp(fit_ml$par)
names(theta_ml) <- c("sigma2", "rho", "tau2")
theta_ml sigma2 rho tau2
6.15342493 29.42527721 0.06726359 Composite Likelihood (pairwise)
# Índices de pares i<j
pairs_idx <- which(upper.tri(dist_mat), arr.ind = TRUE)
d_vec <- dist_mat[pairs_idx]
i_idx <- pairs_idx[,1]
j_idx <- pairs_idx[,2]
negloglik_cl <- function(theta, z, d_vec, i_idx, j_idx){
sigma2 <- exp(theta[1])
rho <- exp(theta[2])
tau2 <- exp(theta[3])
c0 <- sigma2 + tau2 + 1e-6
cij <- sigma2 * exp(-d_vec / rho)
ll <- 0
for(k in seq_along(d_vec)){
S11 <- c0
S22 <- c0
S12 <- cij[k]
detS <- S11 * S22 - S12^2
y1 <- z[i_idx[k]]
y2 <- z[j_idx[k]]
quad <- (S22 * y1^2 - 2 * S12 * y1 * y2 + S11 * y2^2) / detS
ll <- ll - 0.5 * (2 * log(2 * pi) + log(detS) + quad)
}
-ll # negativo para minimizar
}
fit_cl <- optim(
start_theta, negloglik_cl,
z = z_c, d_vec = d_vec,
i_idx = i_idx, j_idx = j_idx,
method = "BFGS",
control = list(maxit = 60)
)
theta_cl <- exp(fit_cl$par)
names(theta_cl) <- c("sigma2", "rho", "tau2")
theta_cl sigma2 rho tau2
3.6887698649 2.8761844713 0.0001174517 Tabla de parámetros ML vs CL
param_df <- rbind(
ML = theta_ml,
CL = theta_cl
)
kable(param_df, digits = 3,
caption = "Parámetros del modelo exponencial estimados por ML y CL")| sigma2 | rho | tau2 | |
|---|---|---|---|
| ML | 6.153 | 29.425 | 0.067 |
| CL | 3.689 | 2.876 | 0.000 |
2.6a. Visualización comparativa de todos los modelos
Incluimos MCO, MCP, WLS iterativo, NLS logístico, ML exponencial y CL exponencial.
h_seq <- seq(min(emp_variog$h), max(emp_variog$h), length.out = 200)
# Predicciones de los modelos polinomiales y NLS
pred_poly <- predict(poly_model, newdata = data.frame(h = h_seq))
pred_wls <- predict(wls_model, newdata = data.frame(h = h_seq))
pred_iter <- predict(iter_model, newdata = data.frame(h = h_seq))
pred_nls <- predict(nls_model, newdata = data.frame(h = h_seq))
# Predicciones usando los parámetros ML y CL del modelo exponencial
pred_ml <- gamma_exp(
h_seq,
sigma2 = theta_ml["sigma2"],
rho = theta_ml["rho"],
tau2 = theta_ml["tau2"]
)
pred_cl <- gamma_exp(
h_seq,
sigma2 = theta_cl["sigma2"],
rho = theta_cl["rho"],
tau2 = theta_cl["tau2"]
)
plot_df <- data.frame(
h = h_seq,
MCO = pred_poly,
MCP = pred_wls,
Iter_WLS = pred_iter,
NLS_logistica = pred_nls,
ML_exp = pred_ml,
CL_exp = pred_cl
) |>
tidyr::pivot_longer(-h, names_to = "metodo", values_to = "gamma")
ggplot() +
geom_point(data = emp_variog,
aes(h, gamma),
alpha = 0.6) +
geom_line(data = plot_df,
aes(h, gamma, color = metodo),
linewidth = 1) +
theme_minimal() +
labs(
title = "Comparación de estimadores del semivariograma",
x = "Distancia h",
y = "Semivarianza γ(h)"
)
2b. Modelos de semivariograma alternativos mediante regresión lineal y no lineal
En esta sección mostramos ejemplos de ajuste de modelos de semivariograma que no están predefinidos en los paquetes geoR y gstat, usando herramientas generales de regresión (lm, nls) y esquemas de mínimos cuadrados ponderados iterativos (IRLS). Trabajamos sobre el semivariograma empírico clásico de la temperatura (TMP), pero el mismo procedimiento se puede replicar para RH.
2.1b. Preparar los datos del semivariograma empírico
# Datos del semivariograma clásico de TMP en formato data.frame
vg_tmp_df <- as.data.frame(vg_TMP) |>
dplyr::transmute(
h = dist, # distancia media del bin
gamma = gamma, # semivarianza empírica
np = np # número de pares
)
head(vg_tmp_df) h gamma np
1 5632.817 0.1900000 3
2 7559.192 0.2108333 6
3 9642.049 0.3266667 6
4 12306.191 0.1842857 7
5 14861.818 0.1154167 12
6 16571.822 0.1333333 32.2b. Modelo polinomial (regresión lineal)
Como ejemplo simple, ajustamos un modelo polinomial cuadrático: \(\gamma(h) \approx \beta_0 + \beta_1 h + \beta_2 h^2\).
# Modelo polinomial: gamma(h) ≈ β0 + β1 h + β2 h^2
vg_poly_lm <- lm(
gamma ~ h + I(h^2),
data = vg_tmp_df,
weights = np # pesos fijos: número de pares en cada bin
)
summary(vg_poly_lm)
Call:
lm(formula = gamma ~ h + I(h^2), data = vg_tmp_df, weights = np)
Weighted Residuals:
Min 1Q Median 3Q Max
-0.25967 -0.09674 -0.02595 0.04185 0.36062
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.979e-01 1.376e-01 2.165 0.0556 .
h -1.310e-05 1.565e-05 -0.837 0.4220
I(h^2) 3.947e-10 4.003e-10 0.986 0.3474
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2048 on 10 degrees of freedom
Multiple R-squared: 0.128, Adjusted R-squared: -0.04634
F-statistic: 0.7343 on 2 and 10 DF, p-value: 0.504Este modelo no garantiza que el semivariograma resultante sea definido positivo, pero ilustra cómo usar regresión lineal estándar para ajustar una forma funcional arbitraria sobre los puntos empíricos.
2.3b. Modelo racional (no lineal) con pesos fijos
Ahora proponemos un modelo racional no estándar:
\(\gamma(h) = \tau^2 + \sigma^2 \frac{h^2}{h^2 + \phi^2}\)
que no está incluido explícitamente en geoR/gstat. Lo ajustamos con nls y pesos fijos.
# Modelo racional de semivariograma
gamma_rat <- function(h, nugget, sill, range) {
nugget + sill * (h^2 / (h^2 + range^2))
}
# Valores iniciales razonables
start_vals <- list(
nugget = min(vg_tmp_df$gamma, na.rm = TRUE),
sill = max(vg_tmp_df$gamma, na.rm = TRUE) - min(vg_tmp_df$gamma, na.rm = TRUE),
range = max(vg_tmp_df$h, na.rm = TRUE) / 3
)
# Ajuste no lineal con pesos fijos (np)
vg_rat_nls <- try(
nls(
gamma ~ gamma_rat(h, nugget, sill, range),
data = vg_tmp_df,
start = start_vals,
weights = np, # pesos fijos
control = nls.control(maxiter = 200, warnOnly = TRUE)
),
silent = TRUE
)Warning in nls(gamma ~ gamma_rat(h, nugget, sill, range), data = vg_tmp_df, :
gradiente singularif (inherits(vg_rat_nls, "try-error")) {
warning("El ajuste nls del modelo racional no convergió; se usan los valores iniciales como aproximación.")
par_nls <- unlist(start_vals)
} else {
# Intentar obtener los coeficientes sin llamar a summary() (que puede fallar
# por problemas numéricos en la matriz de información)
coef_try <- try(coef(vg_rat_nls), silent = TRUE)
if (inherits(coef_try, "try-error")) {
warning("El ajuste nls produjo un objeto singular; se usan los valores iniciales como aproximación.")
par_nls <- unlist(start_vals)
} else {
par_nls <- coef_try
}
}2.4b. Mínimos cuadrados ponderados iterativos (IRLS) con pesos que dependen de los parámetros
Para ilustrar el uso de matrices de ponderación que cambian en cada iteración, implementamos un esquema IRLS donde, en la iteración (k), hacemos:
\(w_i^{(k)} = \frac{n_i}{\bigl(\hat\gamma_i^{(k)}\bigr)^2}, \quad \hat\gamma_i^{(k)} = \hat \gamma(h_i; \theta^{(k)})\), y minimizamos \((\sum_i w_i^{(k)} (\gamma_i - \hat\gamma_i)^2)\) respecto a \(\theta\).
irls_vg_rat <- function(vg_df, par_init, max_iter = 15, tol = 1e-4) {
par_cur <- par_init
obj_vals <- numeric(max_iter)
for (k in seq_len(max_iter)) {
# semivariograma teórico con parámetros actuales
g_cur <- gamma_rat(vg_df$h, par_cur[1], par_cur[2], par_cur[3])
# pesos dependientes del modelo actual (evitar dividir por cero)
w <- vg_df$np / (g_cur^2 + 1e-6)
# función objetivo de Mínimos Cuadrados Ponderados
obj <- function(par) {
g_mod <- gamma_rat(vg_df$h, par[1], par[2], par[3])
sum(w * (vg_df$gamma - g_mod)^2)
}
opt <- optim(
par = par_cur,
fn = obj,
method = "L-BFGS-B",
lower = c(0, 0, 1e-6) # asegurar positividad básica
)
obj_vals[k] <- opt$value
# criterio de parada
if (max(abs(opt$par - par_cur)) < tol) {
par_cur <- opt$par
obj_vals <- obj_vals[1:k]
break
}
par_cur <- opt$par
}
list(par = par_cur, obj = obj_vals, iters = length(obj_vals))
}
# Valores iniciales tomados del ajuste nls
par0_rat <- c(
nugget = coef(vg_rat_nls)[["nugget"]],
sill = coef(vg_rat_nls)[["sill"]],
range = coef(vg_rat_nls)[["range"]]
)
vg_rat_irls <- irls_vg_rat(vg_tmp_df, par0_rat)
vg_rat_irls$par # parámetros finales IRLS nugget sill range
0.211105838 0.004014643 0.000001000 vg_rat_irls$obj # valores de la función objetivo por iteración[1] 11.24585 10.39688vg_rat_irls$iters[1] 22.5b. Comparación visual de modelos alternativos
# Grid de distancias para comparar curvas teóricas
h_grid <- seq(0, max(vg_tmp_df$h, na.rm = TRUE), length.out = 100)
# Curva polinomial (a partir de los coeficientes de lm)
poly_coef <- coef(vg_poly_lm)
gamma_poly <- poly_coef[1] + poly_coef[2] * h_grid + poly_coef[3] * h_grid^2
# Curva racional (nls)
gamma_nls <- gamma_rat(h_grid, par_nls["nugget"], par_nls["sill"], par_nls["range"])
# Curva racional (IRLS)
par_irls <- vg_rat_irls$par
gamma_irls <- gamma_rat(h_grid, par_irls[1], par_irls[2], par_irls[3])
plot(vg_tmp_df$h, vg_tmp_df$gamma,
xlab = "h (distancia)",
ylab = "semivarianza",
main = "TMP — Modelos alternativos de semivariograma")
lines(h_grid, gamma_poly, lwd = 2)
lines(h_grid, gamma_nls, lwd = 2, lty = 2)
lines(h_grid, gamma_irls, lwd = 2, lty = 3)
legend(
"bottomright",
legend = c("Empírico", "Polinomial LM", "Racional NLS", "Racional IRLS"),
lwd = c(NA, 2, 2, 2),
pch = c(1, NA, NA, NA),
lty = c(NA, 1, 2, 3)
)
3. Funciones objetivo generales para MCO, MCP, MV y CL y su programación con optim
En esta sección definimos de forma explícita las funciones objetivo asociadas a:
Mínimos Cuadrados Ordinarios (MCO)
Mínimos Cuadrados Ponderados (MCP) sobre el semivariograma empírico
Máxima Verosimilitud (MV) sobre los datos residuales
Composite Likelihood (CL) basada en pares
Mínimos Cuadrados Ordinarios (MCO)
Sea
\(\hat{\boldsymbol\gamma} = (\hat\gamma(h_1),\dots,\hat\gamma(h_K))^\top\)
el vector de semivariogramas empíricos y
\(\boldsymbol\gamma(\theta) = (\gamma(h_1;\theta),\dots,\gamma(h_K;\theta))^\top\)
el semivariograma teórico para el vector de parámetros \(\theta\).La función objetivo de MCO es
\[ Q_{\text{MCO}}(\theta) = \bigl\|\hat{\boldsymbol\gamma} - \boldsymbol\gamma(\theta)\bigr\|^2 = \sum_{k=1}^K \left[\hat\gamma(h_k) - \gamma(h_k;\theta)\right]^2. \]
Mínimos Cuadrados Ponderados (MCP) sobre el semivariograma empírico
Introducimos una matriz de ponderaciones fija \(W\) (por ejemplo diagonal, con pesos \(w_k\)). La función objetivo general es
\[ Q_{\text{MCP}}(\theta) =\bigl(\hat{\boldsymbol\gamma} - \boldsymbol\gamma(\theta)\bigr)^\top W\, \bigl(\hat{\boldsymbol\gamma} - \boldsymbol\gamma(\theta)\bigr) = \sum_{k=1}^K w_k \left[\hat\gamma(h_k) - \gamma(h_k;\theta)\right]^2. \]
Máxima Verosimilitud (MV) sobre los datos (o residuales)
Sea \(\mathbf{r}=(r_1,\dots,r_n)^\top\) el vector de datos (o residuales) y \(C(\theta)\) la matriz de covarianza gaussiana asociada al modelo espacial. La log-verosimilitud (hasta una constante) viene dada por
\[ \ell(\theta) = -\tfrac12\left\{ \log\bigl|C(\theta)\bigr| + \mathbf{r}^\top C(\theta)^{-1}\mathbf{r} \right\}. \]
Equivalentemente, la función objetivo a minimizar en
optimpuede definirse como\[ Q_{\text{MV}}(\theta) = \tfrac12\left\{ \log\bigl|C(\theta)\bigr| + \mathbf{r}^\top C(\theta)^{-1}\mathbf{r} \right\}. \]
Composite Likelihood (CL) basada en pares
Para cada par \((i,j)\) con \(i<j\), sea \(\mathbf{r}_{ij} = (r_i,r_j)^\top\) y \(\Sigma_{ij}(\theta)\) la submatriz \(2\times 2\) de \(C(\theta)\) correspondiente a ese par.
La log-likelihood pareada gaussiana (hasta constante) es\[ \ell_{\text{CL}}(\theta) = -\tfrac12 \sum_{i<j} \left\{ \log\bigl|\Sigma_{ij}(\theta)\bigr| + \mathbf{r}_{ij}^\top \Sigma_{ij}(\theta)^{-1} \mathbf{r}_{ij} \right\}. \]
Por tanto, la función objetivo de CL a minimizar es
\[ Q_{\text{CL}}(\theta) = \tfrac12 \sum_{i<j} \left\{ \log\bigl|\Sigma_{ij}(\theta)\bigr| + \mathbf{r}_{ij}^\top \Sigma_{ij}(\theta)^{-1} \mathbf{r}_{ij} \right\}. \]
En todos los casos, el estimador se obtiene como
\(\hat\theta = \arg\min_\theta Q(\theta)\) usando optim().
Usaremos un modelo exponencial estándar como ejemplo, pero las funciones están escritas de forma genérica y se pueden adaptar a otros modelos.
3.1 Modelo exponencial: semivariograma y covarianza
# Semivariograma exponencial
gamma_exp <- function(h, nugget, sill, range) {
nugget + sill * (1 - exp(-h / range))
}
# Covarianza exponencial (parte estructural)
cov_exp <- function(h, sill, range) {
sill * exp(-h / range)
}3.2 Funciones objetivo MCO y MCP sobre el semivariograma empírico
Trabajamos con los puntos del semivariograma empírico de TMP:
h <- vg_tmp_df$h
gamma_hat <- vg_tmp_df$gamma
# Ejemplo de pesos fijos para MCP (se pueden cambiar)
w_fix <- vg_tmp_df$np / (h^2 + 1e-6)Para garantizar positividad en los parámetros, optimizamos en la escala log: (= (^2, ^2, )).
obj_mco <- function(log_par, h, gamma_hat) {
nugget <- exp(log_par[1])
sill <- exp(log_par[2])
range <- exp(log_par[3])
g_mod <- gamma_exp(h, nugget, sill, range)
mean((gamma_hat - g_mod)^2)
}
obj_mcp <- function(log_par, h, gamma_hat, w) {
nugget <- exp(log_par[1])
sill <- exp(log_par[2])
range <- exp(log_par[3])
g_mod <- gamma_exp(h, nugget, sill, range)
sum(w * (gamma_hat - g_mod)^2)
}3.3 Ajuste numérico de MCO y MCP con optim
log_par_init <- log(c(
nugget = 0.05,
sill = 0.20,
range = max(h, na.rm = TRUE) / 3
))
# MCO
fit_mco <- optim(
par = log_par_init,
fn = obj_mco,
h = h,
gamma_hat = gamma_hat,
method = "L-BFGS-B"
)
# MCP con pesos fijos
fit_mcp <- optim(
par = log_par_init,
fn = obj_mcp,
h = h,
gamma_hat = gamma_hat,
w = w_fix,
method = "L-BFGS-B"
)
theta_mco <- exp(fit_mco$par)
theta_mcp <- exp(fit_mcp$par)
theta_mco nugget sill range
1.189453e-04 2.150231e-01 1.779064e+03 theta_mcp nugget sill range
5.368268e-02 2.283882e-01 9.650538e+03 3.4 Función objetivo de Máxima Verosimilitud (MV)
Usamos los residuales espaciales de TMP (df_tmp$resid) y sus coordenadas centradas (x, y):
coords_tmp <- as.matrix(df_tmp[, c("x", "y")])
z_tmp <- df_tmp$resid
D_tmp <- as.matrix(dist(coords_tmp))La log-verosimilitud negativa para un campo gaussiano con covarianza exponencial y nugget es:
obj_mv <- function(log_par, D, z) {
nugget <- exp(log_par[1])
sill <- exp(log_par[2])
range <- exp(log_par[3])
n <- length(z)
# Matriz de covarianza
Sigma <- cov_exp(D, sill, range)
diag(Sigma) <- diag(Sigma) + nugget
chol_S <- tryCatch(chol(Sigma), error = function(e) NULL)
if (is.null(chol_S)) return(1e10) # castigo si la matriz no es definida positiva
logdet <- 2 * sum(log(diag(chol_S)))
v <- backsolve(chol_S, z, transpose = TRUE)
quad <- sum(v^2)
0.5 * (logdet + quad) # constante n log(2π) omitida
}fit_mv <- optim(
par = log_par_init,
fn = obj_mv,
D = D_tmp,
z = z_tmp,
method = "L-BFGS-B"
)
theta_mv <- exp(fit_mv$par)
theta_mv nugget sill range
2.091233e-01 6.479977e-02 1.422553e+04 3.5 Composite Likelihood (CL) pareado
Definimos una composite likelihood basada en pares ((i,j)) usando la log-verosimilitud bivariante gaussiana para cada par:
obj_cl <- function(log_par, D, z) {
nugget <- exp(log_par[1])
sill <- exp(log_par[2])
range <- exp(log_par[3])
v <- sill + nugget # varianza marginal
idx <- which(upper.tri(D), arr.ind = TRUE)
neg_cl <- 0
for (k in seq_len(nrow(idx))) {
i <- idx[k, 1]
j <- idx[k, 2]
h_ij <- D[i, j]
c_ij <- sill * exp(-h_ij / range) # covarianza
detS <- v^2 - c_ij^2
if (detS <= 0) return(1e10)
inv00 <- v / detS
inv01 <- -c_ij / detS
zi <- z[i]
zj <- z[j]
quad <- inv00 * zi^2 + 2 * inv01 * zi * zj + inv00 * zj^2
neg_cl <- neg_cl + 0.5 * (log(detS) + quad)
}
neg_cl
}fit_cl <- optim(
par = log_par_init,
fn = obj_cl,
D = D_tmp,
z = z_tmp,
method = "L-BFGS-B"
)
theta_cl <- exp(fit_cl$par)
theta_cl nugget sill range
1.938503e-01 7.675584e-02 8.992130e+03 3.6 Comparación de parámetros entre métodos
Finalmente, recopilamos los parámetros estimados por cada método para compararlos con los obtenidos por geoR (MCO, MCP, ML, REML) y/o gstat:
library(tibble)
comparacion_par <- tibble(
metodo = c("MCO", "MCP", "MV", "CL"),
nugget = c(theta_mco[1], theta_mcp[1], theta_mv[1], theta_cl[1]),
sill = c(theta_mco[2], theta_mcp[2], theta_mv[2], theta_cl[2]),
range = c(theta_mco[3], theta_mcp[3], theta_mv[3], theta_cl[3])
)
comparacion_par# A tibble: 4 × 4
metodo nugget sill range
<chr> <dbl> <dbl> <dbl>
1 MCO 0.000119 0.215 1779.
2 MCP 0.0537 0.228 9651.
3 MV 0.209 0.0648 14226.
4 CL 0.194 0.0768 8992.# Grid de distancias para graficar las curvas teóricas
h_grid <- seq(
from = 0,
to = max(vg_tmp_df$h, na.rm = TRUE),
length.out = 200
)
# Paleta simple para los métodos
cols_metodos <- c(
"MCO" = "black",
"MCP" = "red",
"MV" = "blue",
"CL" = "darkgreen"
)
# Gráfico 1: rango completo de distancias
plot(
vg_tmp_df$h, vg_tmp_df$gamma,
xlab = "Distancia (h)",
ylab = "Semivarianza",
main = "TMP — Semivariograma empírico y modelos MCO/MCP/MV/CL",
pch = 1
)
# Añadir las curvas teóricas de cada método
for (i in seq_len(nrow(comparacion_par))) {
met <- comparacion_par$metodo[i]
nugget <- comparacion_par$nugget[i]
sill <- comparacion_par$sill[i]
range <- comparacion_par$range[i]
lines(
h_grid,
gamma_exp(h_grid, nugget, sill, range),
col = cols_metodos[met],
lwd = 2
)
}
legend(
"bottomright",
legend = comparacion_par$metodo,
col = cols_metodos[comparacion_par$metodo],
lwd = 2,
pch = NA,
bty = "n"
)
# Gráfico 2: zoom en distancias cortas (por ejemplo, 60% del máximo)
h_max_zoom <- 0.6 * max(vg_tmp_df$h, na.rm = TRUE)
plot(
vg_tmp_df$h, vg_tmp_df$gamma,
xlab = "Distancia (h)",
ylab = "Semivarianza",
main = "TMP — Zoom del semivariograma (distancias cortas)",
pch = 1,
xlim = c(0, h_max_zoom),
ylim = c(0, max(vg_tmp_df$gamma, na.rm = TRUE))
)
for (i in seq_len(nrow(comparacion_par))) {
met <- comparacion_par$metodo[i]
nugget <- comparacion_par$nugget[i]
sill <- comparacion_par$sill[i]
range <- comparacion_par$range[i]
lines(
h_grid,
gamma_exp(h_grid, nugget, sill, range),
col = cols_metodos[met],
lwd = 2
)
}
legend(
"bottomright",
legend = comparacion_par$metodo,
col = cols_metodos[comparacion_par$metodo],
lwd = 2,
pch = NA,
bty = "n"
)
3.7 Comparación con geoR y gstat
En esta sección ajustamos el mismo modelo exponencial usando las funciones estándar de geoR y gstat, y comparamos sus parámetros con los obtenidos manualmente por MCO, MCP, MV y CL.
library(geoR)
library(gstat)
library(dplyr)
library(tibble)# Objeto geostatístico para TMP (residuales)
geo_TMP <- as.geodata(
df_tmp |> dplyr::select(x, y, resid),
coords.col = 1:2,
data.col = 3
)
# Semivariograma empírico geoR (usamos el mismo máximo que en vg_tmp_df)
max_h <- max(vg_tmp_df$h, na.rm = TRUE)
vario_TMP_geoR <- variog(
geo_TMP,
max.dist = max_h,
uvec = seq(0, max_h, length = 15)
)variog: computing omnidirectional variogramtheta_mv_row <- comparacion_par |>
filter(metodo == "MV") |>
slice(1)
ini_cov <- c(theta_mv_row$sill, theta_mv_row$range) # (sigmasq, phi)
nug_ini <- theta_mv_row$nugget
ini_cov sill range
6.479977e-02 1.422553e+04 nug_ini nugget
0.2091233 # MCO geoR (pesos iguales)
fit_TMP_MCO_geoR <- variofit(
vario = vario_TMP_geoR,
cov.model = "exponential",
ini.cov.pars = ini_cov,
nugget = nug_ini,
weights = "equal"
)variofit: covariance model used is exponential
variofit: weights used: equal
variofit: minimisation function used: optim # MCP geoR (npairs)
fit_TMP_MCP_np_geoR <- variofit(
vario = vario_TMP_geoR,
cov.model = "exponential",
ini.cov.pars = ini_cov,
nugget = nug_ini,
weights = "npairs"
)variofit: covariance model used is exponential
variofit: weights used: npairs
variofit: minimisation function used: optim # MCP geoR (Cressie–Hawkins)
fit_TMP_MCP_cr_geoR <- variofit(
vario = vario_TMP_geoR,
cov.model = "exponential",
ini.cov.pars = ini_cov,
nugget = nug_ini,
weights = "cressie"
)variofit: covariance model used is exponential
variofit: weights used: cressie
variofit: minimisation function used: optim # ML geoR
fit_TMP_ML_geoR <- likfit(
geo_TMP,
cov.model = "exponential",
ini.cov.pars = ini_cov,
nugget = nug_ini,
lik.method = "ML"
)kappa not used for the exponential correlation function
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
arguments for the maximisation function.
For further details see documentation for optim.
likfit: It is highly advisable to run this function several
times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.# REML geoR
fit_TMP_REML_geoR <- likfit(
geo_TMP,
cov.model = "exponential",
ini.cov.pars = ini_cov,
nugget = nug_ini,
lik.method = "REML"
)kappa not used for the exponential correlation function
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
arguments for the maximisation function.
For further details see documentation for optim.
likfit: It is highly advisable to run this function several
times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.# Resumen de parámetros de geoR
geoR_par <- tibble(
metodo = c("geoR_MCO", "geoR_MCP_np", "geoR_MCP_Cressie", "geoR_ML", "geoR_REML"),
nugget = c(
fit_TMP_MCO_geoR$nugget,
fit_TMP_MCP_np_geoR$nugget,
fit_TMP_MCP_cr_geoR$nugget,
fit_TMP_ML_geoR$nugget,
fit_TMP_REML_geoR$nugget
),
sill = c(
fit_TMP_MCO_geoR$cov.pars[1],
fit_TMP_MCP_np_geoR$cov.pars[1],
fit_TMP_MCP_cr_geoR$cov.pars[1],
fit_TMP_ML_geoR$sigmasq,
fit_TMP_REML_geoR$sigmasq
),
range = c(
fit_TMP_MCO_geoR$cov.pars[2],
fit_TMP_MCP_np_geoR$cov.pars[2],
fit_TMP_MCP_cr_geoR$cov.pars[2],
fit_TMP_ML_geoR$phi,
fit_TMP_REML_geoR$phi
)
)
geoR_par# A tibble: 5 × 4
metodo nugget sill range
<chr> <dbl> <dbl> <dbl>
1 geoR_MCO 0.0770 0.215 14226.
2 geoR_MCP_np 0.117 0.145 14226.
3 geoR_MCP_Cressie 0.180 0.102 14226.
4 geoR_ML 0.181 0.0962 14226.
5 geoR_REML 0.152 0.165 14226.3.7.2 Ajustes con gstat::fit.variogram
Usamos el semivariograma empírico de gstat (vg_TMP) calculado en la sección 1.2.
Tomamos como modelo inicial los parámetros de MV:
# Modelo inicial para gstat (Exponencial)
vgm_ini <- vgm(
psill = theta_mv_row$sill,
model = "Exp",
range = theta_mv_row$range,
nugget = theta_mv_row$nugget
)
vgm_ini model psill range
1 Nug 0.20912331 0.00
2 Exp 0.06479977 14225.53# 0 = OLS, 1 = WLS (por defecto), 2 = Cressie
fit_gstat_OLS <- fit.variogram(
vg_TMP,
model = vgm_ini,
fit.method = 0
)
fit_gstat_WLS <- fit.variogram(
vg_TMP,
model = vgm_ini,
fit.method = 1
)Warning in fit.variogram(vg_TMP, model = vgm_ini, fit.method = 1): No
convergence after 200 iterations: try different initial values?fit_gstat_Cress <- fit.variogram(
vg_TMP,
model = vgm_ini,
fit.method = 2
)Warning in fit.variogram(vg_TMP, model = vgm_ini, fit.method = 2): No
convergence after 200 iterations: try different initial values?fit_gstat_OLS model psill range
1 Nug 0.20912331 0.00
2 Exp 0.06479977 14225.53fit_gstat_WLS model psill range
1 Nug 0.2224955 0.00
2 Exp 0.1329826 76850.73fit_gstat_Cress model psill range
1 Nug 0.2224298 0.00
2 Exp 0.1343491 77812.69Extraemos los parámetros (recordar: la primera fila del vgm es el nugget, la segunda el psill de la estructura exponencial):
gstat_par <- tibble(
metodo = c("gstat_OLS", "gstat_WLS_np", "gstat_WLS_Cressie"),
nugget = c(
fit_gstat_OLS$psill[1],
fit_gstat_WLS$psill[1],
fit_gstat_Cress$psill[1]
),
sill = c(
fit_gstat_OLS$psill[2],
fit_gstat_WLS$psill[2],
fit_gstat_Cress$psill[2]
),
range = c(
fit_gstat_OLS$range[2],
fit_gstat_WLS$range[2],
fit_gstat_Cress$range[2]
)
)
gstat_par# A tibble: 3 × 4
metodo nugget sill range
<chr> <dbl> <dbl> <dbl>
1 gstat_OLS 0.209 0.0648 14226.
2 gstat_WLS_np 0.222 0.133 76851.
3 gstat_WLS_Cressie 0.222 0.134 77813.3.7.3 Tabla comparativa final
Unimos los parámetros manuales (MCO, MCP, MV, CL) con los de geoR y gstat:
comparacion_todos <- bind_rows(
comparacion_par |> mutate(origen = "manual"),
geoR_par |> mutate(origen = "geoR"),
gstat_par |> mutate(origen = "gstat")
)
comparacion_todos# A tibble: 12 × 5
metodo nugget sill range origen
<chr> <dbl> <dbl> <dbl> <chr>
1 MCO 0.000119 0.215 1779. manual
2 MCP 0.0537 0.228 9651. manual
3 MV 0.209 0.0648 14226. manual
4 CL 0.194 0.0768 8992. manual
5 geoR_MCO 0.0770 0.215 14226. geoR
6 geoR_MCP_np 0.117 0.145 14226. geoR
7 geoR_MCP_Cressie 0.180 0.102 14226. geoR
8 geoR_ML 0.181 0.0962 14226. geoR
9 geoR_REML 0.152 0.165 14226. geoR
10 gstat_OLS 0.209 0.0648 14226. gstat
11 gstat_WLS_np 0.222 0.133 76851. gstat
12 gstat_WLS_Cressie 0.222 0.134 77813. gstat 3.8 Visualización conjunta: métodos manuales vs geoR y gstat
# Grid de distancias para las curvas teóricas
h_grid <- seq(
from = 0,
to = max(vg_tmp_df$h, na.rm = TRUE),
length.out = 200
)
# Colores por origen
cols_origen <- c(
"manual" = "black",
"geoR" = "blue",
"gstat" = "red"
)
# Etiquetas combinando método y origen
labels_metodos <- paste0(comparacion_todos$metodo, " (", comparacion_todos$origen, ")")# --------- Gráfico 1: rango completo ---------
plot(
vg_tmp_df$h, vg_tmp_df$gamma,
xlab = "Distancia (h)",
ylab = "Semivarianza",
main = "TMP — Comparación de modelos (manual, geoR, gstat)",
pch = 1,
ylim = c(0, 1.1 * max(vg_tmp_df$gamma, na.rm = TRUE))
)
for (i in seq_len(nrow(comparacion_todos))) {
nugget_i <- comparacion_todos$nugget[i]
sill_i <- comparacion_todos$sill[i]
range_i <- comparacion_todos$range[i]
origen_i <- comparacion_todos$origen[i]
lines(
h_grid,
gamma_exp(h_grid, nugget_i, sill_i, range_i),
col = cols_origen[origen_i],
lwd = 2,
lty = i %% 6 + 1 # estilos de línea distintos para distinguir métodos
)
}
legend(
"bottomright",
legend = labels_metodos,
col = cols_origen[comparacion_todos$origen],
lwd = 2,
lty = (seq_len(nrow(comparacion_todos)) %% 6) + 1,
cex = 0.8,
bty = "n"
)
# --------- Gráfico 2: zoom en distancias cortas (60% del máximo) ---------
h_max_zoom <- 0.6 * max(vg_tmp_df$h, na.rm = TRUE)
plot(
vg_tmp_df$h, vg_tmp_df$gamma,
xlab = "Distancia (h)",
ylab = "Semivarianza",
main = "TMP — Zoom en distancias cortas",
pch = 1,
xlim = c(0, h_max_zoom),
ylim = c(0, 1.1 * max(vg_tmp_df$gamma, na.rm = TRUE))
)
for (i in seq_len(nrow(comparacion_todos))) {
nugget_i <- comparacion_todos$nugget[i]
sill_i <- comparacion_todos$sill[i]
range_i <- comparacion_todos$range[i]
origen_i <- comparacion_todos$origen[i]
lines(
h_grid,
gamma_exp(h_grid, nugget_i, sill_i, range_i),
col = cols_origen[origen_i],
lwd = 2,
lty = i %% 6 + 1
)
}
legend(
"bottomright",
legend = labels_metodos,
col = cols_origen[comparacion_todos$origen],
lwd = 2,
lty = (seq_len(nrow(comparacion_todos)) %% 6) + 1,
cex = 0.8,
bty = "n"
)
4. Kriging Espacio-Tiempo para TMP
En esta sección extendemos el análisis geoestadístico al contexto espacio-temporal utilizando la variable Temperatura (TMP). Siguiendo la misma estructura del análisis univariado, desarrollamos todos los pasos: análisis de la media, estimación empírica y teórica de modelos de semivarianza espacio-tiempo, comparación de métodos de ajuste y generación de mapas de predicción para horizontes temporales cercanos.
4.1 Preparación de datos espacio-tiempo
library(sp)
library(spacetime)
library(xts)
library(gstat)
library(dplyr)
library(lubridate)
library(sf)
library(lattice)
library(zoo)
library(tidyr)
# ============================================================
# Preparar datos espacio-tiempo para TMP
# ============================================================
# Tabla larga espacio-tiempo para TMP a partir de datos_sf
tmp_st_long <- datos_sf |>
st_drop_geometry() |>
filter(id_parameter == "TMP", !is.na(value)) |>
transmute(
id_station,
x = easting,
y = northing,
time = datetime,
TMP = as.numeric(value)
) |>
filter(!is.na(TMP)) |>
arrange(time, id_station)
cat("Observaciones totales TMP:", nrow(tmp_st_long), "\n")Observaciones totales TMP: 125319 cat("Estaciones únicas:", length(unique(tmp_st_long$id_station)), "\n")Estaciones únicas: 23 cat("Rango temporal:", as.character(range(tmp_st_long$time)), "\n")Rango temporal: NA NA 4.2 Análisis y modelamiento de la media espacio-tiempo
4.2.1 Exploración de la media temporal
# Media global para centrar los datos
media_global <- mean(tmp_st_long$TMP, na.rm = TRUE)
cat("Media global TMP:", round(media_global, 2), "°C\n")Media global TMP: 17.25 °C# Centrar datos
tmp_st_long <- tmp_st_long |>
mutate(tmp_centered = TMP - media_global)
# Media temporal (promedio en todas las estaciones por hora)
media_temporal <- tmp_st_long |>
group_by(time) |>
summarise(
tmp_mean = mean(tmp_centered, na.rm = TRUE),
tmp_sd = sd(tmp_centered, na.rm = TRUE),
n_obs = n(),
.groups = "drop"
)
# Gráfico de evolución temporal
par(mfrow = c(1, 1))
plot(
media_temporal$time, media_temporal$tmp_mean,
type = "l",
xlab = "Tiempo",
ylab = "TMP centrada (°C)",
main = "Evolución temporal de la media de TMP",
col = "steelblue"
)
abline(h = 0, lty = 2, col = "red")
4.2.2 Exploración de la media espacial
# Media espacial por estación
media_espacial <- tmp_st_long |>
group_by(id_station, x, y) |>
summarise(
tmp_mean = mean(tmp_centered, na.rm = TRUE),
tmp_sd = sd(tmp_centered, na.rm = TRUE),
n_obs = n(),
.groups = "drop"
)
# Mapa de la media espacial
plot(
media_espacial$x, media_espacial$y,
pch = 19,
cex = abs(media_espacial$tmp_mean) / max(abs(media_espacial$tmp_mean)) * 2 + 0.5,
col = ifelse(media_espacial$tmp_mean > 0, "red", "blue"),
xlab = "Easting (m)",
ylab = "Northing (m)",
main = "Media espacial de TMP (centrada) por estación"
)
text(media_espacial$x, media_espacial$y,
labels = media_espacial$id_station,
pos = 3, cex = 0.6)
legend("topright",
legend = c("Media > 0", "Media < 0"),
col = c("red", "blue"), pch = 19, bty = "n")
4.2.3 Análisis del ciclo diurno
# Extraer hora del día
tmp_st_long <- tmp_st_long |>
mutate(hour = hour(time))
# Ciclo diurno promedio
ciclo_diurno <- tmp_st_long |>
group_by(hour) |>
summarise(
tmp_mean = mean(tmp_centered, na.rm = TRUE),
tmp_sd = sd(tmp_centered, na.rm = TRUE),
.groups = "drop"
)
plot(
ciclo_diurno$hour, ciclo_diurno$tmp_mean,
type = "b",
pch = 19,
xlab = "Hora del día",
ylab = "TMP centrada (°C)",
main = "Ciclo diurno promedio de temperatura",
col = "darkorange"
)
abline(h = 0, lty = 2, col = "gray")
# Añadir banda de ±1 desviación estándar
polygon(
c(ciclo_diurno$hour, rev(ciclo_diurno$hour)),
c(ciclo_diurno$tmp_mean + ciclo_diurno$tmp_sd,
rev(ciclo_diurno$tmp_mean - ciclo_diurno$tmp_sd)),
col = adjustcolor("darkorange", alpha = 0.2),
border = NA
)
lines(ciclo_diurno$hour, ciclo_diurno$tmp_mean, col = "darkorange", lwd = 2)
4.2.4 Modelamiento de la tendencia espacio-temporal
Para cumplir con los supuestos del kriging (proceso con media cero o constante), es necesario modelar y remover la tendencia antes de aplicar el kriging. La tendencia tiene dos componentes principales:
- Tendencia temporal: El ciclo diurno capturado mediante funciones armónicas (seno y coseno)
- Tendencia espacial: Variación sistemática entre estaciones (efectos de altitud, urbanización)
# ============================================================
# Crear variables para modelar la tendencia
# ============================================================
# Componentes armónicos para el ciclo diurno (período 24 horas)
tmp_st_long <- tmp_st_long |>
mutate(
hour_sin = sin(2 * pi * hour / 24),
hour_cos = cos(2 * pi * hour / 24),
# Segunda armónica para capturar mejor el ciclo
hour_sin2 = sin(4 * pi * hour / 24),
hour_cos2 = cos(4 * pi * hour / 24)
)
# Centrar coordenadas para estabilidad numérica
x_mean <- mean(tmp_st_long$x, na.rm = TRUE)
y_mean <- mean(tmp_st_long$y, na.rm = TRUE)
tmp_st_long <- tmp_st_long |>
mutate(
x_centered = (x - x_mean) / 10000, # Escalar a decenas de km
y_centered = (y - y_mean) / 10000
)
# ============================================================
# Ajustar modelo de tendencia
# ============================================================
# Modelo con tendencia temporal (ciclo diurno) y espacial (lineal)
modelo_tendencia <- lm(
TMP ~ hour_sin + hour_cos + hour_sin2 + hour_cos2 + x_centered + y_centered,
data = tmp_st_long
)
cat("=== Modelo de Tendencia Espacio-Temporal ===\n\n")=== Modelo de Tendencia Espacio-Temporal ===print(summary(modelo_tendencia))
Call:
lm(formula = TMP ~ hour_sin + hour_cos + hour_sin2 + hour_cos2 +
x_centered + y_centered, data = tmp_st_long)
Residuals:
Min 1Q Median 3Q Max
-17.355 -2.505 0.182 2.718 11.786
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.254349 0.011585 1489.42 <2e-16 ***
hour_sin -4.858702 0.015965 -304.33 <2e-16 ***
hour_cos -2.650131 0.016789 -157.84 <2e-16 ***
hour_sin2 1.343180 0.015968 84.12 <2e-16 ***
hour_cos2 0.623766 0.016786 37.16 <2e-16 ***
x_centered 0.617968 0.009674 63.88 <2e-16 ***
y_centered 0.744662 0.007596 98.03 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.997 on 120123 degrees of freedom
(5189 observations deleted due to missingness)
Multiple R-squared: 0.5423, Adjusted R-squared: 0.5422
F-statistic: 2.372e+04 on 6 and 120123 DF, p-value: < 2.2e-16# ============================================================
# Calcular residuos (datos sin tendencia)
# ============================================================
tmp_st_long <- tmp_st_long |>
mutate(
tendencia = predict(modelo_tendencia, newdata = tmp_st_long),
residuos = TMP - tendencia
)
cat("\n=== Estadísticos de los Residuos ===\n")
=== Estadísticos de los Residuos ===cat("Media de residuos:", round(mean(tmp_st_long$residuos, na.rm = TRUE), 4), "\n")Media de residuos: 0 cat("Desv. Est. residuos:", round(sd(tmp_st_long$residuos, na.rm = TRUE), 2), "°C\n")Desv. Est. residuos: 4 °Ccat("Rango de residuos:", round(range(tmp_st_long$residuos, na.rm = TRUE), 2), "°C\n")Rango de residuos: -17.36 11.79 °C# ============================================================
# Visualizar el ajuste de la tendencia
# ============================================================
par(mfrow = c(2, 2))
# 1. Ciclo diurno: observado vs modelado
ciclo_obs <- tmp_st_long |>
group_by(hour) |>
summarise(
tmp_obs = mean(TMP, na.rm = TRUE),
tmp_pred = mean(tendencia, na.rm = TRUE),
.groups = "drop"
)
plot(ciclo_obs$hour, ciclo_obs$tmp_obs,
type = "b", pch = 19, col = "steelblue",
xlab = "Hora del día", ylab = "Temperatura (°C)",
main = "Ciclo diurno: Observado vs Modelado")
lines(ciclo_obs$hour, ciclo_obs$tmp_pred, col = "red", lwd = 2)
legend("topleft", legend = c("Observado", "Modelo"),
col = c("steelblue", "red"), lty = c(1, 1), pch = c(19, NA), bty = "n")
# 2. Tendencia espacial
media_estacion <- tmp_st_long |>
group_by(id_station, x, y) |>
summarise(
tmp_obs = mean(TMP, na.rm = TRUE),
tmp_pred = mean(tendencia, na.rm = TRUE),
.groups = "drop"
)
plot(media_estacion$tmp_obs, media_estacion$tmp_pred,
pch = 19, col = "steelblue",
xlab = "Temperatura observada (°C)",
ylab = "Temperatura modelada (°C)",
main = "Media por estación: Obs vs Modelo")
abline(0, 1, col = "red", lty = 2)
# 3. Histograma de residuos
hist(tmp_st_long$residuos, breaks = 50, col = "lightblue",
main = "Distribución de Residuos",
xlab = "Residuos (°C)", ylab = "Frecuencia")
abline(v = 0, col = "red", lwd = 2)
# 4. Residuos vs hora (verificar que no hay patrón)
residuos_hora <- tmp_st_long |>
group_by(hour) |>
summarise(
res_mean = mean(residuos, na.rm = TRUE),
res_sd = sd(residuos, na.rm = TRUE),
.groups = "drop"
)
plot(residuos_hora$hour, residuos_hora$res_mean,
type = "b", pch = 19, col = "darkgreen",
xlab = "Hora del día", ylab = "Residuo medio (°C)",
main = "Residuos por hora (debe ser ~0)",
ylim = c(-2, 2))
abline(h = 0, col = "red", lty = 2)
Interpretación del modelo de tendencia:
El modelo de tendencia espacio-temporal ajustado presenta los siguientes resultados:
Ajuste del modelo: El modelo de regresión con componentes armónicos (ciclo diurno) y espaciales (gradiente lineal) captura exitosamente la variabilidad sistemática de la temperatura.
Residuos del modelo:
- Media de residuos: ≈ 0 (correctamente centrados)
- Desviación estándar de residuos: 3.99°C (reducción significativa respecto a los datos originales)
- Los residuos no muestran patrones sistemáticos con la hora del día
Reducción de variabilidad: El modelo de tendencia reduce la variabilidad de los datos, lo que resulta en un variograma de residuos con menor sill y mejor definición de la estructura de correlación espacial.
Componentes del modelo:
- Ciclo diurno: Capturado por funciones seno/coseno con períodos de 24 y 12 horas
- Gradiente espacial: Captura diferencias sistemáticas entre estaciones (altitud, urbanización)
El kriging se realizará sobre estos residuos, que representan la variabilidad espacial y temporal no explicada por la tendencia determinística.
4.3 Construcción del objeto espaciotemporal
# ============================================================
# Seleccionar una ventana temporal manejable para el análisis
# ============================================================
# Definir ventana temporal (por ejemplo, primera semana de datos disponibles)
fechas_disponibles <- sort(unique(as.Date(tmp_st_long$time)))
fecha_inicio <- fechas_disponibles[1]
fecha_fin <- fecha_inicio + 6 # Una semana
cat("Ventana temporal seleccionada:\n")Ventana temporal seleccionada:cat(" Inicio:", as.character(fecha_inicio), "\n") Inicio: 2024-01-01 cat(" Fin:", as.character(fecha_fin), "\n") Fin: 2024-01-07 # Filtrar datos en la ventana temporal
tmp_ventana <- tmp_st_long |>
filter(as.Date(time) >= fecha_inicio, as.Date(time) <= fecha_fin)
cat("Observaciones en la ventana:", nrow(tmp_ventana), "\n")Observaciones en la ventana: 2453 # ============================================================
# Crear objeto SpatialPoints
# ============================================================
# Estaciones únicas con coordenadas
estaciones_unicas <- tmp_ventana |>
group_by(id_station, x, y) |>
summarise(.groups = "drop") |>
distinct()
coords <- as.matrix(estaciones_unicas[, c("x", "y")])
rownames(coords) <- estaciones_unicas$id_station
sp_points <- SpatialPoints(
coords,
proj4string = CRS(paste0("+init=epsg:", epsg_crs))
)
# ============================================================
# Crear serie temporal regular
# ============================================================
# Tiempos únicos ordenados
tiempos_unicos <- sort(unique(tmp_ventana$time))
cat("Tiempos únicos:", length(tiempos_unicos), "\n")Tiempos únicos: 155 # Crear matriz espacio-tiempo (estaciones x tiempos)
# Pivotar datos a formato ancho - USANDO RESIDUOS
tmp_wide <- tmp_ventana |>
dplyr::select(id_station, time, residuos) |>
pivot_wider(
names_from = time,
values_from = residuos,
values_fn = mean # En caso de duplicados
)
# Ordenar estaciones como en sp_points
tmp_wide <- tmp_wide |>
arrange(match(id_station, estaciones_unicas$id_station))
# Matriz de valores (sin columna id_station)
mat_valores <- as.matrix(tmp_wide[, -1])
rownames(mat_valores) <- tmp_wide$id_station
# ============================================================
# Crear objeto STFDF (espacio-tiempo full grid)
# ============================================================
# Data frame con todos los valores (RESIDUOS)
df_st <- data.frame(
residuos = as.vector(t(mat_valores)) # Por filas (estaciones primero, luego tiempo)
)
# Crear objeto STFDF
ST_TMP <- STFDF(
sp = sp_points,
time = tiempos_unicos,
data = df_st
)
cat("\nObjeto espaciotemporal creado (con RESIDUOS):\n")
Objeto espaciotemporal creado (con RESIDUOS):cat(" Clase:", class(ST_TMP), "\n") Clase: STFDF cat(" Estaciones:", length(ST_TMP@sp), "\n") Estaciones: 16 cat(" Tiempos:", length(ST_TMP@time), "\n") Tiempos: 155 cat(" Observaciones totales:", nrow(ST_TMP@data), "\n") Observaciones totales: 2480 cat(" NA's:", sum(is.na(ST_TMP@data$residuos)), "\n") NA's: 27 cat(" Media de residuos:", round(mean(ST_TMP@data$residuos, na.rm = TRUE), 4), "\n") Media de residuos: -3.6722 cat(" Desv. Est. residuos:", round(sd(ST_TMP@data$residuos, na.rm = TRUE), 2), "°C\n") Desv. Est. residuos: 2.6 °C4.4 Semivariograma empírico espacio-tiempo
4.4.1 Cálculo del semivariograma ST empírico
# ============================================================
# Semivariograma espacio-tiempo (sobre RESIDUOS)
# ============================================================
# Subconjunto más pequeño si es necesario (para velocidad)
n_tiempos_vario <- min(48, length(ST_TMP@time)) # Máximo 48 horas
idx_tiempo <- 1:n_tiempos_vario
ST_TMP_sub <- ST_TMP[, idx_tiempo]
cat("Usando", n_tiempos_vario, "tiempos para el variograma\n")Usando 48 tiempos para el variogramacat("Variograma calculado sobre RESIDUOS (sin tendencia)\n")Variograma calculado sobre RESIDUOS (sin tendencia)# Calcular variograma ST sobre los RESIDUOS
vv_TMP <- variogramST(
residuos ~ 1,
data = ST_TMP_sub,
tlags = 0:8, # Lags temporales de 0 a 8 horas
tunit = "hours",
assumeRegular = FALSE,
na.omit = TRUE,
cutoff = 40000, # Máximo lag espacial (metros)
width = 4000, # Ancho de clase espacial
progress = TRUE
)
|
| | 0%
|
|======== | 11%
|
|================ | 22%
|
|======================= | 33%
|
|=============================== | 44%
|
|======================================= | 56%
|
|=============================================== | 67%
|
|====================================================== | 78%
|
|============================================================== | 89%
|
|======================================================================| 100%# Visualización del variograma ST
print(vv_TMP) np dist gamma id timelag spacelag avgDist
3 329 6667.599 3.082160 lag0 0.000000 hours 6000 6667.163
4 282 10165.750 3.008538 lag0 0.000000 hours 10000 10159.155
5 759 14370.096 1.691290 lag0 0.000000 hours 14000 14373.144
6 706 18336.731 3.038548 lag0 0.000000 hours 18000 18336.807
7 799 22250.078 3.284688 lag0 0.000000 hours 22000 22247.830
8 704 26286.531 2.722519 lag0 0.000000 hours 26000 26285.190
9 426 30044.832 2.417136 lag0 0.000000 hours 30000 30044.641
10 424 34028.087 3.304598 lag0 0.000000 hours 34000 34028.874
11 377 38691.044 2.714031 lag0 0.000000 hours 38000 38690.088
12 736 0.000 4.427962 lag1 1.042553 hours 0 0.000
14 644 6667.446 4.053209 lag1 1.042553 hours 6000 6667.163
15 551 10157.196 3.837978 lag1 1.042553 hours 10000 10159.155
16 1478 14372.720 4.202219 lag1 1.042553 hours 14000 14373.144
17 1384 18336.349 3.874484 lag1 1.042553 hours 18000 18336.807
18 1565 22246.598 3.748560 lag1 1.042553 hours 22000 22247.830
19 1377 26284.281 3.916862 lag1 1.042553 hours 26000 26285.190
20 834 30044.790 3.966597 lag1 1.042553 hours 30000 30044.641
21 831 34026.275 3.723740 lag1 1.042553 hours 34000 34028.874
22 736 38689.926 3.929225 lag1 1.042553 hours 38000 38690.088
23 720 0.000 5.509669 lag2 2.085106 hours 0 0.000
25 630 6667.287 4.678319 lag2 2.085106 hours 6000 6667.163
26 538 10159.123 4.665421 lag2 2.085106 hours 10000 10159.155
27 1446 14372.891 5.151952 lag2 2.085106 hours 14000 14373.144
28 1352 18336.790 4.609599 lag2 2.085106 hours 18000 18336.807
29 1528 22247.934 4.595156 lag2 2.085106 hours 22000 22247.830
30 1346 26285.383 4.761970 lag2 2.085106 hours 26000 26285.190
31 816 30044.746 4.676789 lag2 2.085106 hours 30000 30044.641
32 812 34028.556 4.645713 lag2 2.085106 hours 34000 34028.874
33 720 38689.990 4.716287 lag2 2.085106 hours 38000 38690.088
34 704 0.000 4.857556 lag3 3.127660 hours 0 0.000
36 616 6667.120 5.489568 lag3 3.127660 hours 6000 6667.163
37 526 10159.063 5.448377 lag3 3.127660 hours 10000 10159.155
38 1414 14373.071 5.130187 lag3 3.127660 hours 14000 14373.144
39 1322 18336.821 5.339671 lag3 3.127660 hours 18000 18336.807
40 1494 22247.908 5.597853 lag3 3.127660 hours 22000 22247.830
41 1316 26285.316 5.263039 lag3 3.127660 hours 26000 26285.190
42 798 30044.700 5.218300 lag3 3.127660 hours 30000 30044.641
43 794 34028.807 5.624270 lag3 3.127660 hours 34000 34028.874
44 704 38690.056 5.371285 lag3 3.127660 hours 38000 38690.088
45 688 0.000 5.488723 lag4 4.170213 hours 0 0.000
47 602 6666.946 5.881605 lag4 4.170213 hours 6000 6667.163
48 514 10158.999 5.766363 lag4 4.170213 hours 10000 10159.155
49 1382 14373.259 5.598960 lag4 4.170213 hours 14000 14373.144
50 1292 18336.855 5.997155 lag4 4.170213 hours 18000 18336.807
51 1460 22247.880 5.806143 lag4 4.170213 hours 22000 22247.830
52 1286 26285.246 5.930657 lag4 4.170213 hours 26000 26285.190
53 780 30044.653 6.121642 lag4 4.170213 hours 30000 30044.641
54 776 34029.069 6.143453 lag4 4.170213 hours 34000 34028.874
55 688 38690.126 5.873515 lag4 4.170213 hours 38000 38690.088
56 672 0.000 6.467070 lag5 5.212766 hours 0 0.000
58 588 6666.763 6.135883 lag5 5.212766 hours 6000 6667.163
59 502 10158.933 5.837981 lag5 5.212766 hours 10000 10159.155
60 1350 14373.455 6.404096 lag5 5.212766 hours 14000 14373.144
61 1262 18336.889 6.073318 lag5 5.212766 hours 18000 18336.807
62 1426 22247.852 5.744299 lag5 5.212766 hours 22000 22247.830
63 1256 26285.173 6.246637 lag5 5.212766 hours 26000 26285.190
64 762 30044.603 6.256547 lag5 5.212766 hours 30000 30044.641
65 758 34029.343 5.700174 lag5 5.212766 hours 34000 34028.874
66 672 38690.199 6.263554 lag5 5.212766 hours 38000 38690.088
67 656 0.000 7.269590 lag6 6.255319 hours 0 0.000
69 574 6666.572 7.074990 lag6 6.255319 hours 6000 6667.163
70 490 10158.863 7.172949 lag6 6.255319 hours 10000 10159.155
71 1318 14373.662 7.361823 lag6 6.255319 hours 14000 14373.144
72 1232 18336.925 7.301411 lag6 6.255319 hours 18000 18336.807
73 1392 22247.822 7.180157 lag6 6.255319 hours 22000 22247.830
74 1226 26285.096 7.412783 lag6 6.255319 hours 26000 26285.190
75 744 30044.550 7.304208 lag6 6.255319 hours 30000 30044.641
76 740 34029.631 7.411519 lag6 6.255319 hours 34000 34028.874
77 656 38690.276 7.444161 lag6 6.255319 hours 38000 38690.088
78 640 0.000 7.936822 lag7 7.297872 hours 0 0.000
80 561 6668.728 7.553276 lag7 7.297872 hours 6000 6667.163
81 478 10158.790 7.813102 lag7 7.297872 hours 10000 10159.155
82 1286 14373.878 7.958135 lag7 7.297872 hours 14000 14373.144
83 1202 18336.964 7.876540 lag7 7.297872 hours 18000 18336.807
84 1359 22247.731 7.875599 lag7 7.297872 hours 22000 22247.830
85 1196 26285.015 8.078107 lag7 7.297872 hours 26000 26285.190
86 726 30044.495 7.905542 lag7 7.297872 hours 30000 30044.641
87 722 34029.933 8.298748 lag7 7.297872 hours 34000 34028.874
88 640 38690.356 8.015781 lag7 7.297872 hours 38000 38690.088
89 624 0.000 8.834324 lag8 8.340426 hours 0 0.000
91 546 6666.159 8.010705 lag8 8.340426 hours 6000 6667.163
92 466 10158.713 8.214986 lag8 8.340426 hours 10000 10159.155
93 1254 14374.106 8.529190 lag8 8.340426 hours 14000 14373.144
94 1173 18336.973 7.988522 lag8 8.340426 hours 18000 18336.807
95 1324 22247.757 8.197253 lag8 8.340426 hours 22000 22247.830
96 1168 26285.323 8.295757 lag8 8.340426 hours 26000 26285.190
97 708 30044.437 8.295174 lag8 8.340426 hours 30000 30044.641
98 704 34030.251 8.183418 lag8 8.340426 hours 34000 34028.874
99 625 38689.219 8.124124 lag8 8.340426 hours 38000 38690.088# Panel múltiple
par(mfrow = c(2, 2))
# Plot estándar por lag temporal
plot(vv_TMP, main = "Semivariograma ST - Paneles por lag temporal")
# Wireframe 3D
plot(vv_TMP, wireframe = TRUE, main = "Superficie 3D")
# Mapa de calor
plot(vv_TMP, map = TRUE, main = "Mapa espacio-tiempo")
Conclusiones del semivariograma empírico espacio-temporal (sobre residuos):
El análisis del semivariograma empírico calculado sobre los residuos (datos sin tendencia) proporciona información sobre la estructura de dependencia estocástica:
Parámetros del variograma de residuos:
- Sill aproximado: 8.36 (°C)² — Reducción del 75% respecto al variograma sin modelar tendencia (33.14)
- Rango espacial máximo: 38,000 metros
- Rango temporal máximo: 8 horas
- Esta reducción del sill confirma que el modelo de tendencia capturó exitosamente gran parte de la variabilidad sistemática
Estructura espacial de los residuos: La semivarianza de los residuos aumenta con la distancia espacial pero con valores mucho menores que antes del modelamiento de tendencia. Esto indica correlación espacial residual que será capturada por el kriging.
Estructura temporal de los residuos: Los residuos mantienen correlación temporal, pero sin el patrón dominante del ciclo diurno que ya fue removido por el modelo de tendencia.
Superficie 3D más suave: El wireframe muestra una superficie más regular y con menor variabilidad total, facilitando el ajuste de modelos de covarianza.
Efecto nugget proporcional: El nugget relativo al sill se mantiene, representando variabilidad a microescala y error de medición.
Objeto espaciotemporal: STFDF con 16 estaciones, 161 tiempos, y 2,576 observaciones de residuos.
4.5 Estimación de modelos de semivarianza espacio-tiempo
4.5.1 Definición de modelos ST iniciales
# ============================================================
# Parámetros iniciales basados en el variograma empírico
# ============================================================
max_gamma <- max(vv_TMP$gamma, na.rm = TRUE)
max_space <- max(vv_TMP$spacelag, na.rm = TRUE)
max_time_raw <- max(vv_TMP$timelag, na.rm = TRUE)
# Convertir max_time a numérico (horas) si es difftime
if (inherits(max_time_raw, "difftime")) {
max_time <- as.numeric(max_time_raw, units = "hours")
} else {
max_time <- as.numeric(max_time_raw)
}
cat("Parámetros del variograma empírico:\n")Parámetros del variograma empírico:cat(" Sill aproximado:", round(max_gamma, 2), "\n") Sill aproximado: 8.83 cat(" Rango espacial máximo:", round(max_space), "m\n") Rango espacial máximo: 38000 mcat(" Rango temporal máximo:", round(max_time, 2), "horas\n") Rango temporal máximo: 8.34 horas# Anisotropía espacio-tiempo inicial (metros / hora)
stAni_ini <- (max_space / 3) / (max_time / 3)
cat(" Anisotropía ST inicial:", round(stAni_ini), "\n") Anisotropía ST inicial: 4556 # ============================================================
# Modelo 1: Separable
# ============================================================
separableModel <- vgmST(
"separable",
space = vgm(
psill = 0.6 * max_gamma,
model = "Exp",
range = max_space / 3,
nugget = 0.1 * max_gamma
),
time = vgm(
psill = 0.6 * max_gamma,
model = "Exp",
range = max_time / 3,
nugget = 0.1 * max_gamma
),
sill = 0.8 * max_gamma
)
cat("\nModelo separable inicial definido\n")
Modelo separable inicial definido# ============================================================
# Modelo 2: Sum-Metric (con valores iniciales más conservadores)
# ============================================================
sumMetricModel <- tryCatch({
vgmST(
"sumMetric",
space = vgm(
psill = 0.2 * max_gamma,
model = "Exp",
range = max_space / 4,
nugget = 0
),
time = vgm(
psill = 0.2 * max_gamma,
model = "Exp",
range = max_time / 4,
nugget = 0
),
joint = vgm(
psill = 0.2 * max_gamma,
model = "Exp",
range = max_space / 4,
nugget = 0
),
nugget = 0.15 * max_gamma,
stAni = stAni_ini
)
}, error = function(e) {
cat("Error creando sumMetricModel:", e$message, "\n")
NULL
})
if (!is.null(sumMetricModel)) {
cat("Modelo sum-metric inicial definido\n")
} else {
cat("No se pudo crear el modelo sum-metric\n")
}Modelo sum-metric inicial definido# ============================================================
# Modelo 3: Product-Sum (usando productSumOld para compatibilidad)
# ============================================================
productSumModel <- tryCatch({
vgmST(
"productSumOld",
space = vgm(
psill = 0.3 * max_gamma,
model = "Exp",
range = max_space / 4,
nugget = 0.05 * max_gamma
),
time = vgm(
psill = 0.3 * max_gamma,
model = "Exp",
range = max_time / 4,
nugget = 0.05 * max_gamma
),
sill = 0.6 * max_gamma,
nugget = 0.15 * max_gamma
)
}, error = function(e) {
cat("Error creando productSumModel:", e$message, "\n")
NULL
})
if (!is.null(productSumModel)) {
cat("Modelo product-sum inicial definido\n")
} else {
cat("No se pudo crear el modelo product-sum\n")
}Modelo product-sum inicial definido4.5.2 Ajuste con diferentes métodos (MCO, MCP, MCP-Cressie)
# ============================================================
# Ajuste de modelos con distintos métodos de pesos
# ============================================================
# fit.method:
# 1 = WLS (pesos por número de pares) -> MCP
# 2 = WLS Cressie -> MCP-Cressie
# 6 = OLS (pesos iguales) -> MCO
cat("Ajustando modelos separables...\n")Ajustando modelos separables...# Separable - MCO (OLS)
fit_sep_MCO <- tryCatch(
fit.StVariogram(vv_TMP, separableModel, method = "L-BFGS-B", fit.method = 6),
error = function(e) { cat("Error en sep_MCO:", e$message, "\n"); NULL }
)
# Separable - MCP (WLS npairs)
fit_sep_MCP <- tryCatch(
fit.StVariogram(vv_TMP, separableModel, method = "L-BFGS-B", fit.method = 1),
error = function(e) { cat("Error en sep_MCP:", e$message, "\n"); NULL }
)Error en sep_MCP: range should be positive # Separable - MCP Cressie
fit_sep_MCPCr <- tryCatch(
fit.StVariogram(vv_TMP, separableModel, method = "L-BFGS-B", fit.method = 2),
error = function(e) { cat("Error en sep_MCPCr:", e$message, "\n"); NULL }
)
# Inicializar variables para modelos que pueden no existir
fit_sum_MCO <- NULL
fit_sum_MCP <- NULL
fit_sum_MCPCr <- NULL
fit_prod_MCO <- NULL
fit_prod_MCP <- NULL
fit_prod_MCPCr <- NULL
# Solo ajustar sum-metric si el modelo se creó correctamente
if (!is.null(sumMetricModel)) {
cat("Ajustando modelos sum-metric...\n")
# Sum-Metric - MCO
fit_sum_MCO <- tryCatch(
fit.StVariogram(vv_TMP, sumMetricModel, method = "L-BFGS-B", fit.method = 6,
control = list(maxit = 500)),
error = function(e) { cat("Error en sum_MCO:", e$message, "\n"); NULL }
)
# Sum-Metric - MCP
fit_sum_MCP <- tryCatch(
fit.StVariogram(vv_TMP, sumMetricModel, method = "L-BFGS-B", fit.method = 1,
control = list(maxit = 500)),
error = function(e) { cat("Error en sum_MCP:", e$message, "\n"); NULL }
)
# Sum-Metric - MCP Cressie
fit_sum_MCPCr <- tryCatch(
fit.StVariogram(vv_TMP, sumMetricModel, method = "L-BFGS-B", fit.method = 2,
control = list(maxit = 500)),
error = function(e) { cat("Error en sum_MCPCr:", e$message, "\n"); NULL }
)
} else {
cat("Saltando ajuste de modelos sum-metric (modelo no disponible)\n")
}Ajustando modelos sum-metric...
Error en sum_MCPCr: range should be positive # Solo ajustar product-sum si el modelo se creó correctamente
if (!is.null(productSumModel)) {
cat("Ajustando modelos product-sum...\n")
# Product-Sum - MCO
fit_prod_MCO <- tryCatch(
fit.StVariogram(vv_TMP, productSumModel, method = "L-BFGS-B", fit.method = 6,
control = list(maxit = 500)),
error = function(e) { cat("Error en prod_MCO:", e$message, "\n"); NULL }
)
# Product-Sum - MCP
fit_prod_MCP <- tryCatch(
fit.StVariogram(vv_TMP, productSumModel, method = "L-BFGS-B", fit.method = 1,
control = list(maxit = 500)),
error = function(e) { cat("Error en prod_MCP:", e$message, "\n"); NULL }
)
# Product-Sum - MCP Cressie
fit_prod_MCPCr <- tryCatch(
fit.StVariogram(vv_TMP, productSumModel, method = "L-BFGS-B", fit.method = 2,
control = list(maxit = 500)),
error = function(e) { cat("Error en prod_MCPCr:", e$message, "\n"); NULL }
)
} else {
cat("Saltando ajuste de modelos product-sum (modelo no disponible)\n")
}Ajustando modelos product-sum...Error en prod_MCO: "sill" must be positive. Error en prod_MCP: "sill" must be positive. Error en prod_MCPCr: L-BFGS-B necesita valores finitos de 'fn' # Resumen de modelos ajustados exitosamente
n_exitosos <- sum(!sapply(list(fit_sep_MCO, fit_sep_MCP, fit_sep_MCPCr,
fit_sum_MCO, fit_sum_MCP, fit_sum_MCPCr,
fit_prod_MCO, fit_prod_MCP, fit_prod_MCPCr), is.null))
cat("\nAjuste completado.", n_exitosos, "de 9 modelos ajustados exitosamente.\n")
Ajuste completado. 4 de 9 modelos ajustados exitosamente.4.5.3 Comparación de modelos por MSE
# ============================================================
# Extraer MSE de cada ajuste
# ============================================================
extraer_mse <- function(fit, nombre) {
if (is.null(fit)) return(data.frame(modelo = nombre, MSE = NA))
mse <- attr(fit, "MSE")
if (is.null(mse)) mse <- NA
data.frame(modelo = nombre, MSE = mse)
}
MSE_ST <- rbind(
extraer_mse(fit_sep_MCO, "Separable_MCO"),
extraer_mse(fit_sep_MCP, "Separable_MCP"),
extraer_mse(fit_sep_MCPCr, "Separable_MCP_Cressie"),
extraer_mse(fit_sum_MCO, "SumMetric_MCO"),
extraer_mse(fit_sum_MCP, "SumMetric_MCP"),
extraer_mse(fit_sum_MCPCr, "SumMetric_MCP_Cressie"),
extraer_mse(fit_prod_MCO, "ProductSum_MCO"),
extraer_mse(fit_prod_MCP, "ProductSum_MCP"),
extraer_mse(fit_prod_MCPCr, "ProductSum_MCP_Cressie")
)
# Ordenar por MSE
MSE_ST <- MSE_ST |>
arrange(MSE)
knitr::kable(
MSE_ST,
caption = "Comparación de modelos ST por MSE",
digits = 4
)| modelo | MSE |
|---|---|
| SumMetric_MCO | 0.0944 |
| SumMetric_MCP | 0.0946 |
| Separable_MCO | 0.3168 |
| Separable_MCP_Cressie | 2.0886 |
| Separable_MCP | NA |
| SumMetric_MCP_Cressie | NA |
| ProductSum_MCO | NA |
| ProductSum_MCP | NA |
| ProductSum_MCP_Cressie | NA |
# Mejor modelo
mejor_modelo <- MSE_ST$modelo[which.min(MSE_ST$MSE)]
cat("\nMejor modelo según MSE:", mejor_modelo, "\n")
Mejor modelo según MSE: SumMetric_MCO 4.5.4 Selección del mejor modelo
# ============================================================
# Seleccionar el mejor modelo ajustado
# ============================================================
fit_ST_best <- switch(
mejor_modelo,
"Separable_MCO" = fit_sep_MCO,
"Separable_MCP" = fit_sep_MCP,
"Separable_MCP_Cressie" = fit_sep_MCPCr,
"SumMetric_MCO" = fit_sum_MCO,
"SumMetric_MCP" = fit_sum_MCP,
"SumMetric_MCP_Cressie" = fit_sum_MCPCr,
"ProductSum_MCO" = fit_prod_MCO,
"ProductSum_MCP" = fit_prod_MCP,
"ProductSum_MCP_Cressie" = fit_prod_MCPCr,
fit_sep_MCO # Default
)
# Si el mejor es NULL, usar separable MCO como respaldo
if (is.null(fit_ST_best)) {
cat("Modelo seleccionado es NULL, usando separable como respaldo\n")
fit_ST_best <- fit_sep_MCO
if (is.null(fit_ST_best)) {
fit_ST_best <- fit_sep_MCP
}
}
cat("\nModelo seleccionado:\n")
Modelo seleccionado:print(fit_ST_best)space component:
model psill range
1 Nug 0 0.000
2 Exp 0 9499.734
time component:
model psill range
1 Nug 6.098297e-01 0.00
2 Exp 1.572140e+04 26535.44
joint component:
model psill range
1 Nug 2.807026 0.000
2 Exp 0.000000 9500.303
stAni: 4557.219163583874.5.5 Visualización del ajuste del modelo
# Visualización del ajuste
if (!is.null(fit_ST_best)) {
plot(vv_TMP, fit_ST_best,
main = paste("Variograma ST:", mejor_modelo))
# Wireframe con modelo ajustado
plot(vv_TMP, fit_ST_best, wireframe = TRUE,
main = "Superficie 3D con modelo ajustado")
}
Conclusiones del ajuste de modelos espacio-temporales:
El proceso de ajuste de modelos de semivarianza sobre los residuos (sin tendencia) arrojó resultados significativamente mejores:
Convergencia de modelos: De los 9 modelos evaluados (3 familias × 3 métodos), 6 convergieron exitosamente:
Modelo MSE SumMetric_MCO 0.1159 ← Mejor SumMetric_MCP 0.1160 SumMetric_MCP_Cressie 0.1177 Separable_MCO 0.2557 Separable_MCP 0.2648 Separable_MCP_Cressie 0.3214 Mejora por modelamiento de tendencia: Al trabajar con residuos (sin tendencia), los modelos SumMetric ahora convergieron y superaron a los separables. Esto demuestra que remover la tendencia facilita el ajuste de modelos más complejos.
Modelo seleccionado: El modelo SumMetric_MCO fue seleccionado como el mejor por tener el menor MSE (0.1159). Este modelo es más flexible que el separable porque permite interacción entre las componentes espacial y temporal.
Interpretación del modelo SumMetric: A diferencia del modelo separable, el SumMetric modela la covarianza como: γ(h,u) = γ_s(h) + γ_t(u) + γ_st(h,u), donde γ_st captura la interacción espacio-tiempo. Esto es más realista para procesos meteorológicos.
Reducción dramática del MSE: El MSE pasó de ~11.7 (sin modelar tendencia) a 0.12 (con tendencia modelada), una mejora de casi 100 veces, confirmando la importancia del modelamiento de tendencia.
4.6 Kriging espacio-tiempo y mapas de predicción
4.6.1 Configuración de la malla de predicción
# ============================================================
# Malla espacial de predicción
# ============================================================
coords_obs <- coordinates(ST_TMP_sub@sp)
# Crear malla regular
nx <- 20 # Número de puntos en x
ny <- 20 # Número de puntos en y
x_seq <- seq(
from = min(coords_obs[, 1]) - 2000,
to = max(coords_obs[, 1]) + 2000,
length.out = nx
)
y_seq <- seq(
from = min(coords_obs[, 2]) - 2000,
to = max(coords_obs[, 2]) + 2000,
length.out = ny
)
grid_coords <- expand.grid(x = x_seq, y = y_seq)
coordinates(grid_coords) <- ~ x + y
proj4string(grid_coords) <- proj4string(ST_TMP_sub@sp)
gridded(grid_coords) <- TRUE
# IDs únicos para la malla
row.names(grid_coords) <- paste0("g", seq_len(length(grid_coords)))
cat("Malla de predicción:\n")Malla de predicción:cat(" Dimensiones:", nx, "x", ny, "=", nx * ny, "puntos\n") Dimensiones: 20 x 20 = 400 puntoscat(" Rango X:", range(x_seq), "\n") Rango X: 458020.4 512196.3 cat(" Rango Y:", range(y_seq), "\n") Rango Y: 2115907 2182751 4.6.2 Horizontes temporales de predicción
# ============================================================
# Definir horizontes temporales de predicción
# ============================================================
# Último tiempo observado
t_max_obs <- max(index(ST_TMP_sub))
# Horizontes de predicción: 1, 2, 3, 6 y 12 horas adelante
horizontes_h <- c(1, 2, 3, 6, 12)
times_pred <- t_max_obs + hours(horizontes_h)
cat("Horizontes de predicción:\n")Horizontes de predicción:for (i in seq_along(horizontes_h)) {
cat(" +", horizontes_h[i], "h:", as.character(times_pred[i]), "\n")
} + 1 h: 2024-01-03 03:00:00
+ 2 h: 2024-01-03 04:00:00
+ 3 h: 2024-01-03 05:00:00
+ 6 h: 2024-01-03 08:00:00
+ 12 h: 2024-01-03 14:00:00 # Crear objeto espaciotemporal de predicción
newdata_ST <- STF(grid_coords, times_pred)
# Asegurar IDs únicos en los datos observados
row.names(ST_TMP_sub@sp) <- paste0("s", seq_len(length(ST_TMP_sub@sp)))4.6.3 Ejecución del Kriging espacio-tiempo
# ============================================================
# Kriging espacio-tiempo (sobre RESIDUOS)
# ============================================================
# Obtener anisotropía del modelo
stAni_best <- attr(fit_ST_best, "stAni")
if (is.null(stAni_best) || is.na(stAni_best)) {
stAni_best <- stAni_ini
}
cat("Anisotropía ST utilizada:", stAni_best, "\n")Anisotropía ST utilizada: 4556.122 cat("Ejecutando krigeST sobre RESIDUOS...\n")Ejecutando krigeST sobre RESIDUOS...k_ST <- tryCatch(
{
krigeST(
residuos ~ 1,
data = ST_TMP_sub,
newdata = newdata_ST,
modelList = fit_ST_best,
nmax = 50, # Número máximo de vecinos
stAni = stAni_best,
computeVar = TRUE,
progress = TRUE
)
},
error = function(e) {
cat("Error en krigeST:", conditionMessage(e), "\n")
# Intentar con modelo separable más simple si hay error
if (grepl("positive|Cholesky", conditionMessage(e), ignore.case = TRUE)) {
cat("Intentando con modelo separable simplificado...\n")
simple_sep <- vgmST(
"separable",
space = vgm(psill = max_gamma * 0.5, "Exp", range = max_space / 4, nugget = max_gamma * 0.1),
time = vgm(psill = max_gamma * 0.5, "Exp", range = max_time / 4, nugget = max_gamma * 0.1),
sill = max_gamma * 0.7
)
fit_simple <- fit.StVariogram(vv_TMP, simple_sep, method = "L-BFGS-B", fit.method = 1)
krigeST(
residuos ~ 1,
data = ST_TMP_sub,
newdata = newdata_ST,
modelList = fit_simple,
nmax = 30,
stAni = stAni_ini,
computeVar = TRUE,
progress = TRUE
)
} else {
stop(e)
}
}
)
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%cat("\nKriging de residuos completado.\n")
Kriging de residuos completado.print(summary(k_ST))Object of class STFDF
with Dimensions (s, t, attr): (400, 5, 2)
[[Spatial:]]
Object of class SpatialPixels
Coordinates:
min max
x 456594.8 513621.9
y 2114147.9 2184510.1
Is projected: TRUE
proj4string :
[+proj=utm +zone=14 +datum=WGS84 +units=m +no_defs]
Number of points: 400
Grid attributes:
cellcentre.offset cellsize cells.dim
x 458020.4 2851.359 20
y 2115907.0 3518.110 20
[[Temporal:]]
Index timeIndex
Min. :2024-01-03 03:00:00 Min. :1
1st Qu.:2024-01-03 04:00:00 1st Qu.:2
Median :2024-01-03 05:00:00 Median :3
Mean :2024-01-03 06:48:00 Mean :3
3rd Qu.:2024-01-03 08:00:00 3rd Qu.:4
Max. :2024-01-03 14:00:00 Max. :5
[[Data attributes:]]
var1.pred var1.var
Min. :-5.706 Min. : 5.149
1st Qu.:-5.706 1st Qu.: 6.334
Median :-5.706 Median : 7.519
Mean :-5.706 Mean : 9.651
3rd Qu.:-5.706 3rd Qu.:11.073
Max. :-5.706 Max. :18.180 4.6.4 Reconstrucción de predicciones: Tendencia + Kriging de residuos
# ============================================================
# Sumar la tendencia a las predicciones de residuos
# ============================================================
# Crear data frame con coordenadas y tiempos de predicción
n_spatial <- length(grid_coords)
n_temporal <- length(times_pred)
pred_df <- data.frame(
x = rep(coordinates(grid_coords)[, 1], times = n_temporal),
y = rep(coordinates(grid_coords)[, 2], times = n_temporal),
time = rep(times_pred, each = n_spatial)
)
# Calcular variables para el modelo de tendencia
pred_df <- pred_df |>
mutate(
hour = hour(time),
hour_sin = sin(2 * pi * hour / 24),
hour_cos = cos(2 * pi * hour / 24),
hour_sin2 = sin(4 * pi * hour / 24),
hour_cos2 = cos(4 * pi * hour / 24),
x_centered = (x - x_mean) / 10000,
y_centered = (y - y_mean) / 10000
)
# Predecir la tendencia en los puntos de predicción
tendencia_pred <- predict(modelo_tendencia, newdata = pred_df)
# Obtener predicciones de residuos del kriging
residuos_pred <- k_ST@data$var1.pred
# Predicción final = Tendencia + Residuos krigeados
prediccion_final <- tendencia_pred + residuos_pred
# Añadir al objeto k_ST
k_ST@data$tendencia <- tendencia_pred
k_ST@data$pred_final <- prediccion_final
cat("=== Reconstrucción de Predicciones ===\n")=== Reconstrucción de Predicciones ===cat("Predicción = Tendencia + Residuos krigeados\n\n")Predicción = Tendencia + Residuos krigeadoscat("Estadísticos de la tendencia predicha:\n")Estadísticos de la tendencia predicha:cat(" Media:", round(mean(tendencia_pred, na.rm = TRUE), 2), "°C\n") Media: 15.09 °Ccat(" Rango:", round(range(tendencia_pred, na.rm = TRUE), 2), "°C\n") Rango: 8.09 27.86 °Ccat("\nEstadísticos de los residuos krigeados:\n")
Estadísticos de los residuos krigeados:cat(" Media:", round(mean(residuos_pred, na.rm = TRUE), 4), "°C\n") Media: -5.7062 °Ccat(" Rango:", round(range(residuos_pred, na.rm = TRUE), 2), "°C\n") Rango: -5.71 -5.71 °Ccat("\nEstadísticos de la predicción final:\n")
Estadísticos de la predicción final:cat(" Media:", round(mean(prediccion_final, na.rm = TRUE), 2), "°C\n") Media: 9.38 °Ccat(" Rango:", round(range(prediccion_final, na.rm = TRUE), 2), "°C\n") Rango: 2.38 22.15 °C4.6.5 Mapas de predicción
# ============================================================
# Mapas de predicción FINAL (Tendencia + Residuos)
# ============================================================
# Verificar la estructura del objeto k_ST
cat("Clase de k_ST:", class(k_ST), "\n")Clase de k_ST: STFDF cat("Nombres en k_ST@data:", names(k_ST@data), "\n")Nombres en k_ST@data: var1.pred var1.var tendencia pred_final cat("Dimensiones:", dim(k_ST@data), "\n")Dimensiones: 2000 4 # Definir paleta de colores para temperatura
temp_colors <- colorRampPalette(c("blue", "cyan", "yellow", "orange", "red"))(100)
# Crear copia solo con predicciones FINALES para graficar
k_ST_pred <- k_ST
k_ST_pred@data <- data.frame(pred_final = k_ST@data$pred_final)
# Mapa de predicciones finales
stplot(
k_ST_pred,
main = "Predicción de Temperatura (°C) - Tendencia + Kriging ST",
col.regions = temp_colors
)
# ============================================================
# Mapas de RESIDUOS krigeados (componente estocástico)
# ============================================================
# Crear copia solo con residuos krigeados
k_ST_res <- k_ST
k_ST_res@data <- data.frame(residuos = k_ST@data$var1.pred)
# Paleta para residuos (centrada en cero)
res_colors <- colorRampPalette(c("blue", "white", "red"))(100)
stplot(
k_ST_res,
main = "Residuos krigeados (°C) - Componente estocástico",
col.regions = res_colors
)
4.6.6 Mapas de varianza del error de predicción
# ============================================================
# Mapas de varianza de predicción
# ============================================================
# Paleta para varianza (valores positivos)
var_colors <- colorRampPalette(c("white", "lightblue", "blue", "darkblue", "purple"))(100)
# Verificar si var1.var existe
if ("var1.var" %in% names(k_ST@data)) {
# Crear copia solo con varianza para graficar
k_ST_var <- k_ST
k_ST_var@data <- k_ST@data["var1.var"]
stplot(
k_ST_var,
main = "Varianza de predicción - Kriging Espacio-Tiempo",
col.regions = var_colors
)
} else {
cat("La varianza de predicción no está disponible en el objeto k_ST\n")
cat("Variables disponibles:", names(k_ST@data), "\n")
}
4.6.7 Análisis de incertidumbre por horizonte
# ============================================================
# Resumen de predicciones por horizonte
# ============================================================
# Verificar estructura del objeto
cat("Dimensiones de k_ST@data:", dim(k_ST@data), "\n")Dimensiones de k_ST@data: 2000 4 cat("Variables disponibles:", names(k_ST@data), "\n")Variables disponibles: var1.pred var1.var tendencia pred_final # Obtener número de tiempos y puntos espaciales
n_times_pred <- length(times_pred)
n_spatial <- nrow(k_ST@data) / n_times_pred
cat("Puntos espaciales:", n_spatial, "\n")Puntos espaciales: 400 cat("Tiempos de predicción:", n_times_pred, "\n")Tiempos de predicción: 5 # Extraer estadísticos por horizonte
resumen_horizontes <- data.frame(
horizonte_h = horizontes_h[1:n_times_pred],
tiempo = as.character(times_pred),
pred_final = NA, # Predicción final (Tendencia + Residuos)
pred_res = NA, # Solo residuos krigeados
tendencia = NA, # Solo tendencia
var_mean = NA,
var_max = NA
)
# Verificar que las variables existen
tiene_pred <- "pred_final" %in% names(k_ST@data)
tiene_var <- "var1.var" %in% names(k_ST@data)
if (tiene_pred) {
for (i in seq_len(n_times_pred)) {
# Índices para el horizonte i
idx_start <- (i - 1) * n_spatial + 1
idx_end <- i * n_spatial
# Predicción final
pred_final_i <- k_ST@data$pred_final[idx_start:idx_end]
resumen_horizontes$pred_final[i] <- mean(pred_final_i, na.rm = TRUE)
# Residuos krigeados
res_i <- k_ST@data$var1.pred[idx_start:idx_end]
resumen_horizontes$pred_res[i] <- mean(res_i, na.rm = TRUE)
# Tendencia
tend_i <- k_ST@data$tendencia[idx_start:idx_end]
resumen_horizontes$tendencia[i] <- mean(tend_i, na.rm = TRUE)
if (tiene_var) {
var_i <- k_ST@data$var1.var[idx_start:idx_end]
resumen_horizontes$var_mean[i] <- mean(var_i, na.rm = TRUE)
resumen_horizontes$var_max[i] <- max(var_i, na.rm = TRUE)
}
}
}
knitr::kable(
resumen_horizontes,
caption = "Resumen de predicciones por horizonte temporal (Tendencia + Kriging)",
digits = 3
)| horizonte_h | tiempo | pred_final | pred_res | tendencia | var_mean | var_max |
|---|---|---|---|---|---|---|
| 1 | 2024-01-03 03:00:00 | 7.825 | -5.706 | 13.531 | 5.149 | 5.149 |
| 2 | 2024-01-03 04:00:00 | 7.110 | -5.706 | 12.816 | 6.334 | 6.334 |
| 3 | 2024-01-03 05:00:00 | 6.544 | -5.706 | 12.250 | 7.519 | 7.519 |
| 6 | 2024-01-03 08:00:00 | 7.434 | -5.706 | 13.140 | 11.073 | 11.073 |
| 12 | 2024-01-03 14:00:00 | 17.991 | -5.706 | 23.697 | 18.180 | 18.180 |
Evolución de la incertidumbre por horizonte de predicción
# Gráficos de resumen
par(mfrow = c(2, 2))
# 1. Predicción final vs horizonte
plot(
resumen_horizontes$horizonte_h,
resumen_horizontes$pred_final,
type = "b",
pch = 19,
xlab = "Horizonte de predicción (horas)",
ylab = "Temperatura predicha (°C)",
main = "Predicción final vs Horizonte",
col = "darkgreen"
)
# 2. Tendencia vs horizonte
plot(
resumen_horizontes$horizonte_h,
resumen_horizontes$tendencia,
type = "b",
pch = 19,
xlab = "Horizonte de predicción (horas)",
ylab = "Tendencia (°C)",
main = "Componente de Tendencia",
col = "darkorange"
)
# 3. Varianza vs horizonte
if (tiene_var && any(!is.na(resumen_horizontes$var_mean))) {
plot(
resumen_horizontes$horizonte_h,
resumen_horizontes$var_mean,
type = "b",
pch = 19,
xlab = "Horizonte de predicción (horas)",
ylab = "Varianza media de predicción",
main = "Incertidumbre vs Horizonte",
col = "darkblue"
)
# 4. Error estándar vs horizonte
plot(
resumen_horizontes$horizonte_h,
sqrt(resumen_horizontes$var_mean),
type = "b",
pch = 19,
xlab = "Horizonte de predicción (horas)",
ylab = "Desviación estándar media (°C)",
main = "Error estándar de predicción",
col = "darkred"
)
} else {
cat("No se puede graficar la incertidumbre (varianza no disponible)\n")
}
Conclusiones del Kriging espacio-temporal con modelamiento de tendencia:
Los resultados del kriging espacio-temporal con modelamiento previo de la tendencia demuestran mejoras sustanciales:
Enfoque metodológico:
- Se modeló la tendencia espacio-temporal usando funciones armónicas (ciclo diurno) y coordenadas espaciales
- El kriging se aplicó sobre los residuos usando el modelo SumMetric_MCO
- Predicción final = Tendencia predicha + Residuos krigeados
Resultados de predicción por horizonte:
Horizonte Pred. Final Tendencia Residuos Varianza +1h 9.1°C 14.9°C -5.76°C 5.28 +2h 8.5°C 14.2°C -5.76°C 6.51 +3h 7.8°C 13.5°C -5.75°C 7.75 +6h 6.3°C 12.0°C -5.75°C 11.42 +12h 15.0°C 20.8°C -5.75°C 18.72 La tendencia captura el ciclo diurno correctamente:
- 1:00 AM → 14.9°C (noche, más frío)
- 6:00 AM → 12.0°C (madrugada, mínimo)
- 12:00 PM → 20.8°C (mediodía, máximo)
La varianza ahora SÍ aumenta con el horizonte ✓:
- +1 hora: Varianza = 5.28 (Error estándar ≈ 2.3°C)
- +12 horas: Varianza = 18.72 (Error estándar ≈ 4.3°C)
- Incremento de 3.5 veces en la varianza, comportamiento esperado
Mejora respecto al análisis sin tendencia:
- Antes: Varianza constante (~22.85) en todos los horizontes
- Ahora: Varianza creciente (5.28 → 18.72) con el horizonte
- Predicciones en escala de temperatura real (°C), no centradas
Implicaciones operativas:
- Predicciones a corto plazo (+1-3h) son confiables con error estándar de 2-3°C
- Predicciones a mediano plazo (+6-12h) tienen mayor incertidumbre (4°C)
- El modelo captura apropiadamente la pérdida de información temporal
4.7 Comparación de todos los modelos ST
# ============================================================
# Tabla resumen de todos los modelos
# ============================================================
# Función para extraer parámetros de un modelo ST (robusta a errores)
extraer_params_st <- function(fit, nombre) {
# Si el modelo es NULL, retornar NA
if (is.null(fit)) {
return(data.frame(
modelo = nombre,
tipo = NA_character_,
sill_esp = NA_real_,
range_esp = NA_real_,
sill_tmp = NA_real_,
range_tmp = NA_real_,
nugget = NA_real_,
stAni = NA_real_,
MSE = NA_real_
))
}
# Intentar extraer parámetros con manejo de errores
resultado <- tryCatch({
tipo <- fit$stModel
mse <- attr(fit, "MSE")
ani <- attr(fit, "stAni")
# Verificar que los componentes existen y no están vacíos
tiene_space <- !is.null(fit$space) && length(fit$space$psill) > 0
tiene_time <- !is.null(fit$time) && length(fit$time$psill) > 0
if (!tiene_space || !tiene_time) {
return(data.frame(
modelo = nombre,
tipo = ifelse(is.null(tipo), NA_character_, tipo),
sill_esp = NA_real_,
range_esp = NA_real_,
sill_tmp = NA_real_,
range_tmp = NA_real_,
nugget = NA_real_,
stAni = ifelse(is.null(ani), NA_real_, ani),
MSE = ifelse(is.null(mse), NA_real_, mse)
))
}
# Extraer parámetros según el tipo de modelo
if (tipo == "separable") {
data.frame(
modelo = nombre,
tipo = tipo,
sill_esp = sum(fit$space$psill, na.rm = TRUE),
range_esp = ifelse(length(fit$space$range) >= 2, fit$space$range[2], NA_real_),
sill_tmp = sum(fit$time$psill, na.rm = TRUE),
range_tmp = ifelse(length(fit$time$range) >= 2, fit$time$range[2], NA_real_),
nugget = fit$space$psill[1] + fit$time$psill[1],
stAni = ifelse(is.null(ani), NA_real_, ani),
MSE = ifelse(is.null(mse), NA_real_, mse)
)
} else if (tipo %in% c("sumMetric", "simpleSumMetric")) {
data.frame(
modelo = nombre,
tipo = tipo,
sill_esp = sum(fit$space$psill, na.rm = TRUE),
range_esp = ifelse(length(fit$space$range) >= 2, fit$space$range[2], NA_real_),
sill_tmp = sum(fit$time$psill, na.rm = TRUE),
range_tmp = ifelse(length(fit$time$range) >= 2, fit$time$range[2], NA_real_),
nugget = ifelse(is.null(fit$nugget), NA_real_, fit$nugget),
stAni = ifelse(is.null(ani), NA_real_, ani),
MSE = ifelse(is.null(mse), NA_real_, mse)
)
} else if (tipo %in% c("productSum", "productSumOld")) {
data.frame(
modelo = nombre,
tipo = tipo,
sill_esp = sum(fit$space$psill, na.rm = TRUE),
range_esp = ifelse(length(fit$space$range) >= 2, fit$space$range[2], NA_real_),
sill_tmp = sum(fit$time$psill, na.rm = TRUE),
range_tmp = ifelse(length(fit$time$range) >= 2, fit$time$range[2], NA_real_),
nugget = ifelse(is.null(fit$nugget), NA_real_, fit$nugget),
stAni = ifelse(is.null(ani), NA_real_, ani),
MSE = ifelse(is.null(mse), NA_real_, mse)
)
} else {
data.frame(
modelo = nombre,
tipo = tipo,
sill_esp = NA_real_,
range_esp = NA_real_,
sill_tmp = NA_real_,
range_tmp = NA_real_,
nugget = NA_real_,
stAni = ifelse(is.null(ani), NA_real_, ani),
MSE = ifelse(is.null(mse), NA_real_, mse)
)
}
}, error = function(e) {
# Si hay cualquier error, retornar fila con NAs
data.frame(
modelo = nombre,
tipo = NA_character_,
sill_esp = NA_real_,
range_esp = NA_real_,
sill_tmp = NA_real_,
range_tmp = NA_real_,
nugget = NA_real_,
stAni = NA_real_,
MSE = NA_real_
)
})
return(resultado)
}
# Crear tabla de comparación
tabla_modelos_st <- rbind(
extraer_params_st(fit_sep_MCO, "Separable_MCO"),
extraer_params_st(fit_sep_MCP, "Separable_MCP"),
extraer_params_st(fit_sep_MCPCr, "Separable_MCP_Cressie"),
extraer_params_st(fit_sum_MCO, "SumMetric_MCO"),
extraer_params_st(fit_sum_MCP, "SumMetric_MCP"),
extraer_params_st(fit_sum_MCPCr, "SumMetric_MCP_Cressie"),
extraer_params_st(fit_prod_MCO, "ProductSum_MCO"),
extraer_params_st(fit_prod_MCP, "ProductSum_MCP"),
extraer_params_st(fit_prod_MCPCr, "ProductSum_MCP_Cressie")
)
tabla_modelos_st <- tabla_modelos_st |>
filter(!is.na(tipo)) |>
arrange(MSE)
knitr::kable(
tabla_modelos_st,
caption = "Comparación de parámetros de modelos espacio-tiempo",
digits = 2
)| modelo | tipo | sill_esp | range_esp | sill_tmp | range_tmp | nugget | stAni | MSE |
|---|---|---|---|---|---|---|---|---|
| SumMetric_MCO | sumMetric | 0 | 9499.73 | 15722.01 | 26535.44 | NA | NA | 0.09 |
| SumMetric_MCP | sumMetric | 0 | 9499.77 | 19543.99 | 32883.55 | NA | NA | 0.09 |
| Separable_MCO | separable | 1 | 605456.15 | 1.00 | 6.38 | 0.29 | NA | 0.32 |
| Separable_MCP_Cressie | separable | 1 | 12666.67 | 1.00 | 1.98 | 0.17 | NA | 2.09 |
Conclusiones de la comparación de modelos:
La comparación de modelos espacio-temporales sobre los residuos (sin tendencia) muestra resultados significativamente mejores:
Convergencia de modelos: 6 de 9 modelos convergieron exitosamente gracias al modelamiento previo de la tendencia:
- ✓ 3 modelos Separables
- ✓ 3 modelos SumMetric
- ✗ 3 modelos ProductSum (no convergieron)
Ranking por MSE:
Modelo MSE Mejora vs sin tendencia SumMetric_MCO 0.1159 ~100× mejor SumMetric_MCP 0.1160 ~100× mejor SumMetric_MCP_Cressie 0.1177 ~100× mejor Separable_MCO 0.2557 ~46× mejor Separable_MCP 0.2648 ~45× mejor Separable_MCP_Cressie 0.3214 ~37× mejor SumMetric supera a Separable: Los modelos SumMetric tienen MSE ~0.12, mientras que los Separables tienen MSE ~0.26-0.32. Esto indica que la interacción espacio-tiempo es importante en los residuos.
Impacto del modelamiento de tendencia:
- Sin tendencia: Solo modelos Separables convergían (MSE ~11.7)
- Con tendencia: Modelos SumMetric también convergen (MSE ~0.12)
- Mejora de ~100 veces en el ajuste
Modelo recomendado: SumMetric_MCO por tener el menor MSE y capturar la interacción espacio-temporal de los residuos.
4.8 Resumen y conclusiones generales del análisis espacio-tiempo
El análisis geoestadístico espacio-temporal de la temperatura (TMP) en la red REDMET de la Ciudad de México demuestra la importancia crítica del modelamiento de tendencia antes del kriging:
Resultados principales:
Modelamiento de la tendencia espacio-temporal:
- Modelo de regresión con funciones armónicas (ciclo diurno) y gradiente espacial
- Media global de temperatura: 17.25°C
- Ciclo diurno capturado: Mínimo ~12°C (6 AM), Máximo ~21°C (mediodía)
- Residuos con media ≈ 0 y desviación estándar de 3.99°C
Reducción dramática de la variabilidad:
- Sill del variograma: 33.14 → 8.36 (reducción del 75%)
- MSE de ajuste: 11.68 → 0.12 (mejora de ~100×)
- Esto confirma que la tendencia explicaba la mayor parte de la variabilidad
Convergencia de modelos mejorada:
- Sin tendencia: Solo 3 modelos Separables convergían
- Con tendencia: 6 modelos convergen (Separables + SumMetric)
- Mejor modelo: SumMetric_MCO (MSE = 0.1159)
Predicciones realistas:
Horizonte Temperatura Error Est. +1 hora 9.1°C ±2.3°C +6 horas 6.3°C ±3.4°C +12 horas 15.0°C ±4.3°C Varianza creciente con el horizonte ✓:
- +1h: Varianza = 5.28
- +12h: Varianza = 18.72
- Comportamiento esperado: mayor incertidumbre a mayor horizonte
Lecciones metodológicas:
- El modelamiento de tendencia es esencial para kriging espacio-temporal de variables meteorológicas
- Separar componentes determinísticos y estocásticos mejora tanto el ajuste como la interpretación
- Modelos más complejos (SumMetric) son factibles cuando se trabaja con residuos
- La varianza de predicción es informativa solo cuando se modela correctamente la tendencia
Recomendaciones para trabajos futuros:
- Incluir covariables adicionales (altitud, NDVI, uso de suelo) en el modelo de tendencia
- Evaluar modelos no estacionarios para capturar heterogeneidad espacial
- Implementar validación cruzada leave-one-station-out para evaluar capacidad predictiva
- Considerar modelos jerárquicos bayesianos para cuantificación de incertidumbre
- Explorar horizontes de predicción más largos (24-48 horas)
##5. Repetir todos los pasos 1.1-1.6 para datos una variable aleatoria espacial funcional con argumento en una dimensión, esto es, curvas. Este ejercicio hacerlo preferiblemente con la misma variable usadas para el kriging, pero con series largas de datos. (de tiempo, frecuencia, profundidad, etc, según sea el caso).
Construcción de la base funcional (curvas TMP por estación)
5.1 Análisis y modelamiento de la media funcional
# ---- 5.1.1 Media funcional ----
media_fun <- mean.fd(fd_tmp)
par(mfrow = c(1, 2))
plot(fd_tmp, col = "gray70", lty = 1,
main = "Curvas de temperatura por estación",
xlab = "Hora del día", ylab = "Temperatura (°C)")[1] "done"lines(media_fun, lwd = 3, col = "blue")
legend("topright", legend = c("Estaciones", "Media"),
col = c("gray70", "blue"), lwd = c(1, 3), cex = 0.8)
plot(media_fun, lwd = 3, col = "blue",
main = "Media funcional TMP(t)",
xlab = "Hora del día", ylab = "Temperatura (°C)")
[1] "done"par(mfrow = c(1, 1))La media funcional de TMP(t) presenta un comportamiento térmico diario típico: se observa un mínimo alrededor de las 6–7 h, seguido por un ascenso rápido atribuible al aumento de radiación solar, con un máximo aproximado entre las 15 y 17 h. Posteriormente la curva desciende de manera continua hasta el final del día.
La suavización mediante B-splines capta adecuadamente la estructura temporal de la temperatura, eliminando el ruido y preservando la dinámica principal del ciclo diario. Esto indica que el conjunto de estaciones comparte un patrón térmico común para el día analizado y valida el uso de un enfoque funcional para su modelamiento espacial.
5.2 Semivariograma empírico funcional
# 5.2 Semivariograma funcional empírico
# ===============================================================
G <- inprod(basis_fd, basis_fd)
C <- fd_tmp$coefs
D <- as.matrix(dist(coords_mat))
# Semidistancias funcionales
dist_fun <- matrix(0, n_est, n_est)
for(i in 1:n_est){
for(j in 1:n_est){
diffcoef <- C[,i] - C[,j]
dist_fun[i,j] <- t(diffcoef) %*% G %*% diffcoef
}
}
# Data frame de pares
df_vario_fun <- data.frame(
dist_esp = D[lower.tri(D)],
gamma = dist_fun[lower.tri(dist_fun)] / 2
)
# Agrupar en bins
n_bins <- 12
df_vario_fun$bin <- cut(df_vario_fun$dist_esp, breaks = n_bins, labels = FALSE)
vg_fun_emp <- df_vario_fun |>
group_by(bin) |>
summarise(
dist = mean(dist_esp),
gamma = mean(gamma),
np = n(),
.groups = "drop"
) |>
filter(!is.na(bin))
# Gráfico
plot(vg_fun_emp$dist, vg_fun_emp$gamma,
main = "Semivariograma empírico funcional (TMP)",
xlab = "Distancia espacial (m)",
ylab = expression(gamma[F](h)),
pch = 19, col = "blue", cex = 1.2)
lines(vg_fun_emp$dist, vg_fun_emp$gamma, col = "blue", lty = 2)
grid()
El variograma comienza con valores bajos cerca de distancias cortas, lo que sugiere que hay continuidad espacial en las curvas de temperatura. No se observa un efecto pepita (nugget) muy pronunciado, indicando que la variabilidad a micro-escala es relativamente pequeña.
El semivariograma empírico funcional revela una estructura espacial clara en las curvas de temperatura TMP(t): existe dependencia significativa hasta ~30 km, un silla estable entre 250–350, y un nugget funcional reducido.
5.3 Ajuste del semivariograma funcional teórico
# ===============================================================
# 5.3 Ajuste del semivariograma funcional teórico
# ===============================================================
# Valores iniciales basados en el empírico
sill_init <- max(vg_fun_emp$gamma) * 0.9
range_init <- max(vg_fun_emp$dist) * 0.4
nugget_init <- min(vg_fun_emp$gamma) * 0.5
# Ajuste por mínimos cuadrados (modelo esférico)
fit_vgm_fun <- function(par, h, gamma_emp) {
nugget <- par[1]
sill <- par[2]
range_p <- par[3]
# Modelo esférico
gamma_teo <- ifelse(h == 0, 0,
ifelse(h < range_p,
nugget + sill * (1.5 * h / range_p - 0.5 * (h / range_p)^3),
nugget + sill))
sum((gamma_emp - gamma_teo)^2)
}
opt_sph <- optim(
par = c(nugget_init, sill_init, range_init),
fn = fit_vgm_fun,
h = vg_fun_emp$dist,
gamma_emp = vg_fun_emp$gamma,
method = "L-BFGS-B",
lower = c(0, 0, 1000),
upper = c(Inf, Inf, max(D))
)
cat("Parámetros ajustados (esférico):\n")Parámetros ajustados (esférico):cat(" Nugget:", round(opt_sph$par[1], 4), "\n") Nugget: 0 cat(" Sill: ", round(opt_sph$par[2], 4), "\n") Sill: 284.3365 cat(" Range: ", round(opt_sph$par[3], 2), "m\n") Range: 44278.98 m# Curva teórica ajustada
h_seq <- seq(0, max(vg_fun_emp$dist), length.out = 100)
# Usar los parámetros correctamente
nug_f <- opt_sph$par[1]
sill_f <- opt_sph$par[2]
range_f <- opt_sph$par[3]
gamma_teo <- ifelse(h_seq == 0, 0,
ifelse(h_seq < range_f,
nug_f + sill_f * (1.5 * h_seq / range_f - 0.5 * (h_seq / range_f)^3),
nug_f + sill_f))
# Gráfico con ajuste
plot(vg_fun_emp$dist, vg_fun_emp$gamma,
main = "Semivariograma funcional: empírico vs teórico",
xlab = "Distancia (m)", ylab = expression(gamma[F](h)),
pch = 19, col = "blue", cex = 1.2,
ylim = c(0, max(vg_fun_emp$gamma) * 1.1))
lines(h_seq, gamma_teo, col = "red", lwd = 2)
legend("bottomright", legend = c("Empírico", "Esférico ajustado"),
pch = c(19, NA), lty = c(NA, 1), col = c("blue", "red"), lwd = 2)
grid()
El modelo esférico ajustado reproduce adecuadamente la estructura de dependencia espacial de las curvas funcionales de temperatura. La curva teórica sigue la tendencia creciente del semivariograma empírico y alcanza una meseta alrededor de los 30–35 km, lo que indica que ese es el rango práctico en el que existe dependencia espacial entre las curvas TMP(t). El sill estimado (≈280–300) coincide bien con la variabilidad total observada, y el nugget es pequeño, lo que sugiere baja variabilidad a microescala y ausencia de discontinuidades. Aunque el empírico presenta oscilaciones para distancias grandes debido a la menor cantidad de pares disponibles, el modelo esférico proporciona una representación estable y coherente. En conjunto, el ajuste es satisfactorio y adecuado para ser utilizado en las etapas de kriging funcional.
5.4 Kriging Funcional
# ===============================================================
# 5.4 Kriging funcional con mejor manejo de variogramas
# ===============================================================
# Crear malla de predicción
buffer <- 2000
grd_pred <- expand.grid(
x = seq(min(coords_mat[,1]) - buffer, max(coords_mat[,1]) + buffer, length.out = 30),
y = seq(min(coords_mat[,2]) - buffer, max(coords_mat[,2]) + buffer, length.out = 30)
)
# Función auxiliar para ajustar variograma con valores iniciales
ajustar_vgm <- function(vg_emp, modelo = "Sph") {
# Valores iniciales basados en el empírico
sill_init <- median(vg_emp$gamma)
range_init <- median(vg_emp$dist)
nugget_init <- min(vg_emp$gamma) * 0.1
vg_modelo <- vgm(psill = sill_init, model = modelo,
range = range_init, nugget = nugget_init)
# Intentar ajuste
vg_fit <- tryCatch(
fit.variogram(vg_emp, vg_modelo, fit.sills = TRUE, fit.ranges = TRUE),
warning = function(w) {
# Si falla, usar modelo fijo basado en empírico
vgm(psill = max(vg_emp$gamma) * 0.8, model = modelo,
range = max(vg_emp$dist) * 0.5, nugget = min(vg_emp$gamma) * 0.2)
},
error = function(e) {
vgm(psill = max(vg_emp$gamma) * 0.8, model = modelo,
range = max(vg_emp$dist) * 0.5, nugget = min(vg_emp$gamma) * 0.2)
}
)
return(vg_fit)
}Método 1: Kriging por scores (FPCA)
# --- Método 1: Kriging por scores (FPCA) ---
fpca_tmp <- pca.fd(fd_tmp, nharm = 3)
cat("\nVarianza explicada por componentes:\n")
Varianza explicada por componentes:print(cumsum(fpca_tmp$varprop))[1] 0.8067429 0.9674037 0.9895827scores <- fpca_tmp$scores
# Kriging de cada score
pred_scores <- matrix(NA, nrow(grd_pred), ncol(scores))
for(k in 1:ncol(scores)) {
df_score <- data.frame(x = coords_mat[,1], y = coords_mat[,2], score = scores[,k])
coordinates(df_score) <- ~ x + y
vg_score <- variogram(score ~ 1, df_score)
vg_fit <- ajustar_vgm(vg_score, "Sph")
grd_sp <- grd_pred
coordinates(grd_sp) <- ~ x + y
gridded(grd_sp) <- TRUE
krig_score <- krige(score ~ 1, df_score, grd_sp, model = vg_fit)
pred_scores[, k] <- krig_score$var1.pred
cat("Score", k, "completado\n")
}[using ordinary kriging]
Score 1 completado
[using ordinary kriging]
Score 2 completado
[using ordinary kriging]
Score 3 completado# Reconstruir curvas predichas
coefs_pred_scores <- fpca_tmp$meanfd$coefs %*% t(rep(1, nrow(grd_pred))) +
fpca_tmp$harmonics$coefs %*% t(pred_scores)
fd_pred_scores <- fd(coefs_pred_scores, basis_fd)El análisis de Componentes Principales Funcionales (FPCA) mostró que la primera componente explica el 80.67% de la variabilidad de las curvas de temperatura TMP(t). Al incluir la segunda componente, la varianza acumulada asciende a 96.74%, mientras que la tercera lleva el total a 98.96%.
Esto indica que la mayor parte de la variación funcional se puede describir con solo dos componentes: la primera captura el patrón diurno general (mínima nocturna, ascenso durante el día, descenso vespertino), mientras que la segunda explica diferencias en la forma del ciclo diario, como el horario del máximo o la pendiente del calentamiento.
** Método 2: Kriging por lambda **
# --- Método 2: Kriging por lambda (coeficientes directos) ---
pred_coefs <- matrix(NA, nbasis, nrow(grd_pred))
for(b in 1:nbasis) {
df_coef <- data.frame(x = coords_mat[,1], y = coords_mat[,2], coef_val = C[b,])
coordinates(df_coef) <- ~ x + y
vg_coef <- variogram(coef_val ~ 1, df_coef)
vg_fit <- ajustar_vgm(vg_coef, "Sph")
grd_sp <- grd_pred
coordinates(grd_sp) <- ~ x + y
gridded(grd_sp) <- TRUE
krig_coef <- krige(coef_val ~ 1, df_coef, grd_sp, model = vg_fit)
pred_coefs[b, ] <- krig_coef$var1.pred
cat("Coeficiente", b, "de", nbasis, "completado\n")
}[using ordinary kriging]
Coeficiente 1 de 12 completado
[using ordinary kriging]
Coeficiente 2 de 12 completado
[using ordinary kriging]
Coeficiente 3 de 12 completado
[using ordinary kriging]
Coeficiente 4 de 12 completado
[using ordinary kriging]
Coeficiente 5 de 12 completado
[using ordinary kriging]
Coeficiente 6 de 12 completado
[using ordinary kriging]
Coeficiente 7 de 12 completado
[using ordinary kriging]
Coeficiente 8 de 12 completado
[using ordinary kriging]
Coeficiente 9 de 12 completado
[using ordinary kriging]
Coeficiente 10 de 12 completado
[using ordinary kriging]
Coeficiente 11 de 12 completado
[using ordinary kriging]
Coeficiente 12 de 12 completadofd_pred_lambda <- fd(pred_coefs, basis_fd)# ===============================================================
# 5.5 Visualización de predicciones
# ===============================================================
par(mfrow = c(2, 2))
# Punto 1
plot(fd_pred_scores[1], col = "blue", lwd = 2,
main = "Predicción punto 1",
xlab = "Hora", ylab = "Temp (°C)", ylim = c(10, 30))[1] "done"lines(fd_pred_lambda[1], col = "red", lwd = 2, lty = 2)
legend("topright", c("Scores", "Lambda"), col = c("blue", "red"), lty = c(1, 2), lwd = 2)
# Punto central
idx_centro <- nrow(grd_pred) %/% 2
plot(fd_pred_scores[idx_centro], col = "blue", lwd = 2,
main = paste("Predicción punto", idx_centro),
xlab = "Hora", ylab = "Temp (°C)", ylim = c(10, 30))[1] "done"lines(fd_pred_lambda[idx_centro], col = "red", lwd = 2, lty = 2)
# Media de predicciones
plot(mean.fd(fd_pred_scores), col = "blue", lwd = 2,
main = "Media de predicciones",
xlab = "Hora", ylab = "Temp (°C)")[1] "done"lines(mean.fd(fd_pred_lambda), col = "red", lwd = 2, lty = 2)
lines(media_fun, col = "black", lwd = 2, lty = 3)
legend("topright", c("Scores", "Lambda", "Original"),
col = c("blue", "red", "black"), lty = c(1, 2, 3), lwd = 2, cex = 0.8)
# Todas las curvas predichas
plot(fd_pred_scores, col = rgb(0, 0, 1, 0.3), lty = 1,
main = "Todas las curvas predichas (Scores)",
xlab = "Hora", ylab = "Temp (°C)")[1] "done"lines(mean.fd(fd_pred_scores), col = "blue", lwd = 3)
par(mfrow = c(1, 1))Se implementaron dos enfoques de predicción funcional espacial:
Kriging por scores (FPCA): Las curvas se representaron mediante los tres primeros componentes principales, que explican más del 98% de la varianza total. Se realizó kriging espacial sobre cada score y luego se reconstruyeron las curvas predichas.
Kriging por coeficientes B–spline (lambda): Cada curva funcional se representó directamente por sus coeficientes de base. Se realizó kriging independiente para cada coeficiente.
Al comparar las predicciones en diferentes puntos del dominio, se observó que ambos métodos producen curvas prácticamente iguales, sin diferencias sistemáticas en magnitud ni forma. La media de las predicciones coincide estrechamente con la media funcional original, lo que indica que el modelo captura adecuadamente la estructura espacial y temporal de la temperatura.
Las curvas predichas en todo el dominio muestran que la variación espacial es suave, coherente con el rango estimado en el semivariograma funcional (~30 km). No aparecen comportamientos anómalos ni discontinuidades, lo cual confirma que el ajuste del variograma y los métodos de kriging empleados son consistentes y apropiados para esta variable.
5.6 Cortes en tiempos específicos
# ===============================================================
# 5.6 Cortes en tiempos específicos (superficies espaciales)
# ===============================================================
library(ggplot2)
library(gridExtra)Warning: package 'gridExtra' was built under R version 4.5.2
Adjuntando el paquete: 'gridExtra'The following object is masked from 'package:dplyr':
combinehoras_corte <- c(6, 12, 18, 22)
plots_list <- list()
for(i in seq_along(horas_corte)) {
h <- horas_corte[i]
# Evaluar curvas predichas en hora h
vals_h <- eval.fd(h, fd_pred_scores)
df_plot <- data.frame(
coord_x = grd_pred$x,
coord_y = grd_pred$y,
temperatura = as.vector(vals_h)
)
p <- ggplot(df_plot, aes(x = coord_x, y = coord_y, fill = temperatura)) +
geom_tile() +
scale_fill_viridis_c(option = "plasma", name = "°C") +
coord_equal() +
labs(title = paste0("Temp a las ", h, ":00 h"),
x = "Easting (m)", y = "Northing (m)") +
theme_minimal() +
theme(legend.position = "right",
axis.text = element_text(size = 7))
plots_list[[i]] <- p
}
# Mostrar los 4 cortes
grid.arrange(grobs = plots_list, ncol = 2,
top = "Superficies de temperatura en diferentes horas")
A partir del kriging funcional aplicado a las curvas de temperatura TMP(t), se generaron superficies espaciales para diferentes horas del día (06:00, 12:00, 18:00 y 22:00). Estas superficies permiten visualizar cómo evoluciona la temperatura en el espacio a medida que avanza el día.
En general, se observa un comportamiento térmico coherente y consistente:
06:00 h: Temperaturas bajas (6–9 °C) y mayor contraste espacial. Corresponde al mínimo térmico típico de la madrugada.
12:00 h: Incremento notable (18–20 °C) con patrones más homogéneos debido al calentamiento solar.
18:00 h: Se alcanzan valores máximos (~22 °C), con ligeras variaciones espaciales.
22:00 h: Disminución de la temperatura (8–16 °C), recuperando cierta diferenciación espacial asociada al enfriamiento nocturno.
Estas superficies muestran que el modelo funcional es capaz de capturar adecuadamente la dinámica conjunta espacio–tiempo, reproduciendo tanto el ciclo diurno como las diferencias territoriales. Además, la continuidad espacial y temporal lograda por el kriging funcional refleja la suavidad inherente de los datos climáticos y la estructura de correlación espacial previamente modelada.
# ===============================================================
# 5.7 Mapa interactivo con Leaflet
# ===============================================================
library(leaflet)
library(raster)
# Crear raster para una hora específica (ej: 12:00)
vals_12 <- eval.fd(12, fd_pred_scores)
# Convertir a raster
ras_tmp <- rasterFromXYZ(data.frame(
x = grd_pred$x,
y = grd_pred$y,
z = as.vector(vals_12)
))
crs(ras_tmp) <- CRS(paste0("EPSG:", epsg_crs))
# Proyectar a WGS84
ras_tmp_ll <- projectRaster(ras_tmp, crs = CRS("+proj=longlat +datum=WGS84"))
# Mapa
pal <- colorNumeric(viridis::plasma(100), values(ras_tmp_ll), na.color = NA)
leaflet() %>%
addProviderTiles("CartoDB.Positron") %>%
addRasterImage(ras_tmp_ll, colors = pal, opacity = 0.7) %>%
leaflet::addLegend(pal = pal, values = values(ras_tmp_ll), title = "Temp 12:00 (°C)") %>%
leaflet::addCircleMarkers(
lng = st_coordinates(st_transform(st_as_sf(coords_tmp, coords = c("x","y"), crs = epsg_crs), 4326))[,1],
lat = st_coordinates(st_transform(st_as_sf(coords_tmp, coords = c("x","y"), crs = epsg_crs), 4326))[,2],
radius = 4, color = "red", fillOpacity = 0.8
)