Introducción

El departamento de Caldas pertenece a la región Andina Colombiana, su capital es Manizales. En este departamento es posible encontrar todos los pisos térmicos, desde los cálidos valles del río Magdalena y el río Cauca hasta las nieves perpetuas del Nevado del Ruiz.

Además, el departamento se encuentra completamente atravesado por las cordilleras andinas Central y Occidental en su totalidad. La topografía varía entre los 5.400 y 170 metros sobre el nivel del mar. El punto más alto es el Nevado del Ruiz, mientras que el municipio de La Dorada es el más bajo. Estas altitudes también generan una amplia diversidad en cuanto a clima y paisajes. Además de las grandes elevaciones como el páramo de Letras, el páramo de Sonson y el páramo de San Félix, también se pueden encontrar áreas llanas como el Valle interandino del Magdalena, el Valle del Risaralda, y cañones como el del Cauca.

Lo componen 27 municipios: Manizales, Aguadas, Anserma, Aranzazu, Belalcázar, Chinchiná, Filadelfia, La Dorada, La Merced, Manzanares, Marmato, Marquetalia, Marulanda, Neira, Norcasia, Pácora, Palestina, Pensilvania, Riosucio, Risaralda, Salamina, Samaná, San, José, Supía, Victoria, Villamaría, Viterbo.

Mapa_Caldas

Descripción de los datos de la velocidad del viento en el departamento de Caldas

A continuación, se analizan los datos de la velocidad del viento en el departamento de Caldas, modelando la dependencia espacio-temporal con fines predictivos.

Lectura de los datos

## Warning: package 'readxl' was built under R version 4.2.3

## [1] 162

Resumen de la Variable Velocidad del Viento y Hora

summary(Caldas$Hora)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    18.0    19.0    20.5    20.5    22.0    23.0

summary(Caldas$Velocidad_Viento)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.500   1.000   1.300   1.269   1.500   2.300

mean(Caldas$Velocidad_Viento)

## [1] 1.268519

Como se puede evidenciar el rango horario en el que se encuentran los datos es de 6pm a 11pm del día 12 de febrero del 2023, tiempo a través del cual se registró la velocidad del viento en los 27 municipios del Departamento de Caldas. La mínima velocidad del viento registrada en el presente conjunto de datos correspondió a 0.5 y la máxima a 2.3 medido en metros por segundo(m/s). Donde además la velocidad del viento promedio es de 1.268 m/s.

Comportamiento de la Velocidad del Viento por Hora

par(mfrow=c(1,1))
boxplot(Caldas$Velocidad_Viento~Caldas$Hora,col=c(2:7),
        main='Velocidad del viento por hora',
        ylab='Hora del día',xlab='Velocidad viento',cex=0.5,sub="Distribución de la variable velocidad del viento")
abline(h=mean(Caldas$Velocidad_Viento),col="blue",lty=1)

Cada gráfico de cajas indica la distribución de la variable velocidad del viento en un rango horario de 6pm a 11pm medido por hora,mostrando valores atípicos solo a las 7pm y un comportamiento inusual (en comparación a las demás horas) a las 6pm

boxplot(Caldas$Velocidad_Viento~Caldas$Municipio,col=c(2:7),
        main='Velocidad del viento por Munp.',cex.lab=0.2,
        ylab='Hora del día',xlab='Velocidad viento',cex=0.5,sub="Variable velocidad del viento  por municipio")
abline(h=mean(Caldas$Velocidad_Viento),col="blue",lty=1)

La distribución de la variable velocidad del viento ahora se evidencia espacialmente por municipio, por eso vemos en esta gràfica 27 boxplots, algunos parecen estar alejados de la media de los datos y muchos presentan valores atìpicos. Esto no fue posible verlo en el anterior gráfico

El comportamiento evidenciando en el grafico de boxplot correspondiente a la velocidad del viento por hora se puede corroborar por medio del siguiente gráfico:

par(mfrow=c(2,3))
datos18<-Caldas_utm[Caldas_utm$Hora %in% c(18),]
Este18<-datos18$Este
Norte18<-datos18$Norte
Velocidad_Viento18<-datos18$` Velocidad_Viento`
grillas18 <- interp(Este18,
                    Norte18,
                    Velocidad_Viento18)

persp(grillas18$x,
      grillas18$y,
      grillas18$z,
      xlab = "Este",
      ylab = "Norte",
      zlab = "Velocidad_Viento",
      phi = 30,
      theta = 20,
      col = 2,
      expand = .5,
      ticktype = "detailed",
      main='t=18')

datos19<-Caldas_utm[Caldas_utm$Hora %in% c(19),]
Este19<-datos19$Este
Norte19<-datos19$Norte
Velocidad_Viento19<-datos19$` Velocidad_Viento`
grillas19 <- interp(Este19,
                    Norte19,
                    Velocidad_Viento19)

persp(grillas19$x,
      grillas19$y,
      grillas19$z,
      xlab = "Este",
      ylab = "Norte",
      zlab = "Velocidad_Viento",
      phi = 30,
      theta = 20,
      col =3,
      expand = .5,
      ticktype = "detailed",
      main='t=19')

datos20<-Caldas_utm[Caldas_utm$Hora %in% c(20),]
Este20<-datos20$Este
Norte20<-datos20$Norte
Velocidad_Viento20<-datos20$` Velocidad_Viento`
grillas20 <- interp(Este20,
                    Norte20,
                    Velocidad_Viento20)

persp(grillas20$x,
      grillas20$y,
      grillas20$z,
      xlab = "Este",
      ylab = "Norte",
      zlab = "Velocidad_Viento",
      phi = 30,
      theta = 20,
      col = 4,
      expand = .5,
      ticktype = "detailed",
      main='t=20')

datos21<-Caldas_utm[Caldas_utm$Hora %in% c(21),]
Este21<-datos21$Este
Norte21<-datos21$Norte
Velocidad_Viento21<-datos21$` Velocidad_Viento`
grillas21 <- interp(Este21,
                    Norte21,
                    Velocidad_Viento21)

persp(grillas21$x,
      grillas21$y,
      grillas21$z,
      xlab = "Este",
      ylab = "Norte",
      zlab = "Velocidad_Viento",
      phi = 30,
      theta = 20,
      col = 5,
      expand = .5,
      ticktype = "detailed",
      main='t=21')

datos22<-Caldas_utm[Caldas_utm$Hora %in% c(22),]
Este22<-datos22$Este
Norte22<-datos22$Norte
Velocidad_Viento22<-datos22$` Velocidad_Viento`
grillas22 <- interp(Este22,
                    Norte22,
                    Velocidad_Viento22)

persp(grillas22$x,
      grillas22$y,
      grillas22$z,
      xlab = "Este",
      ylab = "Norte",
      zlab = "Velocidad_Viento",
      phi = 30,
      theta = 20,
      col = 6,
      expand = .5,
      ticktype = "detailed",
      main='t=22')

datos23<-Caldas_utm[Caldas_utm$Hora %in% c(23),]
Este23<-datos23$Este
Norte23<-datos23$Norte
Velocidad_Viento23<-datos23$` Velocidad_Viento`
grillas23 <- interp(Este23,
                    Norte23,
                    Velocidad_Viento23)

persp(grillas23$x,
      grillas23$y,
      grillas23$z,
      xlab = "Este",
      ylab = "Norte",
      zlab = "Velocidad_Viento",
      phi = 30,
      theta = 20,
      col = 7,
      expand = .5,
      ticktype = "detailed",
      main='t=23')

La velocidad del viento a las 6pm tiene un comportamiento un poco inusual en comparación con las demás horas mostradas las cuales si tienen comportamientos bastante similares (esto se verá reflejado en los semivariogramas empíricos marginales).

