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 1000 metros sobre el nivel del mar hasta más de 5000 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)
## Warning: package 'readxl' was built under R version 4.2.3A 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 resultado 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.7 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.85290 0.2429200
## 2 -76.91534 0.2585304
## 3 -76.89973 0.2585304
## 4 -76.88412 0.2585304
## 5 -76.86851 0.2585304
## 6 -76.85290 0.2585304library(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 293784.2 26864.10
## 2 286832.5 28591.46
## 3 288570.5 28591.20
## 4 290308.5 28590.94
## 5 292046.5 28590.69
## 6 293784.4 28590.43Semivariograma 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")Validación cruzada
La validación cruzada en gstat está implementada en la función krige.cv() por lo que se emplea esa función para realizarla. Veamos:
Z<-Velocidad_Viento1
x<-Este1
y<-Norte1
library(gstat)
library(sf)
colnames(Randina_utm)<-c('x','y','Z')
andina_sf <- st_as_sf(Randina_utm, coords = c("x", "y"), remove = FALSE, agr = "constant")
ve.fit1exp = vgm(psill = 0.3726, model = "Exp", range = 327775.3900,
kappa = 0.5)
valcru <-krige.cv(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),locations=andina_sf, model=ve.fit1exp,beta=beta)Realizada la validación es posible calcular estadísticos resumen como el error medio (ME), error cuadrático medio (MSE), media del cociente de la desviación cuadrática (mean squared deviation ratio, MSDR), raíz del error cuadrático medio (RMSE), la RMSE relativa a la media de los observados (RMSE_rel)
ME <- mean(valcru$residual)
MSE <- mean(valcru$residual ^ 2)
MSDR <- mean(valcru$zscore ^ 2)
RMSE <- sqrt(mean(valcru$residual ^ 2))
RMSE_rel <-
sqrt(mean(valcru$residual ^ 2)) /
mean(valcru$observed) * 100
r <- cor(valcru$observed,
valcru$observed - valcru$residual)
tabla <- data.frame(ME, MSE, RMSE,
RMSE_rel, MSDR, r)
tabla## ME MSE RMSE RMSE_rel MSDR r
## 1 -0.003908157 0.02214213 0.1488023 11.03392 1.307571 0.9682361El error promedio es cercano a cero, lo que indica que el modelo tiende a hacer predicciones cercanas al valor real en promedio.
Error cuadrático medio (MSE) es relativamente bajo, lo que sugiere que el modelo tiene un buen rendimiento en términos de ajuste a los datos.
Raíz del error cuadrático medio (RMSE) también es bajo, lo que indica que el modelo tiene un error de predicción aceptable y las predicciones están cerca de los valores reales en general.
El RMSE relativo indica que el modelo tiene un error relativo del 11.03% en comparación con el rango de los valores de los datos originales. Esto implica que el modelo puede tener un rendimiento aceptable, pero el error relativo debe considerarse en el contexto del dominio específico del problema.
Desviación estándar media relativa (MSDR) es aproximadamente 1.3076, lo que indica que los errores del modelo tienen una dispersión relativa moderada en comparación con el rango de los valores de los datos originales.
El valor del coeficiente de correlación es aproximadamente 0.9682, lo que sugiere una correlación positiva fuerte entre las predicciones y los valores reales. Esto indica que el modelo tiene la capacidad de capturar la tendencia general de los datos y realizar predicciones coherentes.
En general, los resultados indican que el modelo tiene un buen rendimiento en la tarea de predicción y está en sintonía con los valores reales.
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.
De acuerdo a la validación cruzada realizada se puede decir que el modelo de semivariograma teórico expuesto permite realizar un kriging con buen rendimiento en la tarea de predicción y está en sintonía con los valores reales.