Estudio geoespacial de la velocidad del viento en la región Andina Colombiana del día 12 de febrero del 2023 a las 6:00pm
Jessica López Mejía y Nicolás Felipe Romero Solano
2023-03-14
Introducción
La estadística espacial es utilizada para analizar patrones de variabilidad climática a lo largo del territorio colombiano. Esto es importante especialmente en nuestro país debido a su gran diversidad geográfica, lo que se traduce en una gran variabilidad climática. Ya que por ejemplo, mientras que en algunas regiones de Colombia hay selvas tropicales y lluvias constantes, en otras regiones hay desiertos y zonas semiáridas.
Por consiguiente, la estadística espacial puede ser utilizada para analizar datos climáticos, como la temperatura, la precipitación, la humedad, la velocidad del viento, entre otros, en diferentes ubicaciones geográficas. Esto permite identificar patrones espaciales de variabilidad climática, como zonas con altas temperaturas, zonas con altas precipitaciones o zonas con fuertes vientos. Estos patrones pueden ser utilizados para entender mejor el clima en Colombia con el propósito de desarrollar estrategias de adaptación y mitigación frente al cambio climático.
Pero, ¿cómo influye la velocidad del viento (variable de interés) en el clima colombiano?
La velocidad del viento afecta la distribución de: la temperatura, pues el viento puede llevar aire frío o caliente a diferentes regiones de Colombia, lo que puede afectar la temperatura local. Por ejemplo, los vientos alisios que soplan desde el este sobre el Océano Pacífico pueden traer aire frío y húmedo a la región pacífica, lo que contribuye a la baja temperatura promedio de esta región; las precipitaciones, ya que la velocidad del viento puede afectar la formación de nubes y la distribución de las precipitaciones. Por ejemplo, los vientos alisios pueden favorecer la formación de nubes y la lluvia en la región pacífica, mientras que en la región Caribe, los vientos pueden impedir la formación de nubes y disminuir la cantidad de precipitaciones; los fenómenos meteorológicos, dado que los vientos fuertes pueden contribuir a la formación y propagación de estos eventos.
En cuanto a la velocidad del viento en la región andina que es la zona donde fueron recolectados los datos se puede decir que en las zonas más altas de los Andes, como en los páramos y las cumbres de las montañas, la velocidad del viento suele ser alta debido a las condiciones topográficas y a la mayor altura. Por ejemplo, en la cima del Nevado del Ruiz, que se encuentra a una altitud de más de 5.000 metros sobre el nivel del mar, la velocidad promedio del viento puede superar los 100 km/h en algunos momentos del año.
En las zonas intermedias de la región andina, donde se encuentran la mayoría de las ciudades y poblaciones, la velocidad del viento suele ser moderada. Sin embargo, también pueden presentarse vientos fuertes en algunas épocas del año, especialmente durante la temporada de lluvias.
Es importante destacar que la velocidad del viento en la región andina puede influir en diferentes aspectos, como la dispersión de contaminantes atmosféricos, la generación de energía eólica, la agricultura y la seguridad de las infraestructuras. Por lo tanto, es importante realizar mediciones y estudios específicos sobre la velocidad del viento en diferentes zonas de la región para entender mejor su comportamiento y sus posibles impactos.
características generales de la Región Andina Colombiana
La región andina colombiana está compuesta por los siguientes departamentos: Antioquia, Boyacá, Caldas, Caquetá, Cauca, Cundinamarca, Huila, Nariño, Norte de Santander, Quindío, Risaralda, Santander, Tolima y Valle del Cauca.
Se caracteriza por tener altitudes elevadas, que van desde los 1.000 metros sobre el nivel del mar hasta más de 5.000 metros en las zonas más altas de los Andes.
Su clima es variado, debido a las diferentes altitudes y ubicaciones geográficas. En general, se pueden identificar dos estaciones: una estación seca y una estación lluviosa.
La vegetación de la región andina es muy diversa, debido a las diferentes altitudes y condiciones climáticas. En las zonas más altas, como en los páramos, se encuentran ecosistemas de alta montaña, con una vegetación adaptada a las bajas temperaturas y a la humedad. En las zonas intermedias, se pueden encontrar bosques y sabanas, mientras que en las zonas más bajas, como en el valle del río Magdalena, se pueden encontrar bosques secos y selvas tropicales.
Los datos obtenidos para el presente trabajo fueron los siguientes (tomado de datos.gov.co)
A continuación se muestran los 15 primeros registros de la tabla de datos cuyos campos corresponde a: Latitud, Longitud, Departamento, Municipio, Velocidad_Viento.
| Latitud | Longitud | Departamento | Municipio | Velocidad_Viento |
|---|---|---|---|---|
| 5.79 | -75.43 | ANTIOQUIA | ABEJORRAL | 1.20 |
| 8.08 | -73.22 | NORTE DE SANTANDER | ABREGO | 2.50 |
| 6.63 | -76.06 | ANTIOQUIA | ABRIAQUI | 1.00 |
| 1.81 | -75.89 | HUILA | ACEVEDO | 1.30 |
| 2.26 | -75.77 | HUILA | AGRADO | 3.10 |
| 4.38 | -74.67 | CUNDINAMARCA | AGUA DE DIOS | 1.00 |
| 8.31 | -73.61 | CESAR | AGUACHICA | 3.00 |
| 6.16 | -73.52 | SANTANDER | AGUADA | 0.80 |
| 5.61 | -75.45 | CALDAS | AGUADAS | 1.00 |
| 3.22 | -75.24 | HUILA | AIPE | 1.10 |
| 4.88 | -74.44 | CUNDINAMARCA | ALBAN | 0.70 |
| 1.47 | -77.08 | NARIÑO | ALBAN | 1.70 |
| 5.76 | -73.91 | SANTANDER | ALBANIA | 0.80 |
| 4.67 | -75.78 | VALLE DEL CAUCA | ALCALA | 0.60 |
| 0.88 | -77.70 | NARIÑO | ALDANA | 2.50 |
Exploración y descripción de los datos
Como se pudo evidenciar en la anterior tabla de datos, estos corresponden a datos climáticos de Colombia en diferentes departamentos alrededor de la región Andina. Como el objetivo de este trabajo es hacer un estudio geoespacial (no temporal) de una sola variable de interés, se seleccionó la variable correspondiente a velocidad del viento,además se filtraron los datos (cuyos campos corresponden a: altitud, longitud y velocidad del aire), teniendo en cuenta una hora (06:00 pm) y un día específico(2023-02-12). Conforme a ello se obtuvo un total de 780 observaciones de la velocidad del viento en diferentes ubicaciones geoespaciales de dicha región.
nrow(Randina)## [1] 780Veamos un resumen de la variable velocidad del viento:
## [,1]
## Media "1.34858974358974"
## Mediana "1.2"
## Varianza "0.344991590138573"
## Desviación_Estandar "0.587359847230447"
## Cv "43.5536344557228"
## MIN "0.3"
## DEPVMIN "SANTANDER"
## MAX "3.3"
## DEPVMAX_1 "ANTIOQUIA"
## DEPVMAX_2 "CUNDINAMARCA"Se observa que de acuerdo al resumen mostrado de la variable velocidad del viento la máxima velocidad registrada se ubicó en los departamentos Cundinamarca y Antioquia y la mínima velocidad del viento se ubicó en Santander.
Se puede evidenciar por medio del diagrama de caja que hay pocos valores atípicos en la variable respuesta.Visualicemos los datos en un mapa para identificar patrones o áreas donde los valores sean particularmente altos o bajos en comparación con el resto del conjunto de datos (como se identificó en la ubicación espacial del valor máximo de velocidad del viento registrada y del mínimo valor del mismo).
## coords.x1 coords.x2 Randina.Velocidad_Viento
## 1 -75.42874 5.789301 1.2
## 2 -73.22173 8.081782 2.5
## 3 -76.06430 6.632282 1.0
## 4 -75.88854 1.805264 1.3
## 5 -75.77190 2.259730 3.1
## 6 -74.66922 4.375315 1.0
Convertir datos a coordenadas
En general se tiene que en el mes de febrero en la región Andina la velocidad del viento puede variar según la ubicación geográfica específica dentro de la región. Sin embargo, en general, febrero es un mes con condiciones de viento moderado a fuerte especialmente en áreas montañosas y altiplanos.
En Bogotá, por ejemplo, la velocidad promedio del viento en febrero es de alrededor de 6,2 km/h, con ráfagas ocasionales que pueden superar los 20 km/h.
Veamos la correlación entre la velocidad del viento con cada una de las coordenadas
## Este Norte Velocidad_Viento
## Este 1.00000000 0.7145658 -0.05683464
## Norte 0.71456581 1.0000000 -0.17121001
## Velocidad_Viento -0.05683464 -0.1712100 1.00000000Se evidencia una muy baja correlación lineal entre la variable velocidad del viento con respecto a cada una de las coordenadas.
Este gráfico nos proporciona más claridad en cuanto al tipo de relación posible en los datos con respecto a las coordenadas cartesianas,pues muestra cierta curvatura, por lo que podría sugerir un modelo polinomial de segundo grado para remover la tendencia.
Otros gráficos descriptivos
Dado que las isolíneas en un mapa de contornos representan los valores de la variable de interés (en este caso Velocidad del viento) en un espacio geográfico, se puede decir que el gráfico muestra que hay una gran variación en la velocidad del viento donde las líneas se ven muy juntas, por lo que en estas zonas de la región Andina la velocidad del viento cambia rápidamente en el espacio por el contrario en las zonas donde se ven las isolíneas bastante separadas indica que la Velocidad del viento en estas zonas está cambiando lentamente en el espacio.
Revisando tendecia en los datos con respecto a cada coordenada
Se evidencia una tendencia no lineal en las observaciones registradas con respecto a cada coordenada.
Corrigiendo la tendencia
Para ello estimamos varios modelos lineales y polinómicos para ver cuál de ellos teniendo en cuenta lo mostrado anteriormente y el principio de parsimonia nos ayuda a remover la tendencia.
Posibles modelos
modelo1<-lm(Z~1)
modelo2<-lm(Z~y+x)
modelo3<-lm(Z~y+x+I(y^2)+I(x^2))
modelo4<-lm(Z~y+x+I(y^2)+I(x^2)+I(x*y)+I(x^2*y^2)+I(x^2*y)+I(y^2*x))Donde en cada modelo obtenemos los siguientes R cuadrados. Para el modelo2, 0.03808; modelo3, 0.1506; modelo4, 0.213.
¿Cuál de los anteriores modelos propuestos nos ayudará a corregir la tendencia? Veamos:
## Number of data points: 780
##
## Coordinates summary
## Este Norte
## min 150838.4 89215.6
## max 819564.9 960688.6
##
## Distance summary
## min max
## 1719.442 1022659.232
##
## Data summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.30000 0.90000 1.20000 1.34859 1.70000 3.30000
Removiendo tendencia con modelo 1
Removiendo tendencia con modelo 2
Removiendo tendencia con modelo 3
Removiendo tendencia con modelo 4
Conforme a lo mostrado y lo indicado al inicio en la descripción de los datos, se elige el modelo polinimial de segundo grado (en R trend=“2nd”) para remover la tendencia.
Estimación del semivariograma empírico removiendo la tendencia
Debido a que tenemos algunos problemas de datos atípicos, el metodo de estimación del semivariograma empírico que eligiremos será el estimador resistente a valores atípicos propuesto por Cressie (“modulus”).
## variog: computing omnidirectional variogram
## variog: computing omnidirectional variogram
## variog: computing omnidirectional variogram
Observamos que la distribución de los residuales mejora bastante cuando se remueve la tendencia con el polinomio de grado 2, por lo que nos quedamos con el semivariograma empírico al cual se le removió esta tendencia.
Explorando Anisotropía
## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
En estadística espacial, la anisotropía se refiere a la propiedad de que las características de los datos varían de manera diferente en diferentes direcciones. Es decir, la estructura espacial de los datos no es uniforme en todas las direcciones, sino que puede ser más fuerte o más débil en diferentes direcciones.
Por ejemplo, si se está analizando la distribución de la temperatura en una región geográfica, es posible que la temperatura sea más alta en una dirección que en otra debido a factores como la topografía, la latitud o la exposición al sol. En este caso, se dice que la temperatura es anisotrópica.
Por lo que de acuerdo a los resultados obtenidos mostrados en los diferentes plots, se puede decir que la velocidad del viento en la región Andina Colombiana para el 12 de febrero del 2023 a las 6:00 pm es Anisotrópica pues la estructura de los semivariogramas varía conforme cambia la dirección.
Estimación del semivariograma teórico
Partiendo del semivariograma empírico se realiza una aproximación al semivariograma teórico desde el cual se obtienen los valores iniciales de los parámetros del modelo que mejor se ajuste por medio de eyefit, el cual correspondió al modelo powered.exponential (disponible en R) con parámetros iniciales: silla=0.31, pepita=0.04, rango=327775.39, kappa=0.5.
Para llegar a este modelo, hicimos primeramente varios intentos con otros modelos y realizamos la comparación calculando el valor Q (mínimos cuadrados ordinarios), el cual corresponde a: \[Q=\sum_{i=1}^{n}(\hat{\gamma}(h)-\gamma(h))^2 \] Donde \(\gamma(h)\) corresponde al modelo teórico del semivariograma con los valores iniciales dados y \(\hat{\gamma}(h)\) corresponde al semivariograma empírico ya calculado antes, \(n\) en este caso es \(13\) y \(h\) es el vector de distancias. Se obtuvo los siguientes resultados:
| cov.model | sigmasq | phi | tausq | kappa | kappa2 | practicalRange |
|---|---|---|---|---|---|---|
| powered.exponential | 0.31 | 327775.39 | 0.04 | 0.5 | 2941591.34573194 | |
| Matern | 0.16 | 106305.53 | 0.16 | 2.03 | 574389.064109716 | |
| spherical | 0.13 | 664409.58 | 0.16 | 664409.58 | ||
| Wave | 0.1 | 106305.53 | 0.19 | 318008.361624559 |
## variog: computing omnidirectional variogram## variog: computing omnidirectional variogram## Q
## powered.exponential 0.1078955
## matern 0.3913899
## spherical 0.4950147
## wave 0.5848878Por lo que finalmente se elige el modelo powered.exponential, debido a que tiene el menor valor Q, con los valores iniciales ya indicados.Este modelo se define de la siguiente manera:\(\gamma(h) =0\) si \(||h||=0\), para \(||h||\neq 0\) queda definido como se muestra:
\[\gamma(h) = c \left[1 - \exp\left(-\left(\frac{h}{a}\right)^p\right)\right]\]
Donde:
\(\gamma(h)\) es la semivarianza para un par de puntos separados por una distancia \(h\). \(c\) es la varianza de la variable regionalizada. \(a\) es el rango espacial, que representa la distancia más allá de la cual la correlación espacial es insignificante. \(p\) es el parámetro de suavizado, que controla la tasa a la que disminuye la semivarianza a medida que aumenta la distancia \(h\).
Se muestra a continuación los diferentes métodos aplicados para realizar el modelo teórico.
## cov.model sigmasq phi tausq kappa kappa2 practicalRange
## 1 powered.exponential 0.31 327775.39 0.04 0.5 <NA> 2941591.34573194Utilizando Mínimos cuadrados ordinarios
fitvar1 <- variofit(v3,
cov.model = "powered.exponential",
ini1,
fix.nugget = TRUE,
nugget =0.04,
wei = "equal") ## variofit: covariance model used is powered.exponential
## variofit: weights used: equal
## variofit: minimisation function used: optimUtilizando la matriz de ponderación “npairs”
fitvar2 <- variofit(v3,
cov.model = "powered.exponential",
ini1,
fix.nugget = TRUE,
nugget = 0.04,
wei = "npairs")## variofit: covariance model used is powered.exponential
## variofit: weights used: npairs
## variofit: minimisation function used: optimUtilizando la matriz de ponderación propuesta por “Cressie(1985)”
fitvar3 <- variofit(v3,
ini1,
fix.nugget = TRUE,
nugget = 0.04,
wei = "cressie")## variofit: covariance model used is powered.exponential
## variofit: weights used: cressie
## variofit: minimisation function used: optimUtilizando el método basado en máxima verosimilitud
fitvar4 <- likfit(DATAg,
coords = DATAg$coords,
data = DATAg$data,
trend = "2nd",
ini.cov.pars = ini1,
fix.nugget = T,
nugget = 0.04,
cov.model = "powered.exponential",
lik.method = "ML") ## ---------------------------------------------------------------
## 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.Gráfico que permite comparar los cuatro semivariogramas teóricos obtenidos
Finalmente, vemos que los modelos teóricos para el semivariograma, obtenidos bajo los cuatro métodos expuestos son muy similares.
## variofit: model parameters estimated by OLS (ordinary least squares):
## covariance model is: powered.exponential with fixed kappa = 0.5
## fixed value for tausq = 0.04
## parameter estimates:
## sigmasq phi
## 0.3726 327775.3900
## Practical Range with cor=0.05 for asymptotic range: 2941591
##
## variofit: minimised sum of squares = 0.0182## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: powered.exponential with fixed kappa = 0.5
## fixed value for tausq = 0.04
## parameter estimates:
## sigmasq phi
## 0.2942 125706.9344
## Practical Range with cor=0.05 for asymptotic range: 1128146
##
## variofit: minimised weighted sum of squares = 78.5135## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: powered.exponential with fixed kappa = 0.5
## fixed value for tausq = 0.04
## parameter estimates:
## sigmasq phi
## 0.3856 327775.3900
## Practical Range with cor=0.05 for asymptotic range: 2941591
##
## variofit: minimised weighted sum of squares = 2440.981## likfit: estimated model parameters:
## beta0 beta1 beta2 beta3 beta4 beta5
## "2.157e+00" "0.000e+00" "0.000e+00" "0.000e+00" "0.000e+00" "0.000e+00"
## sigmasq phi
## "3.434e-01" "3.278e+05"
## Practical Range with cor=0.05 for asymptotic range: 2941591
##
## likfit: maximised log-likelihood = -212.6Vemos que el “minimised sum of squares” mejor es de fitvar1 y sus valores estimados son
# sigmasq phi
# 0.3726 327775.3900Los parámetros óptimos mediante ML son:
# sigmasq phi
# "3.434e-01" "3.278e+05" Los parámetros óptimos mediante REML son:
# sigmasq phi
# "3.434e-01" "3.278e+05" Kriging Simple para predecir la velocidad del viento en un punto no muestreado
Deseamos predecir la velocidad del viento en un lugar no muestreado de la región Andina Colombiana como Palmira Valle del Cauca cuya latitud y longitud son 3.53944, -76.30361 y en coordenadas utm corresponden a: Este: 355203,20 Norte: 391321,20. Tomamos el modelo teórico estimado por el método de MCO modelo exponencial potencial con silla 0.3726 y rango 327775.3900.
\[\gamma(h) = 0.3726 \left[1 - \exp\left(-\left(\frac{||h||}{327775.39}\right)^{0.5}\right)\right]\] Para ello aplicamos Kriging simple dado que podemos dar por conocida la media pues la variable de interés se encuentra centrada (con media 0) y la tendencia de los datos fue inicialmente removida con un modelo polinomial de segundo grado.
# modelo exponencial potencial
# sigmasq phi
# 0.3726 327775.3900
h<-as.matrix(dist(DATAg$coords[,1:2]))
sigma1<-as.matrix(cov.spatial(h,"powered.exponential",
cov.pars = c(0.3726,327775.39), kappa = 0.5))
s0<-c(355203.2,391321.2)
h2<-NULL
for(i in 1:780){
h2[i]<-dist(rbind(s0,DATAg$coords[i,1:2]))
}
k<-cov.spatial(h2,"powered.exponential",
cov.pars = c(0.3726,327775.39), kappa = 0.5)
#Hallamos los pesos (lambda)
lambda<-solve(sigma1)%*%k
head(lambda)## [,1]
## 1 -1.083154e-07
## 2 3.890253e-08
## 3 -4.082297e-06
## 4 -7.931770e-06
## 5 -4.566251e-06
## 6 -6.071453e-07sum(lambda)## [1] 1.000067#Predecimos la velocidad del viento en la ubicación Este: 355203.20
#Norte: 391321.20
a<-s0[1]
b<-s0[2]
#Parametros del modelamiento de la media
modelo4##
## Call:
## lm(formula = Z ~ y + x + I(y^2) + I(x^2) + I(x * y) + I(x^2 *
## y^2) + I(x^2 * y) + I(y^2 * x))
##
## Coefficients:
## (Intercept) y x I(y^2) I(x^2)
## 4.583e+00 8.193e-06 -2.417e-05 -1.149e-11 5.035e-11
## I(x * y) I(x^2 * y^2) I(x^2 * y) I(y^2 * x)
## -1.836e-11 -4.871e-24 -5.131e-17 4.820e-17mu<- 4.583e+00+-2.417e-05*a+8.193e-06*b+5.035e-11*a^2+-1.149e-11*b^2+
-1.836e-11*a*b+-4.871e-24*a^2*b^2+-5.131e-17*a^2*b+ 4.820e-17*b^2*a
#Predicción
zs0<-mu+t(lambda)%*%residuals(modelo4)
zs0## [,1]
## [1,] 1.705354Varianza del error de predicción
sz<-var(Randina$Velocidad_Viento)
G<-t(lambda)%*%k
VEP<-sz-G
VEP## [,1]
## [1,] 0.006857589Vemos que la Varianza del error de predicción es bastante baja (0.00685). Con ello se puede decir que la predicción dada es bastante buena. El reusltado obtenido es que en Palmira Valle del Cauca en la ubicación este: 355203.20 Norte: 391321.20, a las 6:00PM del día 12 de febrero de 2023 tuvo una velocidad del viento de 1.994 metros por segundo.
Kriging simple (espacial) en 100000 puntos no muestreados de la región Andina Colombiana
Se filtran los departamentos de la región Andina
# Filtrar los departamentos de la region andina
#install.packages("stringi")
library("stringi")
departamentos_sf <- st_read("F:/MI PC/Archivos/UNAL/2023-1/Espacial/Parcial1/parcial1/shp/MGN_DPTO_POLITICO.shp",
quiet=TRUE)
departamentos_sf$DPTO_CNMBR<-stringi::stri_trans_general(
c(departamentos_sf$DPTO_CNMBR), "Latin-ASCII")
Randina$Departamento<-stringi::stri_trans_general(
c(Randina$Departamento), "Latin-ASCII")
andina_sf <- departamentos_sf[departamentos_sf$DPTO_CNMBR %in%c(
Randina$Departamento),]Se filtran los municipios de los departamentos de la región Andina
## Deleting layer `andina_sfd' using driver `ESRI Shapefile'
## Writing layer `andina_sfd' to data source
## `andina_sfd.shp' using driver `ESRI Shapefile'
## Writing 791 features with 12 fields and geometry type Multi Polygon.Leyendo el polígono de la región Andina
## class : SpatialPoints
## features : 1
## extent : 452533.3, 452533.3, 639931.4, 639931.4 (xmin, xmax, ymin, ymax)
## crs : NAGrilla con los 100000 puntos no muestreados
muestra = spsample(poligonos, n = 100000, "regular")
muestra1 = as.data.frame(muestra)
names(muestra1) = c("x", "y")
gridded(muestra1) = c("x", "y")## Warning in points2grid(points, tolerance, round): grid has empty column/rows in
## dimension 2plot(muestra1)
Conversión a coordenasas UTM
deci_coord_muestra = SpatialPoints(cbind(muestra1$x, muestra1$y),proj4string = CRS("+proj=longlat"))
muestra_geod <- data.frame(deci_coord_muestra)
head(muestra_geod)## coords.x1 coords.x2
## 1 -76.91395 0.2510043
## 2 -76.89834 0.2510043
## 3 -76.88272 0.2510043
## 4 -76.86711 0.2510043
## 5 -76.85150 0.2510043
## 6 -76.83589 0.2510043library(rgdal)
library(sf)
library(stars)
library(gstat)
CRS_UTM_NY = (c("+proj=utm +zone=18 +datum=NAD83
+units=m +no_defs"))
utm_coord_muestra = spTransform(deci_coord_muestra, CRS(CRS_UTM_NY))
utm_coord_df_muestra = as.data.frame(utm_coord_muestra)
muestra_utm <- data.frame(utm_coord_df_muestra)
head(muestra_utm)## coords.x1 coords.x2
## 1 286987.9 27759.10
## 2 288725.9 27758.85
## 3 290463.9 27758.60
## 4 292201.9 27758.35
## 5 293939.8 27758.10
## 6 295677.8 27757.86Semivariograma teorico
ve.fit1exp = vgm(psill = 0.3726, model = "Exp", range = 327775.3900,
kappa = 0.5)Kriging espacial
Randina_utm2<-Randina_utm
Randina_utm2$z<-residuals(modelo4)
names(Randina_utm2)<-c("x","y","z")
coordinates(Randina_utm2)<-~x+y
names(muestra_utm)<-c("x","y")
coordinates(muestra_utm)<-~x+y
#Suponemos media conocida y constante Mapas de predicciones
#Suponemos media conocida y constante ya que se esta realizando
#el kriging con los residuales
beta=lm(Z~y + x + I(y^2) + I(x^2) + I(x*y) + I(x^2*y^2) + I(x^2*y)+
I(y^2*x))$coef
tprec.ok=krige(Z~y + x + I(y^2) + I(x^2) + I(x*y) + I(x^2*y^2) + I(x^2*y)+
I(y^2*x),Randina_utm2,muestra_utm, ve.fit1exp,beta=beta)## [using simple kriging]li = list("sp.polygons", poligonos)
pts = list("sp.points", Randina_utm2, pch = 20, col = "black", cex = 0.2)
spplot(tprec.ok, c("var1.pred"), as.table = TRUE,
sp.layout = list(li, pts), contour = FALSE, labels = FALSE,
pretty = TRUE, col = "black", col.regions = terrain.colors(100),
main = "Simple kriging predictions") 
spplot(tprec.ok, c("var1.var"), as.table = TRUE,
sp.layout = list(li, pts), contour = FALSE, labels = FALSE,
pretty = TRUE, col = "black", col.regions = terrain.colors(100),
main = "Simple kriging variance")
#writeOGR(tprec.ok, ".", "predicciones2", driver = "ESRI Shapefile")Kriging simple (espacio- tiempo) en un punto no muestreado del departamento de Caldas (perteneciente a la región Andina Colombiana donde la variable de interés sigue siendo la velocidad del viento
Descripción de los datos de la velocidad del viento en el departamento de Caldas
Lectura de los datos
library(readxl)
Caldas <- read_excel("C:/Users/Lenovo/Desktop/Caldas.xlsx")
attach(Caldas)
nrow(Caldas)## [1] 162Resumen 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.0summary(Caldas$Velocidad_Viento)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.500 1.000 1.300 1.269 1.500 2.300mean(Caldas$Velocidad_Viento)## [1] 1.268519Comportamiento 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)
abline(h=mean(Caldas$Velocidad_Viento),col="blue",lty=1)
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)
abline(h=mean(Caldas$Velocidad_Viento),col="blue",lty=1)
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 coportamiento un poco inusual en comparación con las demás horas mostradas las cuales si tienen comportamientos bastante similares.
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 23Caldasg <- 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 locationssummary(Caldasg)## Number of data points: 162
##
## Coordinates summary
## Este Norte
## min 403494.3 551127.1
## max 536730.8 620140.4
##
## Distance summary
## min max
## 0.0 140288.4
##
## Data summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.500000 1.000000 1.300000 1.268519 1.500000 2.300000
##
## Duplicated Coordinates
## dup Este Norte data
## 157 1 403494.3 559655.7 1.1
## 158 2 403494.3 559655.7 1.2
## 159 3 403494.3 559655.7 1.2
## 160 4 403494.3 559655.7 1.1
## 161 5 403494.3 559655.7 1.0
## 162 6 403494.3 559655.7 0.8
## 19 7 409993.4 552036.8 0.9
## 20 8 409993.4 552036.8 0.9
## 21 9 409993.4 552036.8 0.8
## 22 10 409993.4 552036.8 0.7
## 23 11 409993.4 552036.8 0.7
## 24 12 409993.4 552036.8 0.6
## 133 13 412205.1 561816.6 1.0
## 134 14 412205.1 561816.6 1.0
## 135 15 412205.1 561816.6 0.9
## 136 16 412205.1 561816.6 0.8
## 137 17 412205.1 561816.6 0.7
## 138 18 412205.1 561816.6 0.6
## 7 19 413081.6 578858.4 1.2
## 8 20 413081.6 578858.4 1.4
## 9 21 413081.6 578858.4 1.4
## 10 22 413081.6 578858.4 1.4
## 11 23 413081.6 578858.4 1.1
## 12 24 413081.6 578858.4 0.8
## 115 25 414963.2 570865.1 1.1
## 116 26 414963.2 570865.1 1.3
## 117 27 414963.2 570865.1 1.2
## 118 28 414963.2 570865.1 1.1
## 119 29 414963.2 570865.1 1.0
## 120 30 414963.2 570865.1 0.7
## 109 31 422223.2 599541.7 0.9
## 110 32 422223.2 599541.7 1.4
## 111 33 422223.2 599541.7 1.4
## 112 34 422223.2 599541.7 1.3
## 113 35 422223.2 599541.7 1.2
## 114 36 422223.2 599541.7 1.0
## 139 37 428031.8 602095.7 0.8
## 140 38 428031.8 602095.7 1.4
## 141 39 428031.8 602095.7 1.4
## 142 40 428031.8 602095.7 1.3
## 143 41 428031.8 602095.7 1.2
## 144 42 428031.8 602095.7 1.1
## 97 43 430776.0 554669.2 0.8
## 98 44 430776.0 554669.2 1.3
## 99 45 430776.0 554669.2 1.3
## 100 46 430776.0 554669.2 1.3
## 101 47 430776.0 554669.2 1.3
## 102 48 430776.0 554669.2 1.1
## 25 49 432655.3 551127.1 0.8
## 26 50 432655.3 551127.1 1.3
## 27 51 432655.3 551127.1 1.4
## 28 52 432655.3 551127.1 1.3
## 29 53 432655.3 551127.1 1.3
## 30 54 432655.3 551127.1 1.2
## 61 55 433509.0 605093.4 0.8
## 62 56 433509.0 605093.4 1.3
## 63 57 433509.0 605093.4 1.4
## 64 58 433509.0 605093.4 1.3
## 65 59 433509.0 605093.4 1.2
## 66 60 433509.0 605093.4 1.1
## 31 61 437675.7 585533.5 0.8
## 32 62 437675.7 585533.5 1.6
## 33 63 437675.7 585533.5 1.7
## 34 64 437675.7 585533.5 1.6
## 35 65 437675.7 585533.5 1.6
## 36 66 437675.7 585533.5 1.4
## 43 67 439447.1 596506.1 0.7
## 44 68 439447.1 596506.1 1.5
## 45 69 439447.1 596506.1 1.5
## 46 70 439447.1 596506.1 1.5
## 47 71 439447.1 596506.1 1.4
## 48 72 439447.1 596506.1 1.3
## 79 73 442387.2 571115.3 0.7
## 80 74 442387.2 571115.3 1.8
## 81 75 442387.2 571115.3 1.9
## 82 76 442387.2 571115.3 1.9
## 83 77 442387.2 571115.3 1.8
## 84 78 442387.2 571115.3 1.7
## 151 79 443050.4 557529.1 0.8
## 152 80 443050.4 557529.1 1.9
## 153 81 443050.4 557529.1 2.0
## 154 82 443050.4 557529.1 2.0
## 155 83 443050.4 557529.1 2.0
## 156 84 443050.4 557529.1 1.8
## 13 85 445560.8 582663.5 0.7
## 14 86 445560.8 582663.5 1.6
## 15 87 445560.8 582663.5 1.7
## 16 88 445560.8 582663.5 1.7
## 17 89 445560.8 582663.5 1.7
## 18 90 445560.8 582663.5 1.6
## 49 91 445569.4 559061.5 0.8
## 50 92 445569.4 559061.5 1.9
## 51 93 445569.4 559061.5 2.1
## 52 94 445569.4 559061.5 2.1
## 53 95 445569.4 559061.5 2.0
## 54 96 445569.4 559061.5 1.8
## 121 97 446023.5 597235.3 0.7
## 122 98 446023.5 597235.3 1.4
## 123 99 446023.5 597235.3 1.5
## 124 100 446023.5 597235.3 1.4
## 125 101 446023.5 597235.3 1.4
## 126 102 446023.5 597235.3 1.5
## 91 103 449091.1 610957.6 0.8
## 92 104 449091.1 610957.6 1.3
## 93 105 449091.1 610957.6 1.3
## 94 106 449091.1 610957.6 1.2
## 95 107 449091.1 610957.6 1.1
## 96 108 449091.1 610957.6 1.2
## 1 109 449624.5 620140.4 1.0
## 2 110 449624.5 620140.4 1.4
## 3 111 449624.5 620140.4 1.3
## 4 112 449624.5 620140.4 1.3
## 5 113 449624.5 620140.4 1.2
## 6 114 449624.5 620140.4 1.2
## 73 115 471256.8 584100.4 0.7
## 74 116 471256.8 584100.4 0.5
## 75 117 471256.8 584100.4 0.6
## 76 118 471256.8 584100.4 0.5
## 77 119 471256.8 584100.4 0.7
## 78 120 471256.8 584100.4 1.0
## 103 121 482240.8 595034.0 0.8
## 104 122 482240.8 595034.0 0.9
## 105 123 482240.8 595034.0 1.0
## 106 124 482240.8 595034.0 1.0
## 107 125 482240.8 595034.0 1.3
## 108 126 482240.8 595034.0 1.4
## 55 127 483064.1 580932.1 0.9
## 56 128 483064.1 580932.1 0.9
## 57 129 483064.1 580932.1 1.1
## 58 130 483064.1 580932.1 1.1
## 59 131 483064.1 580932.1 1.3
## 60 132 483064.1 580932.1 1.5
## 67 133 494117.0 585552.7 1.0
## 68 134 494117.0 585552.7 1.5
## 69 135 494117.0 585552.7 1.7
## 70 136 494117.0 585552.7 1.8
## 71 137 494117.0 585552.7 1.9
## 72 138 494117.0 585552.7 2.0
## 127 139 500859.2 598325.9 0.7
## 128 140 500859.2 598325.9 1.7
## 129 141 500859.2 598325.9 1.9
## 130 142 500859.2 598325.9 2.1
## 131 143 500859.2 598325.9 2.3
## 132 144 500859.2 598325.9 2.2
## 145 145 509843.7 587752.1 1.2
## 146 146 509843.7 587752.1 1.5
## 147 147 509843.7 587752.1 1.6
## 148 148 509843.7 587752.1 1.6
## 149 149 509843.7 587752.1 1.7
## 150 150 509843.7 587752.1 1.8
## 85 151 512233.2 616202.9 0.7
## 86 152 512233.2 616202.9 1.3
## 87 153 512233.2 616202.9 1.4
## 88 154 512233.2 616202.9 1.5
## 89 155 512233.2 616202.9 1.5
## 90 156 512233.2 616202.9 1.5
## 37 157 536730.8 603574.9 1.2
## 38 158 536730.8 603574.9 1.0
## 39 159 536730.8 603574.9 0.9
## 40 160 536730.8 603574.9 1.0
## 41 161 536730.8 603574.9 0.8
## 42 162 536730.8 603574.9 0.9Explorando 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)
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",
col.grid="lightblue", main="Tendencia", pch=20)
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.0007688par(mfrow=c(2,2))
plot(Este2,residuals(mod1),pch=16,col=2,ylab="residuales",
main="Este vs residuales")
plot(Norte2,residuals(mod1),pch=16,col=3,ylab="residuales",
main="Norte vs residuales")
plot(Hora2, residuals(mod1),pch=16,col=4,ylab="residuales",
main="Hora vs residuales")
hist(residuals(mod1),col=2,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-06par(mfrow=c(2,2))
plot(Este2,residuals(mod2),pch=16,col=2,ylab="residuales",
main="Este vs residuales")
plot(Norte2,residuals(mod2),pch=16,col=3,ylab="residuales",
main="Norte vs residuales")
plot(Hora2, residuals(mod2),pch=16,col=4,ylab="residuales",
main="Hora vs residuales")
hist(residuals(mod2),col=2,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-07par(mfrow=c(2,2))
plot(Este2,residuals(mod3),pch=16,col=2,ylab="residuales",
main="Este vs residuales")
plot(Norte2,residuals(mod3),pch=16,col=3,ylab="residuales",
main="Norte vs residuales")
plot(Hora2, residuals(mod3),pch=16,col=4,ylab="residuales",
main="Hora vs residuales")
hist(residuals(mod3),col=2,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-10par(mfrow=c(2,2))
plot(Este2,residuals(mod4),pch=16,col=2,ylab="residuales",
main="Este vs residuales")
plot(Norte2,residuals(mod4),pch=16,col=3,ylab="residuales",
main="Norte vs residuales")
plot(Hora2, residuals(mod4),pch=16,col=4,ylab="residuales",
main="Hora vs residuales")
hist(residuals(mod4),col=2,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-05par(mfrow=c(2,2))
plot(Este2,residuals(mod5),pch=16,col=2,ylab="residuales",
main="Este vs residuales")
plot(Norte2,residuals(mod5),pch=16,col=3,ylab="residuales",
main="Norte vs residuales")
plot(Hora2, residuals(mod5),pch=16,col=4,ylab="residuales",
main="Hora vs residuales")
hist(residuals(mod5),col=2,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-11par(mfrow=c(2,2))
plot(Este2,residuals(mod6),pch=16,col=2,ylab="residuales",
main="Este vs residuales")
plot(Norte2,residuals(mod6),pch=16,col=3,ylab="residuales",
main="Norte vs residuales")
plot(Hora2, residuals(mod6),pch=16,col=4,ylab="residuales",
main="Hora vs residuales")
hist(residuals(mod6),col=2,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.001389par(mfrow=c(2,2))
plot(Este2,residuals(mod7),pch=16,col=2,ylab="residuales",
main="Este vs residuales")
plot(Norte2,residuals(mod7),pch=16,col=3,ylab="residuales",
main="Norte vs residuales")
plot(Hora2, residuals(mod7),pch=16,col=4,ylab="residuales",
main="Hora vs residuales")
hist(residuals(mod7),col=2,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=2,ylab="residuales",
main="Este vs residuales")
plot(Norte2,residuals(mod8),pch=16,col=3,ylab="residuales",
main="Norte vs residuales")
plot(Hora2, residuals(mod8),pch=16,col=4,ylab="residuales",
main="Hora vs residuales")
hist(residuals(mod8),col=2,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.00173par(mfrow=c(2,2))
plot(Este2,residuals(mod9),pch=16,col=2,ylab="residuales",
main="Este vs residuales")
plot(Norte2,residuals(mod9),pch=16,col=3,ylab="residuales",
main="Norte vs residuales")
plot(Hora2, residuals(mod9),pch=16,col=4,ylab="residuales",
main="Hora vs residuales")
hist(residuals(mod9),col=2,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-11par(mfrow=c(2,2))
plot(Este2,residuals(mod10),pch=16,col=2,ylab="residuales",
main="Este vs residuales")
plot(Norte2,residuals(mod10),pch=16,col=3,ylab="residuales",
main="Norte vs residuales")
plot(Hora2, residuals(mod10),pch=16,col=4,ylab="residuales",
main="Hora vs residuales")
hist(residuals(mod10),col=2,main="Residuales" )
Finalmente escogemos el modelo 10 al considerar que corrije mejor la tendencia.
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 variogramhead(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 TRUEPara 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 variogramhead(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 TRUEPara 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 variogramhead(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 TRUEPara 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 variogramhead(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 TRUEPara 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 variogramhead(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 TRUEPara 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 variogramhead(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 TRUEComparació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")
Kriging espacio tiempo
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]))
head(matDistT)## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
## 1 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4
## 2 1 0 1 2 3 4 1 0 1 2 3 4 1 0 1 2 3 4 1 0 1 2 3 4 1 0 1 2 3
## 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2
## 4 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1
## 5 4 3 2 1 0 1 4 3 2 1 0 1 4 3 2 1 0 1 4 3 2 1 0 1 4 3 2 1 0
## 6 5 4 3 2 1 0 5 4 3 2 1 0 5 4 3 2 1 0 5 4 3 2 1 0 5 4 3 2 1
## 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
## 1 5 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5 0
## 2 4 1 0 1 2 3 4 1 0 1 2 3 4 1 0 1 2 3 4 1 0 1 2 3 4 1
## 3 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2
## 4 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3
## 5 1 4 3 2 1 0 1 4 3 2 1 0 1 4 3 2 1 0 1 4 3 2 1 0 1 4
## 6 0 5 4 3 2 1 0 5 4 3 2 1 0 5 4 3 2 1 0 5 4 3 2 1 0 5
## 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
## 1 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2
## 2 0 1 2 3 4 1 0 1 2 3 4 1 0 1 2 3 4 1 0 1 2 3 4 1 0 1
## 3 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0
## 4 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1
## 5 3 2 1 0 1 4 3 2 1 0 1 4 3 2 1 0 1 4 3 2 1 0 1 4 3 2
## 6 4 3 2 1 0 5 4 3 2 1 0 5 4 3 2 1 0 5 4 3 2 1 0 5 4 3
## 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
## 1 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2
## 2 2 3 4 1 0 1 2 3 4 1 0 1 2 3 4 1 0 1 2 3 4 1 0 1
## 3 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0
## 4 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1
## 5 1 0 1 4 3 2 1 0 1 4 3 2 1 0 1 4 3 2 1 0 1 4 3 2
## 6 2 1 0 5 4 3 2 1 0 5 4 3 2 1 0 5 4 3 2 1 0 5 4 3
## 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
## 1 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3
## 2 2 3 4 1 0 1 2 3 4 1 0 1 2 3 4 1 0 1 2
## 3 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1
## 4 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0
## 5 1 0 1 4 3 2 1 0 1 4 3 2 1 0 1 4 3 2 1
## 6 2 1 0 5 4 3 2 1 0 5 4 3 2 1 0 5 4 3 2
## 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143
## 1 4 5 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4
## 2 3 4 1 0 1 2 3 4 1 0 1 2 3 4 1 0 1 2 3
## 3 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2
## 4 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1
## 5 0 1 4 3 2 1 0 1 4 3 2 1 0 1 4 3 2 1 0
## 6 1 0 5 4 3 2 1 0 5 4 3 2 1 0 5 4 3 2 1
## 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162
## 1 5 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5
## 2 4 1 0 1 2 3 4 1 0 1 2 3 4 1 0 1 2 3 4
## 3 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2 3
## 4 2 3 2 1 0 1 2 3 2 1 0 1 2 3 2 1 0 1 2
## 5 1 4 3 2 1 0 1 4 3 2 1 0 1 4 3 2 1 0 1
## 6 0 5 4 3 2 1 0 5 4 3 2 1 0 5 4 3 2 1 0head(matDistSp)## 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## 1 0 0 0 0 0 0 55132.4 55132.4 55132.4 55132.4 55132.4 55132.4 37696.59 37696.59
## 2 0 0 0 0 0 0 55132.4 55132.4 55132.4 55132.4 55132.4 55132.4 37696.59 37696.59
## 3 0 0 0 0 0 0 55132.4 55132.4 55132.4 55132.4 55132.4 55132.4 37696.59 37696.59
## 4 0 0 0 0 0 0 55132.4 55132.4 55132.4 55132.4 55132.4 55132.4 37696.59 37696.59
## 5 0 0 0 0 0 0 55132.4 55132.4 55132.4 55132.4 55132.4 55132.4 37696.59 37696.59
## 6 0 0 0 0 0 0 55132.4 55132.4 55132.4 55132.4 55132.4 55132.4 37696.59 37696.59
## 15 16 17 18 19 20 21 22
## 1 37696.59 37696.59 37696.59 37696.59 78795.42 78795.42 78795.42 78795.42
## 2 37696.59 37696.59 37696.59 37696.59 78795.42 78795.42 78795.42 78795.42
## 3 37696.59 37696.59 37696.59 37696.59 78795.42 78795.42 78795.42 78795.42
## 4 37696.59 37696.59 37696.59 37696.59 78795.42 78795.42 78795.42 78795.42
## 5 37696.59 37696.59 37696.59 37696.59 78795.42 78795.42 78795.42 78795.42
## 6 37696.59 37696.59 37696.59 37696.59 78795.42 78795.42 78795.42 78795.42
## 23 24 25 26 27 28 29 30
## 1 78795.42 78795.42 71068.87 71068.87 71068.87 71068.87 71068.87 71068.87
## 2 78795.42 78795.42 71068.87 71068.87 71068.87 71068.87 71068.87 71068.87
## 3 78795.42 78795.42 71068.87 71068.87 71068.87 71068.87 71068.87 71068.87
## 4 78795.42 78795.42 71068.87 71068.87 71068.87 71068.87 71068.87 71068.87
## 5 78795.42 78795.42 71068.87 71068.87 71068.87 71068.87 71068.87 71068.87
## 6 78795.42 78795.42 71068.87 71068.87 71068.87 71068.87 71068.87 71068.87
## 31 32 33 34 35 36 37 38 39
## 1 36611.64 36611.64 36611.64 36611.64 36611.64 36611.64 88667.5 88667.5 88667.5
## 2 36611.64 36611.64 36611.64 36611.64 36611.64 36611.64 88667.5 88667.5 88667.5
## 3 36611.64 36611.64 36611.64 36611.64 36611.64 36611.64 88667.5 88667.5 88667.5
## 4 36611.64 36611.64 36611.64 36611.64 36611.64 36611.64 88667.5 88667.5 88667.5
## 5 36611.64 36611.64 36611.64 36611.64 36611.64 36611.64 88667.5 88667.5 88667.5
## 6 36611.64 36611.64 36611.64 36611.64 36611.64 36611.64 88667.5 88667.5 88667.5
## 40 41 42 43 44 45 46 47 48
## 1 88667.5 88667.5 88667.5 25732.45 25732.45 25732.45 25732.45 25732.45 25732.45
## 2 88667.5 88667.5 88667.5 25732.45 25732.45 25732.45 25732.45 25732.45 25732.45
## 3 88667.5 88667.5 88667.5 25732.45 25732.45 25732.45 25732.45 25732.45 25732.45
## 4 88667.5 88667.5 88667.5 25732.45 25732.45 25732.45 25732.45 25732.45 25732.45
## 5 88667.5 88667.5 88667.5 25732.45 25732.45 25732.45 25732.45 25732.45 25732.45
## 6 88667.5 88667.5 88667.5 25732.45 25732.45 25732.45 25732.45 25732.45 25732.45
## 49 50 51 52 53 54 55 56
## 1 61213.34 61213.34 61213.34 61213.34 61213.34 61213.34 51531.56 51531.56
## 2 61213.34 61213.34 61213.34 61213.34 61213.34 61213.34 51531.56 51531.56
## 3 61213.34 61213.34 61213.34 61213.34 61213.34 61213.34 51531.56 51531.56
## 4 61213.34 61213.34 61213.34 61213.34 61213.34 61213.34 51531.56 51531.56
## 5 61213.34 61213.34 61213.34 61213.34 61213.34 61213.34 51531.56 51531.56
## 6 61213.34 61213.34 61213.34 61213.34 61213.34 61213.34 51531.56 51531.56
## 57 58 59 60 61 62 63 64
## 1 51531.56 51531.56 51531.56 51531.56 22048.12 22048.12 22048.12 22048.12
## 2 51531.56 51531.56 51531.56 51531.56 22048.12 22048.12 22048.12 22048.12
## 3 51531.56 51531.56 51531.56 51531.56 22048.12 22048.12 22048.12 22048.12
## 4 51531.56 51531.56 51531.56 51531.56 22048.12 22048.12 22048.12 22048.12
## 5 51531.56 51531.56 51531.56 51531.56 22048.12 22048.12 22048.12 22048.12
## 6 51531.56 51531.56 51531.56 51531.56 22048.12 22048.12 22048.12 22048.12
## 65 66 67 68 69 70 71 72
## 1 22048.12 22048.12 56355.12 56355.12 56355.12 56355.12 56355.12 56355.12
## 2 22048.12 22048.12 56355.12 56355.12 56355.12 56355.12 56355.12 56355.12
## 3 22048.12 22048.12 56355.12 56355.12 56355.12 56355.12 56355.12 56355.12
## 4 22048.12 22048.12 56355.12 56355.12 56355.12 56355.12 56355.12 56355.12
## 5 22048.12 22048.12 56355.12 56355.12 56355.12 56355.12 56355.12 56355.12
## 6 22048.12 22048.12 56355.12 56355.12 56355.12 56355.12 56355.12 56355.12
## 73 74 75 76 77 78 79 80 81
## 1 42033.8 42033.8 42033.8 42033.8 42033.8 42033.8 49556.39 49556.39 49556.39
## 2 42033.8 42033.8 42033.8 42033.8 42033.8 42033.8 49556.39 49556.39 49556.39
## 3 42033.8 42033.8 42033.8 42033.8 42033.8 42033.8 49556.39 49556.39 49556.39
## 4 42033.8 42033.8 42033.8 42033.8 42033.8 42033.8 49556.39 49556.39 49556.39
## 5 42033.8 42033.8 42033.8 42033.8 42033.8 42033.8 49556.39 49556.39 49556.39
## 6 42033.8 42033.8 42033.8 42033.8 42033.8 42033.8 49556.39 49556.39 49556.39
## 82 83 84 85 86 87 88 89
## 1 49556.39 49556.39 49556.39 62732.49 62732.49 62732.49 62732.49 62732.49
## 2 49556.39 49556.39 49556.39 62732.49 62732.49 62732.49 62732.49 62732.49
## 3 49556.39 49556.39 49556.39 62732.49 62732.49 62732.49 62732.49 62732.49
## 4 49556.39 49556.39 49556.39 62732.49 62732.49 62732.49 62732.49 62732.49
## 5 49556.39 49556.39 49556.39 62732.49 62732.49 62732.49 62732.49 62732.49
## 6 49556.39 49556.39 49556.39 62732.49 62732.49 62732.49 62732.49 62732.49
## 90 91 92 93 94 95 96 97
## 1 62732.49 9198.286 9198.286 9198.286 9198.286 9198.286 9198.286 68130.33
## 2 62732.49 9198.286 9198.286 9198.286 9198.286 9198.286 9198.286 68130.33
## 3 62732.49 9198.286 9198.286 9198.286 9198.286 9198.286 9198.286 68130.33
## 4 62732.49 9198.286 9198.286 9198.286 9198.286 9198.286 9198.286 68130.33
## 5 62732.49 9198.286 9198.286 9198.286 9198.286 9198.286 9198.286 68130.33
## 6 62732.49 9198.286 9198.286 9198.286 9198.286 9198.286 9198.286 68130.33
## 98 99 100 101 102 103 104 105
## 1 68130.33 68130.33 68130.33 68130.33 68130.33 41160.12 41160.12 41160.12
## 2 68130.33 68130.33 68130.33 68130.33 68130.33 41160.12 41160.12 41160.12
## 3 68130.33 68130.33 68130.33 68130.33 68130.33 41160.12 41160.12 41160.12
## 4 68130.33 68130.33 68130.33 68130.33 68130.33 41160.12 41160.12 41160.12
## 5 68130.33 68130.33 68130.33 68130.33 68130.33 41160.12 41160.12 41160.12
## 6 68130.33 68130.33 68130.33 68130.33 68130.33 41160.12 41160.12 41160.12
## 106 107 108 109 110 111 112 113
## 1 41160.12 41160.12 41160.12 34280.19 34280.19 34280.19 34280.19 34280.19
## 2 41160.12 41160.12 41160.12 34280.19 34280.19 34280.19 34280.19 34280.19
## 3 41160.12 41160.12 41160.12 34280.19 34280.19 34280.19 34280.19 34280.19
## 4 41160.12 41160.12 41160.12 34280.19 34280.19 34280.19 34280.19 34280.19
## 5 41160.12 41160.12 41160.12 34280.19 34280.19 34280.19 34280.19 34280.19
## 6 41160.12 41160.12 41160.12 34280.19 34280.19 34280.19 34280.19 34280.19
## 114 115 116 117 118 119 120 121
## 1 34280.19 60244.92 60244.92 60244.92 60244.92 60244.92 60244.92 23186.42
## 2 34280.19 60244.92 60244.92 60244.92 60244.92 60244.92 60244.92 23186.42
## 3 34280.19 60244.92 60244.92 60244.92 60244.92 60244.92 60244.92 23186.42
## 4 34280.19 60244.92 60244.92 60244.92 60244.92 60244.92 60244.92 23186.42
## 5 34280.19 60244.92 60244.92 60244.92 60244.92 60244.92 60244.92 23186.42
## 6 34280.19 60244.92 60244.92 60244.92 60244.92 60244.92 60244.92 23186.42
## 122 123 124 125 126 127 128 129
## 1 23186.42 23186.42 23186.42 23186.42 23186.42 55685.47 55685.47 55685.47
## 2 23186.42 23186.42 23186.42 23186.42 23186.42 55685.47 55685.47 55685.47
## 3 23186.42 23186.42 23186.42 23186.42 23186.42 55685.47 55685.47 55685.47
## 4 23186.42 23186.42 23186.42 23186.42 23186.42 55685.47 55685.47 55685.47
## 5 23186.42 23186.42 23186.42 23186.42 23186.42 55685.47 55685.47 55685.47
## 6 23186.42 23186.42 23186.42 23186.42 23186.42 55685.47 55685.47 55685.47
## 130 131 132 133 134 135 136 137
## 1 55685.47 55685.47 55685.47 69295.54 69295.54 69295.54 69295.54 69295.54
## 2 55685.47 55685.47 55685.47 69295.54 69295.54 69295.54 69295.54 69295.54
## 3 55685.47 55685.47 55685.47 69295.54 69295.54 69295.54 69295.54 69295.54
## 4 55685.47 55685.47 55685.47 69295.54 69295.54 69295.54 69295.54 69295.54
## 5 55685.47 55685.47 55685.47 69295.54 69295.54 69295.54 69295.54 69295.54
## 6 55685.47 55685.47 55685.47 69295.54 69295.54 69295.54 69295.54 69295.54
## 138 139 140 141 142 143 144 145
## 1 69295.54 28139.85 28139.85 28139.85 28139.85 28139.85 28139.85 68376.64
## 2 69295.54 28139.85 28139.85 28139.85 28139.85 28139.85 28139.85 68376.64
## 3 69295.54 28139.85 28139.85 28139.85 28139.85 28139.85 28139.85 68376.64
## 4 69295.54 28139.85 28139.85 28139.85 28139.85 28139.85 28139.85 68376.64
## 5 69295.54 28139.85 28139.85 28139.85 28139.85 28139.85 28139.85 68376.64
## 6 69295.54 28139.85 28139.85 28139.85 28139.85 28139.85 28139.85 68376.64
## 146 147 148 149 150 151 152 153
## 1 68376.64 68376.64 68376.64 68376.64 68376.64 62955.46 62955.46 62955.46
## 2 68376.64 68376.64 68376.64 68376.64 68376.64 62955.46 62955.46 62955.46
## 3 68376.64 68376.64 68376.64 68376.64 68376.64 62955.46 62955.46 62955.46
## 4 68376.64 68376.64 68376.64 68376.64 68376.64 62955.46 62955.46 62955.46
## 5 68376.64 68376.64 68376.64 68376.64 68376.64 62955.46 62955.46 62955.46
## 6 68376.64 68376.64 68376.64 68376.64 68376.64 62955.46 62955.46 62955.46
## 154 155 156 157 158 159 160 161
## 1 62955.46 62955.46 62955.46 76068.32 76068.32 76068.32 76068.32 76068.32
## 2 62955.46 62955.46 62955.46 76068.32 76068.32 76068.32 76068.32 76068.32
## 3 62955.46 62955.46 62955.46 76068.32 76068.32 76068.32 76068.32 76068.32
## 4 62955.46 62955.46 62955.46 76068.32 76068.32 76068.32 76068.32 76068.32
## 5 62955.46 62955.46 62955.46 76068.32 76068.32 76068.32 76068.32 76068.32
## 6 62955.46 62955.46 62955.46 76068.32 76068.32 76068.32 76068.32 76068.32
## 162
## 1 76068.32
## 2 76068.32
## 3 76068.32
## 4 76068.32
## 5 76068.32
## 6 76068.32Kriging 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.3cressie1=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(mod6)
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 hessianaCá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=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))
mu1## 1
## 1.363263z_pred0=t(lambda1)%*%datos1[,4]+mu1
z_pred0## [,1]
## [1,] 1.411814Error cuadrático medio
ECM1<-t(sigma0[162,-162])%*%sigma1%*%sigma0[162,-162]
ECM1## [,1]
## [1,] 0.4690119Kriging espacio tiempo con el modelo C R E S S I E - H U A N G (1999)
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 hessianaCá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))
mu2## 1
## 1.363263z_pred02=t(lambda2)%*%datos1[,4]+mu2
z_pred02## [,1]
## [1,] 1.363263Error cuadrático medio
ECM2<-t(sigma02[162,-162]) %*% sigma2 %*% sigma02[162,-162]
ECM2## [,1]
## [1,] 0.001812986Comparación entre las varianzas del error de predicción entre los modelos propuestos
| Modelo | ECM obtenido |
|---|---|
| CRESSIE-HUANG | 0.001812986 |
| 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.
Kriging espacio tiempo con datos funcionales
Los datos funcionales son un tipo de datos en el que cada observación se representa como una función o curva en lugar de un valor único. El análisis de datos funcionales se basa en el estudio de la función que describe la variabilidad de un conjunto de datos en un espacio de muestras1. Esto permite aprovechar la estructura y dependencias presentes en los datos para realizar análisis y predicciones más precisas.
Análisis gráfico de los datos funcionales espaciales
Veamos un resumen de cada variable que compone la base de datos:
## Warning: package 'sqldf' was built under R version 4.2.2## Warning: package 'gsubfn' was built under R version 4.2.2## Warning: package 'proto' was built under R version 4.2.2## Warning: package 'RSQLite' was built under R version 4.2.2## Warning: package 'fda' was built under R version 4.2.3## Warning: package 'fds' was built under R version 4.2.3## Warning: package 'rainbow' was built under R version 4.2.3## Warning: package 'pcaPP' was built under R version 4.2.3## Warning: package 'RCurl' was built under R version 4.2.3## Warning: package 'deSolve' was built under R version 4.2.3## Warning: package 'fda.usc' was built under R version 4.2.3## Warning: package 'nlme' was built under R version 4.2.2## Warning: package 'xtable' was built under R version 4.2.2## Latitud Longitud Region Departamento
## Min. :4.986 Min. :-75.87 Length:2241 Length:2241
## 1st Qu.:5.082 1st Qu.:-75.65 Class :character Class :character
## Median :5.297 Median :-75.51 Mode :character Mode :character
## Mean :5.283 Mean :-75.44
## 3rd Qu.:5.424 3rd Qu.:-75.16
## Max. :5.610 Max. :-74.67
## Municipio Fecha Hora Velocidad_Viento
## Length:2241 Length:2241 Length:2241 Min. :0.0000
## Class :character Class :character Class :character 1st Qu.:0.6000
## Mode :character Mode :character Mode :character Median :1.0000
## Mean :0.9963
## 3rd Qu.:1.3000
## Max. :2.5000Veamos los municipios a los que se les hico la medición de la velocidad del viento desde las 0:00 hasta las 23:00
## [1] "AGUADAS" "ANSERMA" "ARANZAZU" "BELALCAZAR" "CHINCHINA"
## [6] "FILADELFIA" "LA DORADA" "LA MERCED" "MANIZALES" "MANZANARES"
## [11] "MARMATO" "MARQUETALIA" "MARULANDA" "NEIRA" "NORCASIA"
## [16] "PACORA" "PALESTINA" "PENSILVANIA" "RIOSUCIO" "RISARALDA"
## [21] "SALAMINA" "SAMANA" "SAN JOSE" "SUPIA" "VICTORIA"
## [26] "VILLAMARIA" "VITERBO"Gráfico descriptivo de la variable velocidad del viento en el tiempo por cada ubicación espacial
## named integer(0)Método Bag
par(mfrow=c(1,2))
fboxplot(data= TempTs, plot.type = "functional", type = "bag", projmethod = "PCAproj")## Warning in fts(x, notchlow): Please assign column name for the data matrix.## Warning in fts(x, notchupper): Please assign column name for the data matrix.## Warning in fts(x, centercurve): Please assign column name for the data matrix.## Warning in fts(x, outliercurve): Please assign column name for the data matrix.fboxplot(data= TempTs, plot.type = "bivariate", type = "bag", projmethod = "PCAproj",
ylim = c(-3.5,2), xlim = c(-3,3))
Hay curva atipica en el municipio 7
Metodo HDR
par(mfrow=c(1,2))
fboxplot(data= TempTs, plot.type = "functional", type = "hdr", projmethod = "PCAproj")## Warning in fts(x, centercurve): Please assign column name for the data matrix.fboxplot(data= TempTs, plot.type = "bivariate", type = "hdr", projmethod = "PCAproj")
Vemos que se detecta una curva atipica en el municipio 7 y 22.
Método ambos
par(mfrow=c(1,1))
fbplot(fit= C, method = "BD2", xlab = "Hora", ylab = "Velocidad Viento", outliercol = 2)
## $depth
## [1] 0.08760684 0.07763533 0.10968661 0.07407407 0.09971510 0.18803419
## [7] 0.07407407 0.13105413 0.07407407 0.09401709 0.07407407 0.09116809
## [13] 0.07407407 0.07905983 0.11111111 0.14886040 0.13888889 0.08404558
## [19] 0.07763533 0.07407407 0.09259259 0.07407407 0.07407407 0.11111111
## [25] 0.07407407 0.07549858 0.08689459
##
## $outpoint
## integer(0)
##
## $medcurve
## [1] 6Este metodo detecta outliers.
Al comparar la detección de datos atìpicos utilizando las tres herramientas mencionadas (bagplot, hdr plot y fbplot) en la muestra de las 27 curvas de velocidad del viento, se evidencia que el hdr plot resulta ser la herramienta mas estricta o exigente en cuanto a que detecta mejor los outliers.
Otras medidas resumen
El parámetro “trim” especifica la fracción de datos que se eliminaran.Alpha representa el nivel de confianza utilizado en ciertos cálculos de profundidad funcional.Trimmed mean, es una medida de tendencia central que se calcula eliminando un cierto porcentaje de valores extremos o atìpicos en los extremos de un conjunto de datos antes de calcular la media.
par(mfrow=c(1,1))
###Calcula la profundidad de Fraiman and Muniz (2001)
MedFM = fdepth(data = TempTs, type = "FM", trim = 0.1)
plot(MedFM)
Este tipo de profundidad funcional utiliza la distancia de Mahalanobis para calcular la profundidad de un punto de datos en relación con la distribución multivariante general de los datos. Los puntos que tienen una mayor profundidad funcional están más cerca del centro de la distribución y se consideran más “típicos”.
###Calcula la profundidad a traves de proyecciones aletorias de Cuevas et al. (2007)
MedRP = fdepth(data = TempTs, type = "RP", trim = 0.1)
plot(MedRP)
Esta medida de profundidad funcional se basa en proyectar los datos en una dirección aleatoria y evaluar la posición relativa de los puntos de datos en función de la dispersión de las proyecciones.
MedRPD = fdepth(data = TempTs, type = "RPD", trim = 0.1)## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.plot(MedRPD)
Similar al RP depth, pero específico para datos funcionales.
Parametros estimados - Moda
MODA = fdepth(data = TempTs, type = "mode", trim = 0.1)
plot(MODA)
Esta medida de profundidad se basa en el concepto de encontrar los modos de las curvas de contorno generadas a partir de los datos. Los puntos que se encuentran en las regiones más densas de los datos tienen una mayor profundidad funcional.
RADIUS = fdepth(data = TempTs, type = "radius", trim = 0.1, alpha=0.05, weight = "hard")
plot(RADIUS)
RADIUS = fdepth(data = TempTs, type = "radius", trim = 0.1, alpha=0.05, weight = "soft")
plot(RADIUS)
Esta medida calcula la profundidad de un punto en relación con el radio de un contenedor convexo que contiene la mayoría de los datos. Los puntos que se encuentran màs cerca del centro del contenedor convexo tienen una mayor profundidad funcional.
Finalmemnte, a través de procedimientos bootstrap es posible estimar intervalos de confianza para las medidas de resumen.
Boostrap Funcional
Este método se basa en la idea de generar m?ltiples muestras (con reemplazo) a partir de la muestra original y calcular las estimaciones de inter?s en cada una de estas muestras.
BSpl <- create.bspline.basis(rangeval=c(0,3),nbasis=10, norder=4)
Ftemp <- Data2fd(C,basisobj=BSpl,argvals=seq( 0 , 3 , length= 24 ))
plot(Ftemp, main="Datos suavizados de la Velocidad del viento", xlab = "Horas", ylab = "Velocidad del Viento")
## [1] "done"out.boot1=fdata.bootstrap(Ftemp,statistic=func.mean,nb=200,draw=TRUE)
out.boot2=fdata.bootstrap(Ftemp,statistic=func.trim.FM,nb=100,draw=TRUE)
out.boot2=fdata.bootstrap(Ftemp,statistic=func.trim.mode,nb=100,draw=TRUE, trim = 0.1)
out.boot2=fdata.bootstrap(Ftemp,statistic=func.trim.RP,nb=100,draw=TRUE, trim = 0.1)
ACPF
PCA para datos funcionales
El Análisis de Componentes Principales (PCA) es una técnica estadística utilizada para reducir la dimensionalidad de un conjunto de datos y encontrar la variabilidad en el mismo1. En el caso de los datos funcionales, el objetivo del análisis en componentes principales funcionales es reducir la dimensión del espacio, ya que las variables aleatorias funcionales se encuentran en espacios de dimensión infinita2. Esto permite representar los datos de manera más compacta y eficiente, facilitando su análisis y modelado.
par(mfrow=c(3,3))
PCA = pca.fd(Ftemp, 9)
plot(PCA)
Tomaremos una varianza del 90%. Con dos funciones propias se explica este pocentaje de la varianza.
Gráficos de los scores
par(mfrow=c(1,1))
plot(x = PCA$scores[,1], y = PCA$scores[,2])
text(x = PCA$scores[,1], y = PCA$scores[,2] , seq(14,16,1), cex = 1.5)
En la literatura existen varias funciones base que permiten desarrollar la expansión. La selección del tipo de funciones de la base depende de las características que cumpla que fenómeno de estudio. Por ejemplo, las series de Fourier se utilizan para funciones con comportamientos cíclicos o periódicos y las bases monomiales se utilizan cuando existen tendencias simples que pueden ser ajustadas mediante líneas rectas, polinomios cuadráticos, polinomios de orden superior, etc. Las bases más usuales son aquellas basadas en splines debido a que son más flexibles y se ajustan de mejor manera a diferentes comportamientos. Dentro de estas se encuentran B-splines, M-splines, I-splines, y fucniones de potencia truncadas.
ELegir BASES
En general, el número de bases utilizado puede variar según varios factores, como la complejidad del problema, la cantidad de datos disponibles y las limitaciones computacionales. Aquí hay algunas pautas generales que puedes considerar al seleccionar el número de bases.
Complejidad del problema: Si el problema es complejo y se espera que requiera una representación más rica de las funciones, puede ser necesario aumentar el número de bases. Esto permitire capturar mejor las variaciones y las formas más complejas de las funciones que estás modelando.
Cantidad de datos: Si tienes una gran cantidad de datos disponibles, puedes permitirte utilizar un número mayor de bases. Sin embargo, ten en cuenta que un exceso de bases puede llevar a un sobreajuste (overfitting) si no se tiene cuidado.
Limitaciones computacionales: Si tienes restricciones de recursos computacionales, como memoria o tiempo de cálculo limitados, debes considerar limitar el número de bases para evitar problemas de rendimiento.
Bases monomiales
Las bases monomiales requieren el dominio y el número de base. Por ejemplo, la base monomial con K=7 funciones de base definidas en el intervalo [0,3] se puede construir con:
bbasis_obj = create.monomial.basis(rangeval=c(0,3), nbasis = 7)
Ftemp <- Data2fd(argvals=seq(0,3,length= 24 ),C, basisobj=bbasis_obj)
plot(Ftemp, main="Datos suavizados Velocidad_viento Monomiales (k=7)", xlab = "Horas", ylab = "Velocidad Viento")
## [1] "done"Bases de fourier
Para utilizarlas bases de Fourier se requiere que se defina el dominio, el periodo de oscilación y el número de funciones de base.
fbasis_obj <- create.fourier.basis(rangeval=c(0,3), nbasis=7, period = 12)
Ftemp <- Data2fd(seq(0,3,length= 24),C, basisobj=fbasis_obj)
plot(Ftemp, main="Datos suavizados de la Velocidad del viento Fourier(k=7)", xlab = "Horas", ylab = "Velocidad Viento")
## [1] "done"B-splines
Las bases B-spline requieren el dominio, el número de funciones de la base y el orden.
bsbasis_obj <- create.bspline.basis(rangeval=c(0,3),nbasis=7, norder=4)
plot(bsbasis_obj)
Ftemp <- Data2fd(seq(0,3,length= 24),C, basisobj=bsbasis_obj)
plot(Ftemp, main="Datos suavizados de la Velocidad del viento B-splines(4,k=7)", xlab = "Horas", ylab = "Velocidad Viento")
## [1] "done"Estimación parametros
loglam = seq(0,0.05,0.001)
nlam = length(loglam)
dfsave = rep(NA,nlam)
gcvsave = rep(NA,nlam)
for (ilam in 1:nlam) {
lambda = loglam[ilam]
fdParobj = fdPar(bsbasis_obj, Lfdobj=NULL, lambda= lambda)
smoothlist = smooth.basis(seq(0,3,length= 24), C, fdParobj)
dfsave[ilam] = smoothlist$df
gcvsave[ilam] = sum(smoothlist$gcv)
}
par(mfrow=c(1,1))
plot(loglam,gcvsave,xlab=expression(lambda),ylab=expression(GCV(lambda)),
main="Parametros de suavizamiento versus GCV",type="b",cex=0.7)
Tomando como referencia los 27 municipios en la base de datos de interés, se evidencia que el punto que minimiza la estadística GCV se encuentra alrededor de v=0.001 valores superiores antes de 0.001 y después de 0.004. Se destaca que dicho parámetro minimiza la suma del criterio GCV de todas las curvas que se encuentran dentro de la muestra, es decir, el valor de dicho parámetro resulta ser el óptimo en relación con el proceso de suavizamiento de todas las curvas dentro de la muestra.
Estadistica descriptiva
Ftemp es usado con bsplines
meanfdh <- mean.fd(Ftemp)
varfdh <-var.fd(Ftemp)
stdvfdh <- stddev.fd(Ftemp)
plot(Ftemp,col=8, lty=1, xlab = "Hora", ylab = "Velocidad del viento")## [1] "done"lines(meanfdh,col=2,lwd=2)
par(mfrow=c(1,2))
plot(varfdh, main=" Superficie de varianza", xlab = "v", ylab = "s")
plot(stdvfdh, main="Desviacion estándar", xlab = "Hora", ylab = "Velocidad del viento")
## [1] "done"ACP sin las curvas atípicas
CC <- C
CC <- CC[,-7]
CC <- CC[,-21]
par(mfrow=c(1,1))
horas <-0:23
plot(horas, CC[,1], type="b", ylim=c(0, 3), xlab="Horas", ylab="Velocidad Viento",col=1)
for(i in 2:ncol(CC)){
lines(horas, CC[,i], type="b", col=i)}
BSpl <- create.bspline.basis(rangeval=c(0,3),nbasis=7, norder=4)
FtempCC <- Data2fd(CC,argvals=seq( 0 , 3 , length= 24 ), basisobj=BSpl)
plot(FtempCC, main="Datos suavizados de la Velocidad del viento B-splines(4,k=8)", xlab = "Horas", ylab = "Velocidad Viento")
## [1] "done"PCA = pca.fd(FtempCC, 6)
par(mfrow=c(2,3))
plot(PCA)
(PCA$varprop[1]+PCA$varprop[2])*100## [1] 90.49514Tomaremos una varianza del 90%. Con dos funciones propias se explica el 90% de la varianza.
Gráfico de los scores
par(mfrow=c(1,1))
plot(x = PCA$scores[,1], y = PCA$scores[,2], cex=0.001,pch=16)
text(x = PCA$scores[,1], y = PCA$scores[,2] , seq(14,16,1), cex = 1)
Organizando la tabla de datos
CALDAS1 <- sqldf("select Municipio, Latitud, Longitud
from CALDAS
group by Municipio ")
CALDAS1 <- CALDAS1[-7,]
CALDAS1 <- CALDAS1[-21,]
CALDAS1## Municipio Latitud Longitud
## 1 AGUADAS 5.610248 -75.45487
## 2 ANSERMA 5.236476 -75.78434
## 3 ARANZAZU 5.271196 -75.49129
## 4 BELALCAZAR 4.993821 -75.81190
## 5 CHINCHINA 4.985813 -75.60748
## 6 FILADELFIA 5.297098 -75.56247
## 8 LA MERCED 5.396372 -75.54657
## 9 MANIZALES 5.057687 -75.49105
## 10 MANZANARES 5.255708 -75.15284
## 11 MARMATO 5.474002 -75.60025
## 12 MARQUETALIA 5.297525 -75.05309
## 13 MARULANDA 5.284335 -75.25940
## 14 NEIRA 5.166706 -75.51984
## 15 NORCASIA 5.574794 -74.88954
## 16 PACORA 5.527175 -75.45962
## 17 PALESTINA 5.017839 -75.62446
## 18 PENSILVANIA 5.383279 -75.16030
## 19 RIOSUCIO 5.423672 -75.70207
## 20 RISARALDA 5.164192 -75.76727
## 21 SALAMINA 5.403019 -75.48722
## 23 SAN JOSE 5.082311 -75.79206
## 24 SUPIA 5.446834 -75.64966
## 25 VICTORIA 5.317418 -74.91116
## 26 VILLAMARIA 5.043807 -75.51376
## 27 VITERBO 5.062663 -75.87061CALDAS1<-data.frame(CALDAS1,PCA$scores[,1])
CALDAS2<-data.frame(CALDAS1,PCA$scores[,2])
shp_map <- read_sf("C:/Users/Lenovo/Desktop/Caldas/CaldasM.shp")
shp_map1<-shp_map[5]
Cal <- shp_map1
ggplot() +
geom_sf(data=Cal)
CALDAS1<-data.frame(CALDAS1,PCA$scores[,1])
CALDAS2<-data.frame(CALDAS1,PCA$scores[,2])
names(CALDAS2)## [1] "Municipio" "Latitud" "Longitud"
## [4] "PCA.scores...1." "PCA.scores...1..1" "PCA.scores...2."names(CALDAS1)## [1] "Municipio" "Latitud" "Longitud"
## [4] "PCA.scores...1." "PCA.scores...1..1"colnames(CALDAS2)<-c("Municipio","Latitud","Longitud","ScoreS2")
colnames(CALDAS1)<-c("Municipio","Latitud","Longitud","ScoreS1")Mapa del departamento de Caldas con sus respectivos scores
library(rcartocolor)## Warning: package 'rcartocolor' was built under R version 4.2.3library(plotly)
map_variable<- function(datos, variable, map) {
plot1 <- ggplot() +
geom_sf(data = map, size = 0.3) +
geom_point(data = datos, aes_string(x = "Longitud", y = "Latitud",
color = variable)) +
scale_color_viridis_c(direction = -1) +
labs(
x = "",
y = "",
title = paste(variable, "Mapa Dept Caldas")
) +
theme_void() +
theme(
plot.title = element_text(size = 14,
face = "bold",
colour = "black"),
plot.margin = unit(c(1, 1, 1.5, 1.2), "cm"))
return(plot1)
}
pl1 <- map_variable(CALDAS1,
"ScoreS1",
shp_map1)## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation ideoms with `aes()`ggplotly(pl1)pl2 <- map_variable(CALDAS2,
"ScoreS2",
shp_map1)
ggplotly(pl2)Kriging Score 2
Curvas no suavizadas de la velocidad del viento por cada hora en los diferentes municipios del Departamento de Caldas
## Warning: package 'pander' was built under R version 4.2.3## Latitud Longitud Region Departamento
## Min. :4.986 Min. :-75.87 Length:2241 Length:2241
## 1st Qu.:5.082 1st Qu.:-75.65 Class :character Class :character
## Median :5.297 Median :-75.51 Mode :character Mode :character
## Mean :5.283 Mean :-75.44
## 3rd Qu.:5.424 3rd Qu.:-75.16
## Max. :5.610 Max. :-74.67
## Municipio Fecha Hora Velocidad_Viento
## Length:2241 Length:2241 Length:2241 Min. :0.0000
## Class :character Class :character Class :character 1st Qu.:0.6000
## Mode :character Mode :character Mode :character Median :1.0000
## Mean :0.9963
## 3rd Qu.:1.3000
## Max. :2.5000## [1] "AGUADAS" "ANSERMA" "ARANZAZU" "BELALCAZAR" "CHINCHINA"
## [6] "FILADELFIA" "LA DORADA" "LA MERCED" "MANIZALES" "MANZANARES"
## [11] "MARMATO" "MARQUETALIA" "MARULANDA" "NEIRA" "NORCASIA"
## [16] "PACORA" "PALESTINA" "PENSILVANIA" "RIOSUCIO" "RISARALDA"
## [21] "SALAMINA" "SAMANA" "SAN JOSE" "SUPIA" "VICTORIA"
## [26] "VILLAMARIA" "VITERBO"
Veamos un plot de los scores en las dos dimensiones que decidimos retener dado que explican un poco más del 90% de la varianza de los datos
Organizando los datos…
## Municipio Latitud Longitud
## 1 AGUADAS 5.610248 -75.45487
## 2 ANSERMA 5.236476 -75.78434
## 3 ARANZAZU 5.271196 -75.49129
## 4 BELALCAZAR 4.993821 -75.81190
## 5 CHINCHINA 4.985813 -75.60748
## 6 FILADELFIA 5.297098 -75.56247
## 8 LA MERCED 5.396372 -75.54657
## 9 MANIZALES 5.057687 -75.49105
## 10 MANZANARES 5.255708 -75.15284
## 11 MARMATO 5.474002 -75.60025
## 12 MARQUETALIA 5.297525 -75.05309
## 13 MARULANDA 5.284335 -75.25940
## 14 NEIRA 5.166706 -75.51984
## 15 NORCASIA 5.574794 -74.88954
## 16 PACORA 5.527175 -75.45962
## 17 PALESTINA 5.017839 -75.62446
## 18 PENSILVANIA 5.383279 -75.16030
## 19 RIOSUCIO 5.423672 -75.70207
## 20 RISARALDA 5.164192 -75.76727
## 21 SALAMINA 5.403019 -75.48722
## 22 SAMANA 5.413081 -74.99224
## 23 SAN JOSE 5.082311 -75.79206
## 24 SUPIA 5.446834 -75.64966
## 26 VILLAMARIA 5.043807 -75.51376
## 27 VITERBO 5.062663 -75.87061## [1] "Municipio" "Latitud" "Longitud" "PCA.scores...1."## Municipio Latitud Longitud Score2
## 1 AGUADAS 5.610248 -75.45487 -0.096558332
## 2 ANSERMA 5.236476 -75.78434 -0.694643228
## 3 ARANZAZU 5.271196 -75.49129 0.558259462
## 4 BELALCAZAR 4.993821 -75.81190 -0.965254767
## 5 CHINCHINA 4.985813 -75.60748 0.077234525
## 6 FILADELFIA 5.297098 -75.56247 0.146873338
## 8 LA MERCED 5.396372 -75.54657 0.008918755
## 9 MANIZALES 5.057687 -75.49105 0.890670289
## 10 MANZANARES 5.255708 -75.15284 0.326156017
## 11 MARMATO 5.474002 -75.60025 -0.335150220
## 12 MARQUETALIA 5.297525 -75.05309 0.889840807
## 13 MARULANDA 5.284335 -75.25940 -0.158013409
## 14 NEIRA 5.166706 -75.51984 0.643588755
## 15 NORCASIA 5.574794 -74.88954 0.375379514
## 16 PACORA 5.527175 -75.45962 -0.042957556
## 17 PALESTINA 5.017839 -75.62446 -0.043122312
## 18 PENSILVANIA 5.383279 -75.16030 0.338011421
## 19 RIOSUCIO 5.423672 -75.70207 -0.566655582
## 20 RISARALDA 5.164192 -75.76727 -0.805837213
## 21 SALAMINA 5.403019 -75.48722 0.204226519
## 22 SAMANA 5.413081 -74.99224 -0.959743191
## 23 SAN JOSE 5.082311 -75.79206 -0.431079762
## 24 SUPIA 5.446834 -75.64966 0.655434414
## 26 VILLAMARIA 5.043807 -75.51376 0.846367227
## 27 VITERBO 5.062663 -75.87061 -0.861945472## Longitud Latitud Score2
## 1 -75.45487 5.610248 -0.09655833
## 2 -75.78434 5.236476 -0.69464323
## 3 -75.49129 5.271196 0.55825946
## 4 -75.81190 4.993821 -0.96525477
## 5 -75.60748 4.985813 0.07723452
## 6 -75.56247 5.297098 0.14687334
## [1] "Municipio" "Latitud" "Longitud" "Score2" "Longitud.1"
## [6] "Latitud.1"## List of 2
## $ Longitud.1 , Latitud.1 : num [1:25, 1:2] 449624 413082 445561 409993 432655 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:25] "1" "2" "3" "4" ...
## .. ..$ : chr [1:2] "Longitud.1" "Latitud.1"
## $ data : num [1:25] -0.0966 -0.6946 0.5583 -0.9653 0.0772 ...
## - attr(*, "class")= chr "geodata"Explorando anisotropía
## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
Por lo que de acuerdo a los resultados obtenidos mostrados en los diferentes plots, se puede decir que la velocidad del viento en el departamento de Caldas es Anisotrópica pues la estructura de los semivariogramas varía conforme cambia la dirección.
Explorando dos posibles modelos de semovariograma
par(mfrow=c(2,3))
coordinates(CALDAS2) <- ~ Longitud.1 + Latitud.1
g_obj <- gstat(id = "Score2" ,
formula = Score2 ~ 1,
data = CALDAS2)
var_s <- variogram(g_obj)
var_s <- var_s[-1, ]
plot(var_s, pl = T)
show.vgms()
vgm_model1 <- vgm("Wav",
psill = 0.35,
range = 31000)
plot(var_s, vgm_model1, pl = T)
vgm_model2 <- vgm("Gau",
psill = 0.35,
range = 22000)
plot(var_s, vgm_model2, pl = T)
par(mfrow=c(1,2))
fitted_vgm1 <- fit.variogram(object = var_s,
model = vgm_model1)
plot(var_s, fitted_vgm1, pl = T)
fitted_vgm2 <- fit.variogram(object = var_s,
model = vgm_model2)## Warning in fit.variogram(object = var_s, model = vgm_model2): No convergence
## after 200 iterations: try different initial values?plot(var_s, fitted_vgm2, pl = T)
var_s$fitted1 <- variogramLine(fitted_vgm1, dist_vector = var_s$dist)$gamma
var_s$fitted2 <- variogramLine(fitted_vgm2, dist_vector = var_s$dist)$gamma
cmr <- function(obs, fit) {
return(mean((obs - fit)**2))
}
tabla <- data.frame(rbind(fitted_vgm1,
fitted_vgm2
),
CMR =
c(cmr(var_s$gamma, var_s$fitted1),
cmr(var_s$gamma, var_s$fitted2)
))
pander::pander(tabla[c(1:3,10)])| model | psill | range | CMR |
|---|---|---|---|
| Wav | 0.4184 | 31000 | 0.009602 |
| Gau | 0.2392 | 2391 | 0.02428 |
De acuerdo a lo presentado en la tabla resumen escogimos el modelo Wave con silla 0.4184 y rango 31000 para modelar la semivarianza, ya que tiene un CMR mas bajo.
Predicción de los datos funcionales
g_obj <- gstat(id = "Score2",
formula = Score2 ~ 1,
model = fitted_vgm1, # Modelo ajustado 1
data = CALDAS2)
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 fieldsclass(map)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"new <- sp::spsample(map, n = 20000, type = "regular")
CRS_UTM_NY = (c("+proj=utm +zone=18 +datum=NAD83
+units=m +no_defs"))
utm_coord = spTransform(deci_coord, CRS(CRS_UTM_NY))
plot(utm_coord, axes = TRUE, main = "UTM Coordinates", col = "green",cex = 0.7)
utm_coord_df = as.data.frame(utm_coord)
utm_coord2 = spTransform(new, CRS(CRS_UTM_NY))
utm_coord2<-as.data.frame(utm_coord2)
utm_coord2<-st_as_sf(utm_coord2,coords = c("x1","x2"))
predic <- predict(g_obj, newdata = utm_coord2)## [using ordinary kriging]Mapa de predicciones
par(mfrow=c(1,1))
plot1 <- ggplot() +
geom_sf(aes(col = Score2.pred), data = predic) +
scale_color_continuous(na.value = "#111111")
library(viridis)## Warning: package 'viridis' was built under R version 4.2.3plot1 <- ggplot() +
geom_sf(aes(col = Score2.pred), data = predic) +
scale_color_viridis_c(na.value = "#111111")
plot1 <- ggplot() +
geom_sf(aes(col = Score2.pred), data = predic) +
scale_color_continuous(na.value = "#111111")
plot1 <- ggplot() +
geom_sf(aes(col = Score2.pred), data = predic) +
scale_color_gradientn(colors = c("red", "orange", "blue"), na.value = "#111111")
plot1 
Mapa de varianza de las predicciones
plot12<-ggplot()+
geom_sf(aes(col=Score2.var),data=predic)+
scale_fill_viridis_c(option = "magma",
direction = -1) +
theme_void()
plot12 
Y asi, podemos observar los mapas de predicción de los scores y el respectivo mapas con su varianza asociada, con lo cual podemos decir que las predicciones son bastante buenas dado que e general las respectivas varianzas están entre 0.05 y 0.10.
Ahora,llevando estas predicciones de nuevo a los datos funcionales tenemos:
coef_scores=(PCA$harmonics$coefs[,2]%*% t(CALDAS2$Score2))
vec<-PCA$meanfd$coefs
M<-NULL
MM<-NULL
for (i in 1:25) {
M<-coef_scores[,i]+vec
MM<-cbind(MM,M)
}
par(mfrow=c(1,1))
v1=(fd(MM, PCA$harmonics$basis))
plot(v1)
## [1] "done"coef_scores2=PCA$harmonics$coefs[,2] %*%t(predic$Score2.pred)
N<-NULL
NN<-NULL
for (j in 1:ncol(coef_scores2)) {
N<-coef_scores2[,j]+vec
NN<-cbind(NN,N)
}
v2=(fd(NN, PCA$harmonics$basis))
plot(v2)
## [1] "done"Kriging Score 1
Organizando los datos…
## Municipio Latitud Longitud
## 1 AGUADAS 5.610248 -75.45487
## 2 ANSERMA 5.236476 -75.78434
## 3 ARANZAZU 5.271196 -75.49129
## 4 BELALCAZAR 4.993821 -75.81190
## 5 CHINCHINA 4.985813 -75.60748
## 6 FILADELFIA 5.297098 -75.56247
## 8 LA MERCED 5.396372 -75.54657
## 9 MANIZALES 5.057687 -75.49105
## 10 MANZANARES 5.255708 -75.15284
## 11 MARMATO 5.474002 -75.60025
## 12 MARQUETALIA 5.297525 -75.05309
## 13 MARULANDA 5.284335 -75.25940
## 14 NEIRA 5.166706 -75.51984
## 15 NORCASIA 5.574794 -74.88954
## 16 PACORA 5.527175 -75.45962
## 17 PALESTINA 5.017839 -75.62446
## 18 PENSILVANIA 5.383279 -75.16030
## 19 RIOSUCIO 5.423672 -75.70207
## 20 RISARALDA 5.164192 -75.76727
## 21 SALAMINA 5.403019 -75.48722
## 22 SAMANA 5.413081 -74.99224
## 23 SAN JOSE 5.082311 -75.79206
## 24 SUPIA 5.446834 -75.64966
## 26 VILLAMARIA 5.043807 -75.51376
## 27 VITERBO 5.062663 -75.87061## [1] "Municipio" "Latitud" "Longitud" "PCA.scores...1."## Longitud Latitud ScoreS1
## 1 -75.45487 5.610248 -0.09655833
## 2 -75.78434 5.236476 -0.69464323
## 3 -75.49129 5.271196 0.55825946
## 4 -75.81190 4.993821 -0.96525477
## 5 -75.60748 4.985813 0.07723452
## 6 -75.56247 5.297098 0.14687334
## [1] "Municipio" "Latitud" "Longitud" "ScoreS1" "Longitud.1"
## [6] "Latitud.1"## List of 2
## $ Longitud.1 , Latitud.1 : num [1:25, 1:2] 449624 413082 445561 409993 432655 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:25] "1" "2" "3" "4" ...
## .. ..$ : chr [1:2] "Longitud.1" "Latitud.1"
## $ data : num [1:25] -0.0966 -0.6946 0.5583 -0.9653 0.0772 ...
## - attr(*, "class")= chr "geodata"Explorando anisotropía
## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
Por lo que de acuerdo a los resultados obtenidos mostrados en los diferentes plots, se puede decir que la velocidad del viento en el departamento de Caldas es Anisotrópica pues la estructura de los semivariogramas varía conforme cambia la dirección.
Explorando dos posibles modelos de semovariograma
par(mfrow=c(2,3))
coordinates(CALDAS3) <- ~ Longitud.1 + Latitud.1
g_obj <- gstat(id = "ScoreS1" ,
formula = ScoreS1 ~ 1,
data = CALDAS3)
var_s <- variogram(g_obj)
var_s <- var_s[-1, ]
plot(var_s, pl = T)
show.vgms()
vgm_model1 <- vgm("Wav",
psill = 0.35,
range = 31000)
plot(var_s, vgm_model1, pl = T)
vgm_model2 <- vgm("Gau",
psill = 0.35,
range = 22000)
plot(var_s, vgm_model2, pl = T)
par(mfrow=c(1,2))
fitted_vgm1 <- fit.variogram(object = var_s,
model = vgm_model1)
plot(var_s, fitted_vgm1, pl = T)
fitted_vgm2 <- fit.variogram(object = var_s,
model = vgm_model2)## Warning in fit.variogram(object = var_s, model = vgm_model2): No convergence
## after 200 iterations: try different initial values?plot(var_s, fitted_vgm2, pl = T)
var_s$fitted1 <- variogramLine(fitted_vgm1, dist_vector = var_s$dist)$gamma
var_s$fitted2 <- variogramLine(fitted_vgm2, dist_vector = var_s$dist)$gamma
cmr <- function(obs, fit) {
return(mean((obs - fit)**2))
}
tabla <- data.frame(rbind(fitted_vgm1,
fitted_vgm2
),
CMR =
c(cmr(var_s$gamma, var_s$fitted1),
cmr(var_s$gamma, var_s$fitted2)
))
pander::pander(tabla[c(1:3,10)])| model | psill | range | CMR |
|---|---|---|---|
| Wav | 0.4184 | 31000 | 0.009602 |
| Gau | 0.2392 | 2391 | 0.02428 |
De acuerdo a lo presentado en la tabla resumen escogimos el modelo Wave con silla 0.4184 y rango 31000 para modelar la semivarianza, ya que tiene un CMR mas bajo.
Predicción de los datos funcionales
g_obj <- gstat(id = "ScoreS1",
formula = ScoreS1 ~ 1,
model = fitted_vgm1, # Modelo ajustado 1
data = CALDAS3)
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 fieldsclass(map)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"new <- sp::spsample(map, n = 20000, type = "regular")
CRS_UTM_NY = (c("+proj=utm +zone=18 +datum=NAD83
+units=m +no_defs"))
utm_coord = spTransform(deci_coord, CRS(CRS_UTM_NY))
plot(utm_coord, axes = TRUE, main = "UTM Coordinates", col = "green",cex = 0.7)
utm_coord_df = as.data.frame(utm_coord)
utm_coord2 = spTransform(new, CRS(CRS_UTM_NY))
utm_coord2<-as.data.frame(utm_coord2)
utm_coord2<-st_as_sf(utm_coord2,coords = c("x1","x2"))
predic <- predict(g_obj, newdata = utm_coord2)## [using ordinary kriging]Mapa de predicciones
Mapa de varianza de las predicciones
plot12<-ggplot()+
geom_sf(aes(col=ScoreS1.var),data=predic)+
scale_fill_viridis_c(option = "magma",
direction = -1) +
theme_void()
plot12 
Y asi, podemos observar los mapas de predicción de los scores y el respectivo mapas con su varianza asociada, con lo cual podemos decir que las predicciones son bastante buenas dado que e general las respectivas varianzas están entre 0.05 y 0.10.
Ahora,llevando estas predicciones de nuevo a los datos funcionales tenemos:
coef_scores=(PCA$harmonics$coefs[,1]%*% t(CALDAS3$ScoreS1))
vec<-PCA$meanfd$coefs
M<-NULL
MM<-NULL
for (i in 1:25) {
M<-coef_scores[,i]+vec
MM<-cbind(MM,M)
}
par(mfrow=c(1,1))
v1=(fd(MM, PCA$harmonics$basis))
plot(v1)
## [1] "done"coef_scores1=PCA$harmonics$coefs[,1] %*%t(predic$ScoreS1.pred)
N<-NULL
NN<-NULL
for (j in 1:ncol(coef_scores1)) {
N<-coef_scores1[,j]+vec
NN<-cbind(NN,N)
}
v2=(fd(NN, PCA$harmonics$basis))
plot(v2)
## [1] "done"Conclusiones
La semivarianza aumenta con la distancia hasta alcanzar un valor máximo (silla), a partir del cual la semivarianza tiende a estabilizarse y no aumentar significativamente con la distancia.
Un buen estimador empírico del semivariagroma para este caso es un powered exponential con silla 0.31, rango 327775.39, nugget (pepita) 0.04 y kappa 0.5. Y los estimadores teóricos óptimos fue el de mínimos cuadrados ordinarios con con silla 0.3726, rango 327775.3900, nugget 0.04 y kappa 0.5. Y el max-logverosimil con silla 3.434e-01, rango 3.278e+05, nugget 0.04 y kappa 0.5.Es preciso anotar que se trata de un modelo de semivariograma con estacionariedad de segundo orden.
El modelo powered.exponential no debería ser utilizado en la predicción espacial debido a las inestabilidades numéricas que produce en los algoritmos kriging (especialmente cuando el efecto nugget es grande), sin embargo para este caso se consideró dado que el efecto nugget es muy pequeño (0.04).
Se puede observar además que el Kriging simple realizado para la predicción de la velocidad del viento en un punto no muestreado de la región Andina Colombiana como lo es Palmira Valle del Cauca, el valor obtenido fue de 1.7/s para el día 12 de febrero a las 6pm y cuya varianza del error de predicción obtenido fue bastante baja (0.0068).
Sin embargo también se llevó a cabo Kriging simple para 100000 puntos no muestreados de la región Andina cuyos valores mostrados en el mapa se pueden agrupar 5 intervalos,donde el valor máximo es de 3.95m/s (zonas costeras), valor mínimo es de 0.3092m/s la cual corresponde a la velocidad del viento en zonas más bajas de la región andina. En general se tiene registrado que en las zonas montañosas la velocidad del aire puede ser mayor y persistente debido a la topografía accidentada.Esto justifica la distribución de los colores en el mapa de predicción.
En cuanto al Kriging espacio tiempo se puede decir es más complejo en cuanto a su aplicación por la demanda computacional que tiene debido a la complejidad del manejo de matrices de grandes y vectores de grandes dimensiones de acuerdo a la cantidad de datos que se tienen disponibles y a los rezagos espaciales.Es por eso que en esta área se recomienda el uso de datos funcionales.
El kriging funcional es una herramienta de predicción de la Geoestadística funcional que aprovecha las dependencias espacio-temporales existentes en los registros medioambientales para realizar predicciones en localizaciones e instantes temporales no observados. En el caso del modelamiento de la velocidad del viento en el departamento de Caldas, el uso de datos funcionales y kriging funcional permitirió obtener predicciones más precisas en todas las ubicaciones y en cada hora del día de una manera computacionalmente más eficiente.