Debido a que el análisis espacio tiempo en este caso se hace sobre coordenadas planas, es importante hacer la conversión a coordenadas UTM, dado que los datos datos presentan coordenadas de longitud y latitud.

Conversión a coordenadas UTM

deci_coord3 = SpatialPoints(cbind(Caldas$Longitud, Caldas$Latitud),
                            proj4string = CRS("+proj=longlat"))
Caldas_geod <- data.frame(deci_coord3, Caldas$Velocidad_Viento,Caldas$Fecha_Hora,Caldas$Hora)
head(Caldas_geod)

##   coords.x1 coords.x2 Caldas.Velocidad_Viento   Caldas.Fecha_Hora Caldas.Hora
## 1 -75.45487  5.610248                     1.0 1899-12-31 18:00:00          18
## 2 -75.45487  5.610248                     1.4 1899-12-31 19:00:00          19
## 3 -75.45487  5.610248                     1.3 1899-12-31 20:00:00          20
## 4 -75.45487  5.610248                     1.3 1899-12-31 21:00:00          21
## 5 -75.45487  5.610248                     1.2 1899-12-31 22:00:00          22
## 6 -75.45487  5.610248                     1.2 1899-12-31 23:00:00          23

# UTM
library(rgdal)
CRS_UTM_NY3 = (c("+proj=utm +zone=18 +datum=NAD83 
                  +units=m +no_defs"))
utm_coord3 = spTransform(deci_coord3, CRS(CRS_UTM_NY3))
utm_coord_df3 = as.data.frame(utm_coord3)
Caldas_utm <- data.frame(utm_coord3,  Caldas$Velocidad_Viento,Caldas$Fecha_Hora,Caldas$Hora)
colnames(Caldas_utm)<-c('Este','Norte',' Velocidad_Viento','Fecha_Hora','Hora')
head(Caldas_utm)

##       Este    Norte  Velocidad_Viento          Fecha_Hora Hora
## 1 449624.5 620140.4               1.0 1899-12-31 18:00:00   18
## 2 449624.5 620140.4               1.4 1899-12-31 19:00:00   19
## 3 449624.5 620140.4               1.3 1899-12-31 20:00:00   20
## 4 449624.5 620140.4               1.3 1899-12-31 21:00:00   21
## 5 449624.5 620140.4               1.2 1899-12-31 22:00:00   22
## 6 449624.5 620140.4               1.2 1899-12-31 23:00:00   23

Caldasg <- as.geodata(Caldas_utm)

## as.geodata: 135 replicated data locations found. 
##  Consider using jitterDupCoords() for jittering replicated locations. 
## WARNING: there are data at coincident or very closed locations, some of the geoR's functions may not work.
##  Use function dup.coords() to locate duplicated coordinates.
##  Consider using jitterDupCoords() for jittering replicated locations

Al igual que en el Kriging espacial, un paso importante antes de llegar al Kriging espacio tiempo es la modelación de la media, la cual va incluir no sólo las coordenadas Este y Norte, sino además el tiempo.

Para ello es necesario primero explorar cómo se comporta la velocidad del viento con respecto no sólo a las coordenadas sino también con respecto al tiempo, ya que de esta maneta tendremos una idea del tipo de tendencia presente.

Explorando tendencia

Explorando tendencia con coordenada Este, Velocidad del viento y Hora

par(mfrow=c(1,1))
scatterplot3d(Caldas_utm[,c(1,3,4)], highlight.3d=TRUE, col.axis="blue",
              col.grid="lightblue", main="Tendencia", pch=20,xlab='Este',
              ylab='Hora',zlab='Velocidad del viento',sub='Tendenciacon respecto a Este')

Explorando tendencia con coordenada Este, Velocidad del viento y Hora

Explorando tendencia con coordenadas Este y Norte

scatterplot3d(Caldas_utm[,c(1,2,3)], highlight.3d=TRUE, col.axis="blue",xlab='Este',
              ylab='Velocidad del viento',zlab='Norte',
              col.grid="lightblue",
              main="Tendencia", pch=20,sub='Tendencia con respecto a ambas coordenadas')

Los gráficos parecen mostrar un comportamiento cíclico con respecto a cada una de las coordenadas sin embargo no es muy claro.

A continuación se proponen 10 modelos de regresión, a través de los cuales se explorará el comportamieto de los residualespara determinar qué tanto nos permite cada uno corregir la tendencia de los datos.

Modelando y corrigiendo la tendencia

Velocidad_Viento2<-Caldas_utm[,3]
Este2<-Caldas_utm[,1]
Norte2<-Caldas_utm[,2]
Hora2<-Caldas_utm[,5]
Fecha_Hora<-Caldas_utm[,4]

Modelo 1

mod1<-lm(Velocidad_Viento2~I(Este2)+I(Norte2)+I(Hora2))
summary(mod1)

## 
## Call:
## lm(formula = Velocidad_Viento2 ~ I(Este2) + I(Norte2) + I(Hora2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.84757 -0.28318  0.01159  0.25639  0.84400 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -5.015e-02  9.494e-01  -0.053  0.95794   
## I(Este2)     2.510e-06  9.905e-07   2.534  0.01226 * 
## I(Norte2)   -1.759e-06  1.709e-06  -1.029  0.30519   
## I(Hora2)     5.915e-02  1.765e-02   3.351  0.00101 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3837 on 158 degrees of freedom
## Multiple R-squared:  0.1007, Adjusted R-squared:  0.08362 
## F-statistic: 5.897 on 3 and 158 DF,  p-value: 0.0007688

par(mfrow=c(2,2))
plot(Este2,residuals(mod1),pch=16,col=4,ylab="residuales",
     main="Este vs residuales")
plot(Norte2,residuals(mod1),pch=16,col=4,ylab="residuales",
     main="Norte vs residuales")
plot(Hora2, residuals(mod1),pch=16,col=4,ylab="residuales",
     main="Hora vs residuales")
hist(residuals(mod1),col=4,main="Residuales" )

Modelo 2

mod2<-lm(Velocidad_Viento2~I(Este2)+I(Norte2)+I(Hora2)+I(Este2*Norte2))
summary(mod2)

## 
## Call:
## lm(formula = Velocidad_Viento2 ~ I(Este2) + I(Norte2) + I(Hora2) + 
##     I(Este2 * Norte2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.95272 -0.29947  0.03438  0.25596  0.87565 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -5.528e+01  1.403e+01  -3.941 0.000122 ***
## I(Este2)           1.285e-04  3.195e-05   4.023 8.92e-05 ***
## I(Norte2)          9.156e-05  2.371e-05   3.862 0.000164 ***
## I(Hora2)           5.915e-02  1.689e-02   3.502 0.000601 ***
## I(Este2 * Norte2) -2.128e-10  5.392e-11  -3.946 0.000120 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3672 on 157 degrees of freedom
## Multiple R-squared:  0.1818, Adjusted R-squared:  0.161 
## F-statistic: 8.723 on 4 and 157 DF,  p-value: 2.197e-06

par(mfrow=c(2,2))
plot(Este2,residuals(mod2),pch=16,col=4,ylab="residuales",
     main="Este vs residuales")
plot(Norte2,residuals(mod2),pch=16,col=4,ylab="residuales",
     main="Norte vs residuales")
plot(Hora2, residuals(mod2),pch=16,col=4,ylab="residuales",
     main="Hora vs residuales")
hist(residuals(mod2),col=4,main="Residuales" )

Modelo 3

mod3<-lm(Velocidad_Viento2~I(Este2)+I(Norte2)+I(Hora2)+
           I(Este2*Norte2)+I(Este2*Norte2*Hora2))
summary(mod3)

## 
## Call:
## lm(formula = Velocidad_Viento2 ~ I(Este2) + I(Norte2) + I(Hora2) + 
##     I(Este2 * Norte2) + I(Este2 * Norte2 * Hora2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.96177 -0.28002  0.02071  0.24932  0.78972 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -4.659e+01  1.422e+01  -3.275 0.001301 ** 
## I(Este2)                   1.285e-04  3.143e-05   4.090 6.90e-05 ***
## I(Norte2)                  9.156e-05  2.332e-05   3.927 0.000129 ***
## I(Hora2)                  -3.647e-01  1.700e-01  -2.146 0.033450 *  
## I(Este2 * Norte2)         -2.457e-10  5.464e-11  -4.496 1.34e-05 ***
## I(Este2 * Norte2 * Hora2)  1.605e-12  6.407e-13   2.506 0.013248 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3611 on 156 degrees of freedom
## Multiple R-squared:  0.2135, Adjusted R-squared:  0.1883 
## F-statistic: 8.469 on 5 and 156 DF,  p-value: 4.109e-07

par(mfrow=c(2,2))
plot(Este2,residuals(mod3),pch=16,col=4,ylab="residuales",
     main="Este vs residuales")
plot(Norte2,residuals(mod3),pch=16,col=4,ylab="residuales",
     main="Norte vs residuales")
plot(Hora2, residuals(mod3),pch=16,col=4,ylab="residuales",
     main="Hora vs residuales")
hist(residuals(mod3),col=4,main="Residuales" )

Modelo 4

mod4<-lm(Velocidad_Viento2~I(Este2)+I(Hora2)+I(Este2*Norte2*Hora2)+
       I(Hora2^2)+I(Este2^2))
summary(mod4)

## 
## Call:
## lm(formula = Velocidad_Viento2 ~ I(Este2) + I(Hora2) + I(Este2 * 
##     Norte2 * Hora2) + I(Hora2^2) + I(Este2^2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1114 -0.2198 -0.0048  0.2339  0.8634 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -4.357e+01  6.728e+00  -6.477 1.16e-09 ***
## I(Este2)                   9.040e-05  2.080e-05   4.346 2.49e-05 ***
## I(Hora2)                   2.239e+00  4.476e-01   5.002 1.51e-06 ***
## I(Este2 * Norte2 * Hora2) -1.846e-13  1.651e-13  -1.118    0.265    
## I(Hora2^2)                -5.198e-02  1.086e-02  -4.787 3.90e-06 ***
## I(Este2^2)                -9.207e-11  2.191e-11  -4.202 4.43e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3448 on 156 degrees of freedom
## Multiple R-squared:  0.2832, Adjusted R-squared:  0.2602 
## F-statistic: 12.33 on 5 and 156 DF,  p-value: 4.41e-10

par(mfrow=c(2,2))
plot(Este2,residuals(mod4),pch=16,col=4,ylab="residuales",
     main="Este vs residuales")
plot(Norte2,residuals(mod4),pch=16,col=4,ylab="residuales",
     main="Norte vs residuales")
plot(Hora2, residuals(mod4),pch=16,col=4,ylab="residuales",
     main="Hora vs residuales")
hist(residuals(mod4),col=4,main="Residuales" )

Modelo 5

mod5<-lm(Velocidad_Viento2~I(Este2)+I(Hora2)+I(Este2*Norte2*Hora2)
         +I(Este2*Hora2))
summary(mod5)

## 
## Call:
## lm(formula = Velocidad_Viento2 ~ I(Este2) + I(Hora2) + I(Este2 * 
##     Norte2 * Hora2) + I(Este2 * Hora2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.98566 -0.27385  0.00158  0.26382  0.81489 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                1.180e+01  4.638e+00   2.546  0.01188 * 
## I(Este2)                  -2.602e-05  1.024e-05  -2.541  0.01203 * 
## I(Hora2)                  -5.733e-01  2.257e-01  -2.540  0.01206 * 
## I(Este2 * Norte2 * Hora2) -2.532e-13  1.848e-13  -1.370  0.17265   
## I(Este2 * Hora2)           1.549e-06  5.151e-07   3.007  0.00307 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3751 on 157 degrees of freedom
## Multiple R-squared:  0.146,  Adjusted R-squared:  0.1242 
## F-statistic: 6.708 on 4 and 157 DF,  p-value: 5.206e-05

par(mfrow=c(2,2))
plot(Este2,residuals(mod5),pch=16,col=4,ylab="residuales",
     main="Este vs residuales")
plot(Norte2,residuals(mod5),pch=16,col=4,ylab="residuales",
     main="Norte vs residuales")
plot(Hora2, residuals(mod5),pch=16,col=4,ylab="residuales",
     main="Hora vs residuales")
hist(residuals(mod5),col=4,main="Residuales" )

Modelo 6

mod6<-lm(Velocidad_Viento2~Este2+Norte2+I(Este2*Norte2)+Hora2+
           I(Hora2^2)+I(sin(Hora2))+I(Hora2*Norte2)+I(Hora2*Este2)+
           I(Hora2*Norte2^2)+I(Hora2*Este2^2)+I(cos(Hora2)))
summary(mod6)

## 
## Call:
## lm(formula = Velocidad_Viento2 ~ Este2 + Norte2 + I(Este2 * Norte2) + 
##     Hora2 + I(Hora2^2) + I(sin(Hora2)) + I(Hora2 * Norte2) + 
##     I(Hora2 * Este2) + I(Hora2 * Norte2^2) + I(Hora2 * Este2^2) + 
##     I(cos(Hora2)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0696 -0.1449  0.0190  0.1852  0.7584 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)   
## (Intercept)         -6.973e+01  3.597e+01  -1.938  0.05445 . 
## Este2                2.834e-05  5.485e-05   0.517  0.60605   
## Norte2               5.129e-05  4.357e-05   1.177  0.24094   
## I(Este2 * Norte2)   -9.607e-11  9.101e-11  -1.056  0.29282   
## Hora2                4.206e+00  3.044e+00   1.382  0.16912   
## I(Hora2^2)          -1.223e-01  6.136e-02  -1.994  0.04802 * 
## I(sin(Hora2))       -2.446e-01  2.227e-01  -1.098  0.27378   
## I(Hora2 * Norte2)   -1.192e-06  5.123e-06  -0.233  0.81628   
## I(Hora2 * Este2)     4.818e-06  1.609e-06   2.994  0.00322 **
## I(Hora2 * Norte2^2)  5.312e-13  4.335e-12   0.123  0.90263   
## I(Hora2 * Este2^2)  -3.536e-12  1.645e-12  -2.149  0.03322 * 
## I(cos(Hora2))        1.989e-01  7.294e-02   2.727  0.00715 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3253 on 150 degrees of freedom
## Multiple R-squared:  0.3864, Adjusted R-squared:  0.3414 
## F-statistic: 8.586 on 11 and 150 DF,  p-value: 1.164e-11

par(mfrow=c(2,2))
plot(Este2,residuals(mod6),pch=16,col=4,ylab="residuales",
     main="Este vs residuales")
plot(Norte2,residuals(mod6),pch=16,col=4,ylab="residuales",
     main="Norte vs residuales")
plot(Hora2, residuals(mod6),pch=16,col=4,ylab="residuales",
     main="Hora vs residuales")
hist(residuals(mod6),col=4,main="Residuales" )

Modelo 7

mod7<-lm(Velocidad_Viento2~I(cos(Este2))+I(sin(Norte2))+Hora2)
summary(mod7)

## 
## Call:
## lm(formula = Velocidad_Viento2 ~ I(cos(Este2)) + I(sin(Norte2)) + 
##     Hora2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.84380 -0.29222  0.00458  0.23637  0.97407 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)   
## (Intercept)     0.03184    0.36479   0.087  0.93055   
## I(cos(Este2))   0.02369    0.04221   0.561  0.57539   
## I(sin(Norte2)) -0.11686    0.05213  -2.242  0.02638 * 
## Hora2           0.05915    0.01772   3.338  0.00105 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3852 on 158 degrees of freedom
## Multiple R-squared:  0.09356,    Adjusted R-squared:  0.07635 
## F-statistic: 5.436 on 3 and 158 DF,  p-value: 0.001389

par(mfrow=c(2,2))
plot(Este2,residuals(mod7),pch=16,col=4,ylab="residuales",
     main="Este vs residuales")
plot(Norte2,residuals(mod7),pch=16,col=4,ylab="residuales",
     main="Norte vs residuales")
plot(Hora2, residuals(mod7),pch=16,col=4,ylab="residuales",
     main="Hora vs residuales")
hist(residuals(mod7),col=4,main="Residuales" )

Modelo 8

mod8<-lm(Velocidad_Viento2~I(sin(Norte2))+Hora2+I(Hora2^2)+
           I(sin(Norte2)*Hora2)+
           I(sin(Norte2)^2*Hora2+I(sin(Norte2)*Hora2^2)))
par(mfrow=c(2,2))
plot(Este2,residuals(mod8),pch=16,col=4,ylab="residuales",
     main="Este vs residuales")
plot(Norte2,residuals(mod8),pch=16,col=4,ylab="residuales",
     main="Norte vs residuales")
plot(Hora2, residuals(mod8),pch=16,col=4,ylab="residuales",
     main="Hora vs residuales")
hist(residuals(mod8),col=4,main="Residuales" )

Modelo 9

mod9<-lm(Velocidad_Viento2~I(sin(Norte2))+I(Hora2^2)+I(sin(Norte2)^2*Hora2^2))
summary(mod9)

## 
## Call:
## lm(formula = Velocidad_Viento2 ~ I(sin(Norte2)) + I(Hora2^2) + 
##     I(sin(Norte2)^2 * Hora2^2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.87186 -0.27724 -0.00645  0.24127  0.97769 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 0.6681885  0.1857658   3.597  0.00043 ***
## I(sin(Norte2))             -0.1113875  0.0521409  -2.136  0.03420 *  
## I(Hora2^2)                  0.0012987  0.0004404   2.949  0.00368 ** 
## I(sin(Norte2)^2 * Hora2^2)  0.0002083  0.0002250   0.925  0.35615    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3858 on 158 degrees of freedom
## Multiple R-squared:  0.0909, Adjusted R-squared:  0.07364 
## F-statistic: 5.266 on 3 and 158 DF,  p-value: 0.00173

par(mfrow=c(2,2))
plot(Este2,residuals(mod9),pch=16,col=4,ylab="residuales",
     main="Este vs residuales")
plot(Norte2,residuals(mod9),pch=16,col=4,ylab="residuales",
     main="Norte vs residuales")
plot(Hora2, residuals(mod9),pch=16,col=4,ylab="residuales",
     main="Hora vs residuales")
hist(residuals(mod9),col=4,main="Residuales" )

Modelo 10

mod10<-lm(Velocidad_Viento2 ~ I(Este2) + I(Norte2) + I(Hora2) +
            I(sin(2*pi*Hora2/24)) +I(cos(2*pi*Hora2/24))+I(Este2^2))
summary(mod10)

## 
## Call:
## lm(formula = Velocidad_Viento2 ~ I(Este2) + I(Norte2) + I(Hora2) + 
##     I(sin(2 * pi * Hora2/24)) + I(cos(2 * pi * Hora2/24)) + I(Este2^2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.07536 -0.18982  0.00006  0.21411  0.88197 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -5.586e+01  1.477e+01  -3.781 0.000222 ***
## I(Este2)                   9.207e-05  2.018e-05   4.564 1.02e-05 ***
## I(Norte2)                 -2.842e-06  1.522e-06  -1.867 0.063755 .  
## I(Hora2)                   1.730e+00  6.848e-01   2.527 0.012511 *  
## I(sin(2 * pi * Hora2/24)) -5.367e+00  1.707e+00  -3.145 0.001990 ** 
## I(cos(2 * pi * Hora2/24)) -4.410e+00  2.208e+00  -1.997 0.047529 *  
## I(Este2^2)                -9.592e-11  2.159e-11  -4.443 1.68e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3372 on 155 degrees of freedom
## Multiple R-squared:  0.3186, Adjusted R-squared:  0.2922 
## F-statistic: 12.08 on 6 and 155 DF,  p-value: 4.117e-11

par(mfrow=c(2,2))
plot(Este2,residuals(mod10),pch=16,col=4,ylab="residuales",
     main="Este vs residuales")
plot(Norte2,residuals(mod10),pch=16,col=4,ylab="residuales",
     main="Norte vs residuales")
plot(Hora2, residuals(mod10),pch=16,col=4,ylab="residuales",
     main="Hora vs residuales")
hist(residuals(mod10),col=4,main="Residuales" )

Finalmente escogemos el modelo 10 al considerar que corrige mejor la tendencia dado los gráficos presentados.

Conforme a lo anterior se procede calcular el semivariograma empírico marginal fijando un tiempo específico en todas las ubicaciones espaciales disponibles. Para ello se utiliza el estimador Matheron por medio de la funciòn variog de R.

Semivariograma empírico marginal (sin tendencia)**

Para t=18 (6pm)

Caldas_utm$res<-residuals(mod10)
Caldasg18<-as.geodata(Caldas_utm[Hora==18,])
x18<-(Caldas_utm[Hora==18,])$Este
y18<-(Caldas_utm[Hora==18,])$Norte
Hora18<-(Caldas_utm[Hora==18,])$Hora
v18<-variog(Caldasg18,trend=~x18+y18+Hora18+sin(2*pi*Hora18/24)+
              cos(2*pi*Hora18/24)+x18^2)

## variog: computing omnidirectional variogram

head(v18)

## $u
##  [1]   5395.709  16187.128  26978.547  37769.966  48561.385  59352.803
##  [7]  70144.222  80935.641  91727.060 102518.479 113309.897 124101.316
## [13] 134892.735
## 
## $v
##  [1] 0.00650725 0.01025879 0.03488947 0.03424451 0.01537295 0.02663296
##  [7] 0.03501000 0.02252215 0.04519366 0.06005188 0.05171854 0.02277749
## [13] 0.03284614
## 
## $n
##  [1] 13 40 50 56 55 41 33 21 14 15  6  4  2
## 
## $sd
##  [1] 0.007467097 0.019626239 0.036607486 0.037797890 0.019551864 0.034963646
##  [7] 0.038258757 0.027273783 0.047445571 0.043747878 0.032340939 0.033818872
## [13] 0.019508041
## 
## $bins.lim
##  [1] 1.000000e-12 1.079142e+04 2.158284e+04 3.237426e+04 4.316568e+04
##  [6] 5.395709e+04 6.474851e+04 7.553993e+04 8.633135e+04 9.712277e+04
## [11] 1.079142e+05 1.187056e+05 1.294970e+05 1.402884e+05
## 
## $ind.bin
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Para t=19

Caldasg19<-as.geodata(Caldas_utm[Hora==19,])
x19<-(Caldas_utm[Hora==19,])$Este
y19<-(Caldas_utm[Hora==19,])$Norte
Hora19<-(Caldas_utm[Hora==19,])$Hora
v19<-variog(Caldasg19,trend=~x19+y19+Hora19+sin(2*pi*Hora19/24)+
              cos(2*pi*Hora19/24)+x19^2)

## variog: computing omnidirectional variogram

head(v19)

## $u
##  [1]   5395.709  16187.128  26978.547  37769.966  48561.385  59352.803
##  [7]  70144.222  80935.641  91727.060 102518.479 113309.897 124101.316
## [13] 134892.735
## 
## $v
##  [1] 0.01180141 0.06040763 0.10942160 0.20422567 0.14404455 0.06794910
##  [7] 0.04837115 0.04529205 0.09917010 0.14436957 0.04875351 0.05023684
## [13] 0.01193803
## 
## $n
##  [1] 13 40 50 56 55 41 33 21 14 15  6  4  2
## 
## $sd
##  [1] 0.01607683 0.07787520 0.17260439 0.22793512 0.13444139 0.09920299
##  [7] 0.04817132 0.04771789 0.08395762 0.13932496 0.02142130 0.04526973
## [13] 0.01055133
## 
## $bins.lim
##  [1] 1.000000e-12 1.079142e+04 2.158284e+04 3.237426e+04 4.316568e+04
##  [6] 5.395709e+04 6.474851e+04 7.553993e+04 8.633135e+04 9.712277e+04
## [11] 1.079142e+05 1.187056e+05 1.294970e+05 1.402884e+05
## 
## $ind.bin
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Para t=20

Caldasg20<-as.geodata(Caldas_utm[Hora==20,])
x20<-(Caldas_utm[Hora==20,])$Este
y20<-(Caldas_utm[Hora==20,])$Norte
Hora20<-(Caldas_utm[Hora==20,])$Hora
v20<-variog(Caldasg20,trend=~x20+y20+Hora20+sin(2*pi*Hora20/24)+
              cos(2*pi*Hora20/24)+x20^2)

## variog: computing omnidirectional variogram

head(v20)

## $u
##  [1]   5395.709  16187.128  26978.547  37769.966  48561.385  59352.803
##  [7]  70144.222  80935.641  91727.060 102518.479 113309.897 124101.316
## [13] 134892.735
## 
## $v
##  [1] 0.017050968 0.077610521 0.137691348 0.248295128 0.171301764 0.078813950
##  [7] 0.048700161 0.064778685 0.184000582 0.266530691 0.122400699 0.104854677
## [13] 0.001765365
## 
## $n
##  [1] 13 40 50 56 55 41 33 21 14 15  6  4  2
## 
## $sd
##  [1] 0.026804758 0.098620725 0.192766218 0.264368846 0.149121400 0.106195065
##  [7] 0.047308051 0.077648828 0.143127920 0.249772817 0.041272331 0.078591353
## [13] 0.001654101
## 
## $bins.lim
##  [1] 1.000000e-12 1.079142e+04 2.158284e+04 3.237426e+04 4.316568e+04
##  [6] 5.395709e+04 6.474851e+04 7.553993e+04 8.633135e+04 9.712277e+04
## [11] 1.079142e+05 1.187056e+05 1.294970e+05 1.402884e+05
## 
## $ind.bin
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Para t=21

Caldasg21<-as.geodata(Caldas_utm[Hora==21,])
x21<-(Caldas_utm[Hora==21,])$Este
y21<-(Caldas_utm[Hora==21,])$Norte
Hora21<-(Caldas_utm[Hora==21,])$Hora
v21<-variog(Caldasg21,trend=~x21+y21+Hora21+sin(2*pi*Hora21/24)+
              cos(2*pi*Hora21/24)+x21^2)

## variog: computing omnidirectional variogram

head(v21)

## $u
##  [1]   5395.709  16187.128  26978.547  37769.966  48561.385  59352.803
##  [7]  70144.222  80935.641  91727.060 102518.479 113309.897 124101.316
## [13] 134892.735
## 
## $v
##  [1] 0.021790613 0.097909129 0.157601185 0.285205425 0.197546239 0.092932244
##  [7] 0.067519091 0.080070848 0.212497273 0.292965471 0.124341117 0.131294777
## [13] 0.001695084
## 
## $n
##  [1] 13 40 50 56 55 41 33 21 14 15  6  4  2
## 
## $sd
##  [1] 0.029119589 0.121148703 0.237046785 0.318156049 0.170976661 0.125125282
##  [7] 0.069732924 0.104952305 0.186696802 0.286170715 0.048134032 0.102699002
## [13] 0.001276805
## 
## $bins.lim
##  [1] 1.000000e-12 1.079142e+04 2.158284e+04 3.237426e+04 4.316568e+04
##  [6] 5.395709e+04 6.474851e+04 7.553993e+04 8.633135e+04 9.712277e+04
## [11] 1.079142e+05 1.187056e+05 1.294970e+05 1.402884e+05
## 
## $ind.bin
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Para t=22

Caldasg22<-as.geodata(Caldas_utm[Hora==22,])
x22<-(Caldas_utm[Hora==22,])$Este
y22<-(Caldas_utm[Hora==22,])$Norte
Hora22<-(Caldas_utm[Hora==22,])$Hora
v22<-variog(Caldasg22,trend=~x22+y22+Hora22+sin(2*pi*Hora22/24)+
              cos(2*pi*Hora22/24)+x22^2)

## variog: computing omnidirectional variogram

head(v22)

## $u
##  [1]   5395.709  16187.128  26978.547  37769.966  48561.385  59352.803
##  [7]  70144.222  80935.641  91727.060 102518.479 113309.897 124101.316
## [13] 134892.735
## 
## $v
##  [1] 0.01549740 0.08692558 0.14484948 0.24541773 0.15724051 0.08371613
##  [7] 0.05558598 0.10541218 0.33432641 0.43705481 0.28308769 0.17627272
## [13] 0.05496025
## 
## $n
##  [1] 13 40 50 56 55 41 33 21 14 15  6  4  2
## 
## $sd
##  [1] 0.024101324 0.100387450 0.198541759 0.294691894 0.148808380 0.092873934
##  [7] 0.073334977 0.128384971 0.267309222 0.413556293 0.158641833 0.127664086
## [13] 0.005772763
## 
## $bins.lim
##  [1] 1.000000e-12 1.079142e+04 2.158284e+04 3.237426e+04 4.316568e+04
##  [6] 5.395709e+04 6.474851e+04 7.553993e+04 8.633135e+04 9.712277e+04
## [11] 1.079142e+05 1.187056e+05 1.294970e+05 1.402884e+05
## 
## $ind.bin
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Para t=23

Caldasg23<-as.geodata(Caldas_utm[Hora==23,])
x23<-(Caldas_utm[Hora==23,])$Este
y23<-(Caldas_utm[Hora==23,])$Norte
Hora23<-(Caldas_utm[Hora==23,])$Hora
v23<-variog(Caldasg23,trend=~x23+y23+Hora23+sin(2*pi*Hora23/24)+
              cos(2*pi*Hora23/24)+x23^2)

## variog: computing omnidirectional variogram

head(v23)

## $u
##  [1]   5395.709  16187.128  26978.547  37769.966  48561.385  59352.803
##  [7]  70144.222  80935.641  91727.060 102518.479 113309.897 124101.316
## [13] 134892.735
## 
## $v
##  [1] 0.009185506 0.057488375 0.114980753 0.174931276 0.101096138 0.066198812
##  [7] 0.056391687 0.104221095 0.357286591 0.409387837 0.302186596 0.128804610
## [13] 0.125774648
## 
## $n
##  [1] 13 40 50 56 55 41 33 21 14 15  6  4  2
## 
## $sd
##  [1] 0.012021705 0.070341700 0.143794128 0.198774997 0.136087321 0.099963150
##  [7] 0.053520229 0.107814038 0.262895169 0.402337191 0.203855185 0.126309082
## [13] 0.003426758
## 
## $bins.lim
##  [1] 1.000000e-12 1.079142e+04 2.158284e+04 3.237426e+04 4.316568e+04
##  [6] 5.395709e+04 6.474851e+04 7.553993e+04 8.633135e+04 9.712277e+04
## [11] 1.079142e+05 1.187056e+05 1.294970e+05 1.402884e+05
## 
## $ind.bin
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Comparación de los semivariogramas empíricos marginales

par(mfrow=c(3,2))
plot(v18,pch=16,col=2,main="semivariograma t=6pm")
plot(v19,pch=16,col=3,main="semivariograma t=7pm")
plot(v20,pch=16,col=4,main="semivariograma t=8pm")
plot(v21,pch=16,col=5,main="semivariograma t=9pm")
plot(v22,pch=16,col=6,main="semivariograma t=10pm")
plot(v23,pch=16,col=7,main="semivariograma t=11pm")

Vemos que así como en la parte exploratoria de los datos ya presentada se mostraba un comportamiento inusual en la velocidad del viento a las 6pm, el semivariograma marginal correspondiente a las 6pm también muestra un comportamiento particular con respecto a los demás semivariogramas empíricos marginales.

Kriging espacio tiempo

Veamos los rezagos espaciales presente en los datos así como los rezagos temporales:

Matriz de distancias

x1 <- Caldas_utm$Este
x2 <- Caldas_utm$Norte
t <- Caldas_utm$Hora
grillaSpT=cbind(x1,x2,t)
#matriz de distancias (rezagos) espaciales
matDistSp=as.matrix(dist(grillaSpT[,1:2]))
#matriz de distancias (rezagos) temporales
matDistT=as.matrix(dist(grillaSpT[,3:3]))

Una vez calculadas las matrices de rezagos tango espaciales como temporales es importante explorar modelos de semivariograma que se ajusten lo mejor posible al comportamiento anteriomente mostrado (aunque marginalmente).Esta exloración se hace un poco a ciegas debido a las dificultades que presenta en este sentido el kriging espacio tiempo.

Para esta exploración escogimos el modelo Cressie 1 y el modelo Cressie Huang para determinar por medio del MSE después de la correspondiente optimización cuál modelo sería el más adecuado.

Kriging espacio tiempo con el modelo Cressie1

Definición de la función Cressie1

library(SimDesign)

## Warning: package 'SimDesign' was built under R version 4.2.3

cressie1=function(h,u,p){(p[1]^2/((p[2]^2*u^2+1)))*
    exp(-(p[3]^2*h^2)/(p[2]^2*u^2+1))}
sim1=residuals(mod10)
datos1=cbind(grillaSpT,sim1)
names(datos1)=c("x","y","t","z((x,y),t)")

Planteamiento de la Log-Verosimilitud

LV<-function(p,h,u,modelo,z){0.5*(length(z)*log(2*pi)+
    log(det(modelo(p,h=matDistSp,u=matDistT)))+
      t(z)%*%solve(modelo(p,h=matDistSp,u=matDistT))%*%z)}

Encontrando los parámetros óptimos

p0 <- c(0.8,0.5,0.9) 
lo <- c(1,1,1)    
estima=optim(p0,LV,h=matDistSp,u=matDistT,modelo=cressie1,z=sim1,
             hessian = T,lower = lo,upper = c(Inf,Inf,Inf),
             method = "L-BFGS-B")
# se le puede dar la opcion de la matriz hessiana

Cálculo del semivariograma con los parámetros óptimos

sigma1=cressie1(matDistSp,matDistT,p=estima$par)

Kriging Simple

library(rgdal)
# Coordenadas de latitud y longitud
lat <- 5.3332729249446365
lon <- -75.4857825511143
# Crear objeto SpatialPoints
coords <- data.frame(lon = lon, lat = lat)
coords_sp <- SpatialPoints(coords, proj4string =
                             CRS("+proj=longlat +datum=WGS84"))
# Transformar a UTM
coords_utm <- spTransform(coords_sp, CRS("+proj=utm +zone=18 +datum=WGS84"))
# Obtener las coordenadas UTM
grillaSpT0=rbind(cbind(x1,x2,t),c(446176.4,589525.2,19))
matDistSp0=as.matrix(dist(grillaSpT0[,1:2]))
matDistT0=as.matrix(dist(grillaSpT0[,3:3]))
sigma0=cressie1(matDistSp0,matDistT0,p=estima$par)
#vector de covarianzas entre la coordenada a predecir y las observadas
lambda1=solve(sigma1)%*%sigma0[162,-162]
mu1<-predict(mod10, data.frame(Este2=446176.4,Norte2=589525.2,Hora2=19))
z_pred0=t(lambda1)%*%datos1[,4]+mu1
z_pred0

##          [,1]
## [1,] 1.434491

Error cuadrático medio

ECM1<-t(sigma0[162,-162])%*%sigma1%*%sigma0[162,-162]
ECM1

##           [,1]
## [1,] 0.4690119

Kriging espacio tiempo con el modelo C R E S S I E - H U A N G (1999)

Este método se basa en las representaciones espectrales de las funciones de covarianza; se construyen familias paramétricas de covarianzas estacionarias, encontrando transformadas inversas de Fourier de densidades espectrales. Este procedimiento no impone restricciones de separabilidad, aunque podrían surgir modelos separables como casos particulares. (Bohorquez M, 2010)

Definición de la función CRESSIE-HUANG

library(SimDesign)
CH_1=function(h,u,p){(p[1]^2/((p[2]^2*u^2+1)))*
    exp(-(p[3]^2*h^2)/(p[2]^2*u^2+1))}

Planteamiento de la Log-Verosimilitud

LV<-function(p,h,u,modelo,z){0.5*(length(z)*log(2*pi)+
    log(det(modelo(p,h=matDistSp,u=matDistT)))+
      t(z)%*%solve(modelo(p,h=matDistSp,u=matDistT))%*%z)}

Encontrando los parámetros óptimos

p0 <- c(1,1,1) 
lo <- c(0.001,0.001,0.01)   
estima2=optim(p0,LV,h=matDistSp,u=matDistT,modelo=CH_1,z=sim1,
              hessian = T,method = "BFGS") 
# se le puede dar la opcion de la matriz hessiana

Cálculo del semivariograma con los parámetros óptimos

sigma2=CH_1(matDistSp,matDistT,p=c(estima2$par))

kriging simple

sigma02=CH_1(matDistSp0,matDistT0,p=c(estima2$par))
#vector de covarianzas entre la coordenada a predecir y las observadas
lambda2=solve(sigma2)%*%sigma02[163,-163]
mu2<-predict(mod10, data.frame(Este2=446176.4,Norte2=589525.2,Hora2=19))
z_pred02=t(lambda2)%*%datos1[,4]+mu2
z_pred02

##          [,1]
## [1,] 1.363263

Error cuadrático medio

ECM2<-t(sigma02[162,-162]) %*% sigma2 %*% sigma02[162,-162]
ECM2

##             [,1]
## [1,] 0.002828884

Comparación entre las varianzas del error de predicción entre los modelos propuestos

Modelo	ECM obtenido
CRESSIE-HUANG	0.0028
Cressie1	0.469

Por lo que se puede ver que el modelo CRESSIE-HUANG es mucho más adecuado en este caso para el modelamiento del semivariograma espacio-tiempo donde la variable de interés es la velocidad del viento.Por lo que para hacer la predicción en la grilla completa tomamos como modelo de semivariograma teórico el modelo Cressie-Huang con parámetros: \(d=-0.3268804\), \(b=0.5311734\), \(a=1.0000000\). Donde el parámetro de escala en el tiempo es \(a\), el parámetro de escala en el espacio es \(b\) y el parámetro \(d\) es la dimensión espacial.

Kriging en una grilla de 2000 puntos no muestreados juntos con sus respectivos rezagos temporales:

library(rgdal)
par(mfrow=c(1,1))
municipios_sf <- st_read("F:/MI PC/Archivos/UNAL/2023-1/Espacial/Parcial1/parcial1/shp/MGN_MPIO_POLITICO.shp",quiet=TRUE)

Obtener mapa del departamento de Caldasy crear la grilla con los 2000 puntos no muestreados del departamentode Caldas

Cal <- municipios_sf[municipios_sf$DPTO_CNMBR %in% c("CALDAS"),]

st_write(Cal, "Cal.shp",append=FALSE)

## Deleting layer `Cal' using driver `ESRI Shapefile'
## Writing layer `Cal' to data source `Cal.shp' using driver `ESRI Shapefile'
## Writing 27 features with 12 fields and geometry type Multi Polygon.

map <- readOGR("Cal.shp")

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\Lenovo\Desktop\Cal.shp", layer: "Cal"
## with 27 features
## It has 12 fields

class(map)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

new <- sp::spsample(map, n = 2000, type = "regular")

Conversión a coordenadas UTM

CRS_UTM_NY = (c("+proj=utm +zone=18 +datum=NAD83 
                  +units=m +no_defs"))
utm_coord2 = spTransform(new, CRS(CRS_UTM_NY))
plot(utm_coord2)

utm_coord2<-as.data.frame(utm_coord2)
m<-nrow(utm_coord2)
utm_coord2<-utm_coord2[rep(1:m, each = 6), ]
time<-18:23
utm_coord2$t<-rep(time,m)
utm_coord2<-as.data.frame(cbind(utm_coord2$x1,utm_coord2$x2,utm_coord2$t))
colnames(utm_coord2)<-c('x1','x2','t')

Cálculo de la matriz de distancias

# Obtener las coordenadas UTM
x1 <- Caldas_utm$Este
x2 <- Caldas_utm$Norte
t <- Caldas_utm$Hora
grillaSpT=cbind(x1,x2,t)
grillaSpT0=rbind(cbind(x1,x2,t),utm_coord2)
matDistSp0=as.matrix(dist(grillaSpT0[,1:2]))
matDistT0=as.matrix(dist(grillaSpT0[,3:3]))
sigma03=CH_1(matDistSp0,matDistT0,p=c(estima2$par))

Realización del Kriging espacio tiempo con el modelo Cressie Huang

V<-NULL
PRECT<-NULL
for(i in 1:12000){
  grillaS=rbind(cbind(x1,x2,t),utm_coord2[i,])
  matSp0=as.matrix(dist(grillaS[,1:2]))
  mattT0=as.matrix(dist(grillaS[,3:3]))
  sigma03=CH_1(matSp0,mattT0,p=c(estima2$par))
  lambda=solve(sigma2)%*%sigma03[163,-163]
  mu<-predict(mod10, data.frame(Este2=utm_coord2[i,1],Norte2=utm_coord2[i,2],Hora2=utm_coord2[i,3]))
  z_pred=t(lambda)%*%datos1[,4]+mu
  PRECT[i]<-z_pred
  sz<-var(Caldas$Velocidad_Viento)
  G<-t(lambda)%*%sigma03[163,-163]
  v[i]<-t(sigma02[162,-162]) %*% sigma2 %*% sigma02[162,-162]
}

Mapa de predicción visualizado por tiempo

utm_coord2<-utm_coord2[1:12000,]
PRECT<-as.matrix(PRECT[1:12000])
PRECT18<-cbind((utm_coord2[utm_coord2$t %in% c(18),])[,1:2],PRECT[utm_coord2$t %in% c(18),])
colnames(PRECT18)<-c('Este','Norte','vvpred')
VARP18<-cbind((utm_coord2[utm_coord2$t %in% c(18),])[,1:2],
              tprec.ok$var1.var[1:2000])
colnames(VARP18)<-c('Este','Norte','varpred')
PRECT19<-cbind(utm_coord2[utm_coord2$t %in% c(19),][,1:2],PRECT[utm_coord2$t %in% c(19),])
colnames(PRECT19)<-c('Este','Norte','vvpred')
VARP19<-cbind((utm_coord2[utm_coord2$t %in% c(19),])[,1:2],
              tprec.ok$var1.var[2001:4000])
colnames(VARP19)<-c('Este','Norte','varpred')
PRECT20<-cbind(utm_coord2[utm_coord2$t %in% c(20),][,1:2],PRECT[utm_coord2$t %in% c(20),])
colnames(PRECT20)<-c('Este','Norte','vvpred')
VARP20<-cbind((utm_coord2[utm_coord2$t %in% c(20),])[,1:2],
              tprec.ok$var1.var[4003:6002])
colnames(VARP20)<-c('Este','Norte','varpred')
PRECT21<-cbind(utm_coord2[utm_coord2$t %in% c(21),][,1:2],PRECT[utm_coord2$t %in% c(21),])
colnames(PRECT21)<-c('Este','Norte','vvpred')
VARP21<-cbind((utm_coord2[utm_coord2$t %in% c(21),])[,1:2],
              tprec.ok$var1.var[6006:8005])
colnames(VARP21)<-c('Este','Norte','varpred')
PRECT22<-cbind(utm_coord2[utm_coord2$t %in% c(22),][,1:2],PRECT[utm_coord2$t %in% c(22),])
colnames(PRECT22)<-c('Este','Norte','vvpred')
VARP22<-cbind((utm_coord2[utm_coord2$t %in% c(22),])[,1:2],
              tprec.ok$var1.var[8008:10007])
colnames(VARP22)<-c('Este','Norte','varpred')
PRECT23<-cbind(utm_coord2[utm_coord2$t %in% c(23),][,1:2],PRECT[utm_coord2$t %in% c(23),])
colnames(PRECT23)<-c('Este','Norte','vvpred')
VARP23<-cbind((utm_coord2[utm_coord2$t %in% c(23),])[,1:2],
              tprec.ok$var1.var[10010:12009])
colnames(VARP23)<-c('Este','Norte','varpred')

library(ggplot2)
library(patchwork)

## Warning: package 'patchwork' was built under R version 4.2.3

coordenadas18 <- PRECT18[, c("Este", "Norte")]
ubicaciones18 <- SpatialPointsDataFrame(coords = coordenadas18, data = PRECT18,
                                        proj4string = CRS("+proj=longlat +datum=WGS84"))
ubicaciones18<-as.data.frame(ubicaciones18)
mapa1 <- ggplot() +
  geom_point(data = ubicaciones18, aes(x = Este, y = Norte, color = vvpred),size=3) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "Map.Pred.t=18") +
  xlab("Este") +
  ylab("Norte") +
  theme_minimal()

coordenadas19 <- PRECT19[, c("Este", "Norte")]
ubicaciones19 <- SpatialPointsDataFrame(coords = coordenadas19, data = PRECT19,
                                        proj4string = CRS("+proj=longlat +datum=WGS84"))
ubicaciones19<-as.data.frame(ubicaciones19)
mapa2 <- ggplot() +
  geom_point(data = ubicaciones19, aes(x = Este, y = Norte, color = vvpred),size=3) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "Map.Pred.t=19") +
  xlab("Este") +
  ylab("Norte") +
  theme_minimal()

coordenadas20 <- PRECT20[, c("Este", "Norte")]
ubicaciones20 <- SpatialPointsDataFrame(coords = coordenadas20, data = PRECT20,
                                        proj4string = CRS("+proj=longlat +datum=WGS84"))
ubicaciones20<-as.data.frame(ubicaciones20)
mapa3 <- ggplot() +
  geom_point(data = ubicaciones20, aes(x = Este, y = Norte, color = vvpred),size=3) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "Map.Pred.t=20") +
  xlab("Este") +
  ylab("Norte") +
  theme_minimal()

coordenadas21 <- PRECT21[, c("Este", "Norte")]
ubicaciones21 <- SpatialPointsDataFrame(coords = coordenadas21, data = PRECT21,
                                      proj4string = CRS("+proj=longlat +datum=WGS84"))
ubicaciones21<-as.data.frame(ubicaciones21)
mapa4 <- ggplot() +
  geom_point(data = ubicaciones21, aes(x = Este, y = Norte, color = vvpred),size=3) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "Map.Pred.t=21") +
  xlab("Este") +
  ylab("Norte") +
  theme_minimal()

coordenadas22 <- PRECT22[, c("Este", "Norte")]
ubicaciones22 <- SpatialPointsDataFrame(coords = coordenadas22, data = PRECT22,
                                        proj4string = CRS("+proj=longlat +datum=WGS84"))
ubicaciones22<-as.data.frame(ubicaciones22)
mapa5 <- ggplot() +
  geom_point(data = ubicaciones22, aes(x = Este, y = Norte, color = vvpred),size=3) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "Map.Pred.t=22") +
  xlab("Este") +
  ylab("Norte") +
  theme_minimal()

coordenadas23 <- PRECT23[, c("Este", "Norte")]
ubicaciones23 <- SpatialPointsDataFrame(coords = coordenadas23, data = PRECT23,
                                        proj4string = CRS("+proj=longlat +datum=WGS84"))
ubicaciones23<-as.data.frame(ubicaciones23)
mapa6 <- ggplot() +
  geom_point(data = ubicaciones23, aes(x = Este, y = Norte, color = vvpred),size=3) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "Map.Pred.t=23") +
  xlab("Este") +
  ylab("Norte") +
  theme_minimal()

# Combinamos los gráficos
mapa1 + mapa2 + mapa3 + mapa4 +mapa5+mapa6+
  plot_layout(ncol = 2)

Mapa de varianza de predicción

library(ggplot2)
library(patchwork)

coordenadas18 <- VARP18[, c("Este", "Norte")]
ubicaciones18 <- SpatialPointsDataFrame(coords = coordenadas18, data = VARP18,
                                        proj4string = CRS("+proj=longlat +datum=WGS84"))
ubicaciones18<-as.data.frame(ubicaciones18)
mapa1 <- ggplot() +
  geom_point(data = ubicaciones18, aes(x = Este, y = Norte, color = varpred),size=3) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "Map.Var.Pred.t=18") +
  xlab("Este") +
  ylab("Norte") +
  theme_minimal()

coordenadas19 <- VARP19[, c("Este", "Norte")]
ubicaciones19 <- SpatialPointsDataFrame(coords = coordenadas19, data = VARP19,
                                        proj4string = CRS("+proj=longlat +datum=WGS84"))
ubicaciones19<-as.data.frame(ubicaciones19)
mapa2 <- ggplot() +
  geom_point(data = ubicaciones19, aes(x = Este, y = Norte, color = varpred),size=3) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "Map.var.Pred.t=19") +
  xlab("Este") +
  ylab("Norte") +
  theme_minimal()

coordenadas20 <- VARP20[, c("Este", "Norte")]
ubicaciones20 <- SpatialPointsDataFrame(coords = coordenadas20, data = VARP20,
                                        proj4string = CRS("+proj=longlat +datum=WGS84"))
ubicaciones20<-as.data.frame(ubicaciones20)
mapa3 <- ggplot() +
  geom_point(data = ubicaciones20, aes(x = Este, y = Norte, color = varpred),size=3) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "Map.var.Pred.t=20") +
  xlab("Este") +
  ylab("Norte") +
  theme_minimal()

coordenadas21 <- VARP21[, c("Este", "Norte")]
ubicaciones21 <- SpatialPointsDataFrame(coords = coordenadas21, data = VARP21,
                                      proj4string = CRS("+proj=longlat +datum=WGS84"))
ubicaciones21<-as.data.frame(ubicaciones21)
mapa4 <- ggplot() +
  geom_point(data = ubicaciones21, aes(x = Este, y = Norte, color = varpred),size=3) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "Map.var.Pred.t=21") +
  xlab("Este") +
  ylab("Norte") +
  theme_minimal()

coordenadas22 <- VARP22[, c("Este", "Norte")]
ubicaciones22 <- SpatialPointsDataFrame(coords = coordenadas22, data = VARP22,
                                        proj4string = CRS("+proj=longlat +datum=WGS84"))
ubicaciones22<-as.data.frame(ubicaciones22)
mapa5 <- ggplot() +
  geom_point(data = ubicaciones22, aes(x = Este, y = Norte, color = varpred),size=3) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "Map.var.Pred.t=22") +
  xlab("Este") +
  ylab("Norte") +
  theme_minimal()

coordenadas23 <- PRECT23[, c("Este", "Norte")]
ubicaciones23 <- SpatialPointsDataFrame(coords = coordenadas23, data = VARP23,
                                        proj4string = CRS("+proj=longlat +datum=WGS84"))
ubicaciones23<-as.data.frame(ubicaciones23)
mapa6 <- ggplot() +
  geom_point(data = ubicaciones23, aes(x = Este, y = Norte, color = varpred),size=3) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "Map.var.Pred.t=23") +
  xlab("Este") +
  ylab("Norte") +
  theme_minimal()

# Combinamos los gráficos
mapa1 + mapa2 + mapa3 + mapa4 +mapa5+mapa6+
  plot_layout(ncol = 2)

Conclusiones

Los resultados presentados dan a conocer que la predicción realizada estuvo dentro de los límites naturales en los que varía la velocidad del viento en el departamento de Caldas-región Andina. Somos concientes de la necesidad de un futuro proyecto incluir como covariable a la altura ya que esta es un factor influyente en la variable de estudio. El Kriging espacio tiempo como se desarrolló en este trabajo no es muy recomendable dada las limitaciones computacionales y la gran cantidad de datos, así como la falta de claridad en la doumentación disponible de los paquetes de R disponibles para realizar este tipo de Kriging, ya que esto dificulta y demora en gran medida la obtención de dichas predicciones en una grilla con la mayor cantidad de puntos posibles. En nuestro caso decidimos emplear las funciones “manuales” para este Kriging dadas por la profesora Martha dado que los paquetes de R nos arrojaban errores en las matrices. En la predicción así como el comportamiento particular que se vislumbró desde el inicio con respecto a la velocidad del viento a las 6pm se evidenció una escala de variación de la velocidad del viento mucho menor, pues en las horas de las 19 a las 23 la velocidad del viento alcanza un valor d 1.6 m/s pero a las 18 horas alcanza escasamente un poco más de 1 m/s, lo que nos lleva a plantear la pregunta **¿en el departameto de Caldas la velocidad del viento aumenta a medida que el día va declinando?.

Bibliografía

Bohorquez, Martha (2010). Tesis de Maestría:DIFERENCIABILIDAD DE FUNCIONES DE COVARIANZA ESPACIO TEMPORAL NO SEPARABLES. Recuperado de: https://repositorio.unal.edu.co/bitstream/handle/unal/6750/196875.2010.pdf?sequence=1

Bohorquez,Martha(2023). Notas de Clases: Estadística Espacial. Universidad Nacional de Colombia.

Kriging Espacio Tiempo de la velocidad del viento en el departamento de Caldas-Colombia

Jessica López Mejía y Nicolás Felipe Romero Solano