EVALUACIÓN DE MÉTODOS DE INTERPOLACIÓN PARA GENERAR MAPAS DE PRECIPITACIÓN EN EL DEPARTAMENTO DE ANTIOQUIA, COLOMBIA.
Zuly Olaya Acosta, Ing. topografica. Estudiante maestria en geomática. Universidad Nacional de Colombia
Nadia Luque Sanabria, Ing. Agronoma. Estudiante maestria en geomática. Universidad Nacional de Colombia
30.06.2022
1. Introducción
La producción de información y conocimientos de los fenómenos atmosféricos tales como la precipitación es parte fundamental en la implementación de la política económica, social y ambiental en el país, además de ser un valioso recurso para el sector agrícola, que se basa en estos productos como fuente para establecer tiempos y lugares adecuados para la siembra (Sánchez, 2022).Siendo la precipitación un fenómeno atmosférico importante para varios actores, se resalta la necesidad de elaborar modelos que se ajusten adecuadamente a la realidad de la zona y así poder realizar pronósticos de calidad y tomar decisiones acertadas. Es así como surge la necesidad de evaluar diversos métodos de interpolación y encontrar el que brinde la mejor exactitud temática en el departamento de Antioquia, que se caracteriza por tener precipitaciones que van desde los 1.000 mm a los 6.500 mm.
Los datos de referencia o la información específica de cada zona son vitales para realizar un análisis espacial. Sin embargo, es complejo tener información de todas las áreas de interés. A veces solo se cuenta con información relativamente cercana, por lo que se deben buscar métodos para usar la única información disponible. Los métodos de interpolación se usan para estimar valores en ubicaciones donde no tenemos información y es desconocida, basándose en valores cercanos que se muestran espacialmente. Estos métodos asumen una alta correlación entre la información o los puntos más cercanos que los que se encuentran más lejos (Lixin, Weitian & Reinhard, 2020). Las técnicas de interpolación espacial son esenciales para estimar las variables biofísicas de las ubicaciones no muestreadas. Estas se han desarrollado y aplicado en diversas disciplinas. Factores como el tamaño de la muestra, el diseño de la muestra y las propiedades de los datos afectan las estimaciones de los métodos de interpolación existentes, por eso es tan difícil seleccionar el método de interpolación espacial apropiado para un conjunto de datos especifico. Cada método de interpolación tiene sus supuestos y características específicas, que se pueden discutir como globales versus locales, exactas versus inexactas, determinista versus estocástica y gradual versus abrupta (Li & Heap, 2008). La variedad de métodos disponibles trae consigo dudas acerca de cuál modelo es el más apropiado y eficiente, en cuanto a exactitud temática, de acuerdo con el tipo de evento a modelar, sobre todo porque los estudios previos no convergen hacia ningún método (Zimmerman et. al., 1999).
1.1 Objetivo
El objetivo de este trabajo fue evaluar el método de interpolación determinístico (IDW) y dos métodos de interpolación geoestadísticos (Kriging Simple y Ordinario), para establecer cuál de estos presenta una mayor exactitud temática en los datos de precipitación disponibles para el departamento de Antioquia, Colombia.
1.2 Conceptos teóricos
De los métodos de interpolación existentes se seleccionaron tres, Distancia Inversa Ponderada IDW, Kriging Ordinario y Kriging Simple, debido a que estos métodos son los que se usan con mayor frecuencia en el modelamiento de variables espaciales. El método kriging cuantifica la estructura espacial de los datos mediante el uso de variogramas los cuales predice mediante la interpolación, usando estadística. Se asume que los datos más cercanos a un punto conocido tienen mayor peso o influencia sobre la interpolación, influencia que va disminuyendo conforme se aleja del punto de interés. La medición de la probabilidad, efectuada por los métodos kriging, hace la diferencia con respecto a los métodos determinísticos para interpolaciones espaciales, de los cuales los más usados son el de ponderación de distancias inversas (IDW: Inverse Distance Weighting) (Burrough y McDonnell 1998). El método IDW es similar al kriging ordinario, ya que da más peso a los valores cercanos a un punto, pero posee una menor complejidad del cálculo. El IDW utiliza un algoritmo simple basado en distancias (Lu & Wong, 2008) Ambos modelos, kriging ordinario e IDW, asumen que las predicciones son una combinación lineal de los datos, como lo muestra la ecuación (Abedini & Nasseri, 2009, Li et al, 2018).
γ^’(S_0 ) = ∑_(i=1)^n▒〖λ_i γ(S_i ) 〗 (1)
Donde γ’es el valor estimado en el punto interpolado S0; es el peso dado al valor observado λ(si) en las cercanías del valor S0 Abedini & Nasseri, 2009. Este último parámetro hace la diferencia entre kriging y el IDW.
Kriging Simple
Para abordar el problema de kriging simple se debe suponer que hay una variable regionalizada estacionaria con media m y covarianza conocidas, situación que en casos prácticos se da con poca frecuencia. Entonces para tratar este caso se asume un modelo establecido en el cual el valor de la variable regionalizada será igual a la media más un error aleatorio con media cero. La diferencia con respecto a kriging ordinario es que los errores no son independientes, serán ellos justamente sobre los cuales se establecerá la estructura de correlación (Kleijnen, 2017). Dicho modelo suele expresarse de la siguiente forma:
Z(s)= m(s)+ ε(s) (2)Donde se supone conocido el variograma o el covariograma de E [Z(s)]= m(s) (3) Esto no implica que la media sea constante.
Kriging Ordinario
El kriging ordinario establece una condición adicional al sistema de ecuaciones del kriging simple para filtrar el valor desconocido de la media (Lichtenstern, 2013). Propone que el valor de la variable en el sitio no muestreado puede predecirse como una combinación lineal de las n variables aleatorias, como se muestra a continuación:
Z*(S_0 ) ∑_(i=1)^n▒λ_i Z(S_i ) (4)en donde los λi representan los pesos o ponderaciones de los valores de las variables en los sitios muestreados. Dichos pesos se calculan en función de la distancia entre los puntos muestreados y el punto donde se va a llevar a cabo la correspondiente predicción. La suma de los pesos debe ser igual a uno para que el valor esperado del predictor sea igual al valor esperado de la variable. Esto último se conoce como requisito de insesgamiento (Gómez M, 2005). Se dice que Z*(S0) es el mejor predictor lineal en este caso, porque los pesos se obtienen de tal forma que minimicen la varianza del error de predicción sujeto a que se cumpla el requisito de insesgamiento, es decir, que se minimice la expresión:
V(Z*(S_0 )- Z (S_0 ) ) sujeto a ∑_(i=1)^n▒λ_i =1 (5)La aplicación del método de los multiplicadores de Lagrange como técnica de optimización, en conjunto con la determinación de la matriz de covarianzas a partir de la estructura de autocorrelación espacial, permite la determinación de los pesos λi óptimos (Lichtenstern, 2013).
IDW
En el método IDW, las ubicaciones usadas tienen valores identificados y otras ubicaciones no identificadas tienen valores calculados, así se identifica una vecindad para el punto interpolado y se toma un promedio ponderado dentro de esta vecindad. Este método se usa para pronosticar no identificados para cualquier dato de ubicación geográfica, por ejemplo, precipitación, altura, profundidad, concentración de parámetros químicos, niveles de contaminación, etc. En este caso el usuario define la forma matemática de la función de ponderación y el tamaño del vecindario (expresado como radio o número de puntos) (Paramasivam & Venkatramanan, 2019). El método de interpolación IDW es determinístico, simple e intuitivo. Se basa en el principio de que los valores de las muestras más cercanas a la ubicación de predicción tienen más influencia en el valor de predicción de valores de muestra más separados. La principal desventaja de IDW es el efecto de “ojo de buey” alrededor de los puntos de muestreo (Achilleos, 2011).
2. Datos
2.1 Zona de estudio
El área de estudio corresponde al departamento de Antioquia, Colombia, localizado al noroeste del país sobre la cordillera central de los Andes, entre los 05º 26’ 20” y 08º 52’ 23” de latitud norte, y los 73º 53’ 11” y 77º 07’ 16” de longitud oeste. Cuenta con una superficie de 63.612 km² lo que representa el 5.6 % del territorio nacional. Limita por el Norte con el mar Caribe y los departamentos de Córdoba y Bolívar; por el Este con Bolívar, Santander y Boyacá; por el Sur con Caldas y Risaralda y por el Oeste con el departamento del Chocó, como se muestra en la figura 1.
library(tidyverse)
library(dplyr)
library(ggplot2)
library(sf)
library(RColorBrewer)
library(leaflet)
library(readxl)
library(rgdal)
library(gstat)
library (sf)
library(car)
library(aplpack)
library(scatterplot3d)
library(akima)
library(rgl)
library(lava)
library(dismo)Ant = st_read("D:/Dropbox/Asig_maestria_Geomatica/Geoestadistica/Proyecto_final/Municip_Antio.shp")## Reading layer `Municip_Antio' from data source
## `D:\Dropbox\Asig_maestria_Geomatica\Geoestadistica\Proyecto_final\Municip_Antio.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 125 features and 62 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -8585831 ymin: 604092.3 xmax: -8224427 ymax: 991819.5
## Projected CRS: WGS 84 / Pseudo-Mercatorggplot(Ant)+geom_sf()
Figura 1. Municipios de Antioquia
Antioquia se caracteriza por un relieve variado, representado por áreas planas localizadas en el valle del Magdalena y las zonas próximas al Chocó y el Urabá y una extensa área montañosa que hace parte de las cordilleras Central y Occidental, en donde se resaltan 202 altos importantes, con alturas que oscilan entre los 1.000 y los 4.080 metros sobre el nivel del mar. (Antioquia, 2014).
Ant = st_read("D:/Dropbox/Asig_maestria_Geomatica/Geoestadistica/Proyecto_final/ProyectoFinal_ZOlayaA_NLuqueS/Antioquia.shp")
names(Ant)
Ant.gcs <- st_transform(Ant, "+proj=longlat +datum=WGS84")
st_crs(Ant.gcs)
m<-leaflet() %>%
addTiles() %>%
addPolygons(data = Ant.gcs ,stroke = TRUE, fillOpacity = 0.8,
highlight = highlightOptions(weight = 3,
color = "red",
bringToFront = TRUE))mFigura 2. Ubicación del departamento de Antioquia en Colombia.
2.2 Descripción de datos
La evaluación de los métodos de interpolación se realizó con los datos de precipitación de las 169 estaciones meteorológicas que cubren el departamento de Antioquia pertenecientes a la red hidrológica y meteorológica de Colombia. Esta información fue obtenida de las bases de datos de la red pluviométrica del IDEAM y contiene: el valor promedio mensual de precipitación, código, tipo, nombre, municipio de ubicación, elevación, latitud y longitud. La cartografía base utilizada corresponde al “Límite Departamental” en formato shapefile obtenido del Sistema de Información Geográfica para la Planeación y el Ordenamiento Territorial SIG-OT disponible en http:/ /sigotn. igac. gov. co/ sigotn /default. aspx.
3. Métodos
La generación de la cartografía temática se hizo empleando los tres métodos mencionados anteriormente. Para el IDW se tomaron los datos y se ejecutó el proceso seleccionando la potencia y el número de vecinos mínimo y máximo. Por otro lado, la interpolación mediante métodos geoestadísticos involucró varios procesos previos como el Análisis Exploratorio de Datos Espaciales -AEDE- para la validación del supuesto de normalidad y el Análisis estructural y de tendencia, que se hizo con un nivel de confianza del 95% (Ver Figura 2). Luego se estimaron los parámetros del variograma: pepita, meseta parcial y σ2, y se calculó el variograma en diferentes ángulos con el fin de identificar anisotropía. Finalmente se evaluó el modelo mediante validación cruzada, que consistió en excluir un punto observado y recalcularlo mediante la interpolación. La diferencia entre el valor estimado y el valor real de todos los puntos será un indicador de ajuste del modelo. A continuación, se presenta la descripción de la metodología para la obtención de mapas de predicción.
A caption
4. Resultados
library(gstat)
library(sp)
library(raster)Antioquia<-shapefile("D:/Dropbox/Asig_maestria_Geomatica/Geoestadistica/Proyecto_final/Antio_project.shp")
Preci_Ant <- shapefile('D:/Dropbox/Asig_maestria_Geomatica/Geoestadistica/Proyecto_final/Antio_precip.shp')
#Se leen los datos en formato tabla
library(readxl)
setwd('D:/Dropbox/Asig_maestria_Geomatica/Geoestadistica/Proyecto_final/ProyectoFinal_ZOlayaA_NLuqueS')
Antio <- read_excel("Antio.xlsx", col_types = c("numeric","numeric", "numeric", "numeric", "numeric"))
View(Antio)
data = Antioproj4string(Preci_Ant) <- CRS("+init=epsg:3116")
proj4string(Antioquia) <- CRS("+init=epsg:3116")
Preci_Ant@bbox <-Antioquia@bbox4.1 Distancia Inversa Ponderada - Inverse Distance Weighted Method (IDW)
La distancia inversa ponderada es un método sencillo que realiza una interpolación en función de la distancia inversa elevada a una potencia p, entre más alto sea el valor de la potencia, los valores más cercanos tendrán mayor peso, la Figura 3. muestra el mapa generado mediante el método IDW.
# # Crear una cuadrícula vacía donde n es el número total de celdas
grd <- as.data.frame(spsample(Preci_Ant, "regular", n=319000))
# # Necesita averiguar cuál es el tamaño esperado del grd de salida
names(grd) <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd) <- TRUE # Crear objeto SpatialPixel
fullgrid(grd) <- TRUE # Crear objeto SpatialGrid
# Agregar información de proyección de Precipitacion a la cuadrícula vacía
proj4string(grd) <- proj4string(Preci_Ant)
# Interpolar las celdas de la cuadrícula utilizando un valor de potencia de 2 (idp = 2.0)
P.idw <- gstat::idw(Precip_Anu ~ 1, Preci_Ant, newdata=grd, idp=2.0)r <- raster(P.idw)
rr.m <- raster::mask(r, Antioquia)
r.mlibrary(tmap)## Warning: package 'tmap' was built under R version 4.1.3tm_shape(r.m) +
tm_raster(n=10,palette = "RdBu", auto.palette.mapping = FALSE,
title="Distancia inversa ponderada \nPrecipitación prevista \n(en mm)") +
tm_shape(Preci_Ant) + tm_dots(size=0.2) +
tm_legend(legend.outside=TRUE)## Warning: The argument auto.palette.mapping is deprecated. Please use midpoint
## for numeric data and stretch.palette for categorical data to control the palette
## mapping.
Un inconveniente con este método es que generó islas donde hay puntos
con valores extremos o donde hay muy pocos datos, como se puede ver en
la zona occidental del departamento. Además de mostrar variaciones
abruptas en la zona central. Además de que los valores de precipitación
no concuerdan con la interpolación (analizado visualmente) mostrando así
que es necesario generar otro método que ajuste mejor los valores de
precipitación.
P=Preci_Ant
IDW.out <- vector(length = length(P))
for (i in 1:length(P)) {IDW.out[i] <- gstat::idw(Precip_Anu ~ 1, P[-i,], P[i,], idp=2)$var1.pred}OP <- par(pty="s", mar=c(4,3,0,0))
plot(IDW.out ~ P$Precip_Anu, asp=1, xlab="Observado", ylab="Previsto", pch=16,
col=rgb(0,0,0,0.5))
abline(lm(IDW.out ~ P$Precip_Anu), col="red", lw=2,lty=2)
abline(0,1)
El error cuadrático medio (RMSE) para la interpolcación con este método
es de 726 mm aprox.
par(OP)
# Calcular RMSE
RMSE_IDW = sqrt( sum((IDW.out - P$Precip_Anu)^2) / length(P))
RMSE_IDW4.2 Poligonos de voronoi
Para corroborar la información anterior sobre los resultados del método de IDW con la información de precipitación, el análisis siguiente es generar los polígonos de Voronoi, el cual consiste en un método matemático que se basa en la proximidad que existe en una nube de puntos para regionalizar según el área de influencia del punto. El diagrama de Voronoi obtenido para el departamento de Antioquia (figura, 4), permite comprender que la mayoría de estaciones pluviométricas se encuentran ubicadas en la zona sur y central, ya que se representan polígonos más pequeños. En la zona limite se ubican los polígonos más grandes ya que existen pocas estaciones, dichas áreas deben ser cubiertas por una sola estación. Esto también indica que las zonas que contienen más información tienen a su vez menor error.
library(dismo)
xy=data.frame(data[,2:3])
Voro = voronoi(xy)## Loading required namespace: deldirplot(Voro, col = sf.colors(8))
## 4.3 Análisis exploratorio de datos El análisis exploratorio de datos
mostró una media y mediana diferentes por lo que se presumió que los
datos no se distribuyen normalmente, a esta conclusión se llegó
revisando el coeficiente de asimetría que indica que los datos se
acumulan en los valores bajos. (Ver Tabla 1.)
#Para llamar el script que calcula las estadisticas descriptivas
# Directorio de trabajo
setwd('D:/Dropbox/Asig_maestria_Geomatica/Geoestadistica/Proyecto_final/ProyectoFinal_ZOlayaA_NLuqueS')
source('programasR.txt')
tabla1=tabla.estadisticas.cont(Preci_Ant$Precip_Anu)
tabla1## x
## Mínimo 1041.70
## Máximo 6455.50
## Promedio 2864.55
## Mediana 2688.60
## Desviacion estándar 1005.69
## Desviacion mediana 515.60
## Asimetría 0.99
## Kurtosis 3.69
## Coef. Var. Promedio(%) 35.11
## Coef. Var. Mediana(%) 19.18Tabla 1. Estadística descriptiva de la variable precipitación
Un posterior análisis mediante el histograma y qq-plot evidenció que los valores de precipitación eran claramente asimétricos a la derecha como se ve en la figura 5. También se verificó la normalidad mediante el test de Shapiro Wilk, que arrojó un p-valor de 1.669e-07, indicando que los datos no eran normales por lo que se tuvieron que transformar mediante una función logarítmica.
tabla.estadisticas.cont(Preci_Ant[4:5])## Precip_Anu ELEVACION coords.x1 coords.x2 optional
## Mínimo 1041.70 4.00 -76.86 5.53 1
## Máximo 6455.50 2830.00 -73.94 8.85 1
## Promedio 2864.55 1130.39 -75.54 6.77 1
## Mediana 2688.60 1200.00 -75.52 6.63 1
## Desviacion estándar 1005.69 877.60 0.62 0.75 0
## Desviacion mediana 515.60 850.00 0.40 0.49 0
## Asimetría 0.99 0.09 -0.19 0.55 NaN
## Kurtosis 3.69 1.59 2.55 2.52 NaN
## Coef. Var. Promedio(%) 35.11 77.64 -0.81 11.12 0
## Coef. Var. Mediana(%) 19.18 70.83 -0.53 7.36 0op=par(mfrow=c(1,3))
hist(Preci_Ant$Precip_Anu, freq=F, main="HISTOGRAMA",xlab="Precipitación")
lines(density(Preci_Ant$Precip_Anu))
boxplot(Preci_Ant$Precip_Anu, main="BOX - PLOT")
plot(ecdf(Preci_Ant$Precip_Anu), main="E.C.D.F")
plot(Preci_Ant[,c(2,4)], col="red", xlab="Coordenada X", ylab="Precipitación")
plot(Preci_Ant[,c(3,4)], col="blue", xlab="Coordenada Y", ylab="Precipitación")
points(x=Preci_Ant$LONGITUD,y=Preci_Ant$LATITUD, col="red")
plot((Preci_Ant$LONGITUD)*(Preci_Ant$LATITUD), col="chartreuse4", Preci_Ant$Precip_Anu, xlab="CoordenadaX*CoordenadaY", ylab="Altura de cabeza")
datos1=data[2:5]
coordinates(datos1)=c("LONGITUD","LATITUD")
spplot(datos1, col.regions=rainbow(3))
hist(Preci_Ant$Precip_Anu,nclass=20, xlab="Precipitación", ylab="Frecuencia", main="")
qqnorm(Preci_Ant$Precip_Anu, xlab="", ylab="")
Figura 5. a) Histograma asimétrico y b) qq-plot de datos pluviométricos
sin distribución normal. Fuente: Propia
#Se realiza el test de Shapiro wilk para evaluar normalidad
Shapiro = shapiro.test(Preci_Ant$Precip_Anu)
Shapiro
#Como los datos no son normales se ejecuta la función que indica cúal es la transformación que se debe hacer
powerTransform(Preci_Ant$Precip_Anu)
#El resultado indica que se debe transformar con una función Logaritmica
data_pt=log(Preci_Ant$Precip_Anu)
#Se genera el histograma y el qq plot para verificar normalidad
hist(data_pt,nclass=20, xlab="Precipitación", ylab="Frecuencia", main="")
qqnorm(data_pt, xlab="", ylab="")
#Se ejecuta nuevamente el test de shapiro wilk
ShapiroLog = shapiro.test(data_pt)
ShapiroLogFigura 6. c) Histograma ajustado a una distribución normal con intervalo de confianza del 95% y d) QQ-plot de datos pluviométricos con distribución normal. Fuente: Propia
Además de evaluar la normalidad, se hizo un análisis de tendencia mediante una breve interpolación bilineal en la que no se pudo establecer relación alguna entre la latitud y la longitud (figura 7). También se implementó un modelo lineal que resultó significativo con un p-valor de 0.00006247 pero que explica la variable en tan solo un 9% según el R² ajustado, por lo que se consideró que la precipitación en el departamento de Antioquia no presentaba tendencia respecto a latitud y longitud.
#se convierten los datos a formato geodata para evaluar tendencia
geodata=st_as_sf(data, coords=c(2:3))
#se crea una superficie con un método de interpolación bivariado
x=data$LONGITUD
y=data$LATITUD
z=data_pt
inter=interp(x, y, z, xo=seq(min(x), max(x), length = 1000),yo=seq(min(y), max(y), length = 1000))
#se genera el gráfico
color=c("gray1", "gray2", "gray3", "gray4", "gray5", "gray6", "gray7", "gray8", "gray9", "gray10", "gray11", "gray12", "gray13", "gray14", "gray15", "gray16", "gray17", "gray18", "gray19", "gray20", "gray21", "gray22", "gray23", "gray24", "gray25", "gray26", "gray27", "gray28", "gray29", "gray30", "gray31", "gray32", "gray33", "gray34", "gray35", "gray36", "gray37", "gray38", "gray39", "gray40","gray41","gray42","gray43","gray44","gray45","gray46","gray47","gray48","gray49","gray50","gray51","gray52","gray53","gray54","gray55","gray56","gray57","gray58","gray59","gray60","gray61","gray62","gray63","gray64","gray65","gray66","gray67","gray68","gray69","gray70","gray71","gray72", "gray73", "gray74","gray75","gray76","gray77","gray78","gray79","gray80","gray81","gray82","gray83","gray84","gray85","gray86","gray87","gray88","gray89","gray90","gray91","gray92","gray93","gray94","gray95","gray96","gray97","gray98","gray99","gray100")
image(inter, col =color)4.4 Análisis Estructural
En este punto se realizó un modelo lineal para identificar si los datos presentaban tendencia respecto a la latitud y longitud encontrando que ambas variables resultaban ser significativas dentro del modelo, sin embargo, al revisar el coeficiente de determinación R2 se evidenció que el modelo sólo explica los datos en un 9% por lo que se consideró irrelevante y no se tuvo en cuenta.
Tabla 2. Modelo lineal
#se crea un modelo lineal para evaluar tendencia
lm=lm(data_pt~data$LATITUD+data$LONGITUD)
#se genera un resumen
summary(lm)Luego se realizó la selección y ajuste del semivariograma (figura 8), encontrando que en el punto once se mantiene la asíntota del variograma (sill), hasta el valor de 1.0286 podemos verificar si hay una correlación. El valor del rango es de 1.0286 y el nugget es de 0.0207.
#formato requerido por la librería "gstat" para aplicarle sus funciones
coordinates(data)=~LATITUD+LONGITUD
###Seleccion y ajuste de semivariograma###
Vario.Var=variogram(log(data$Precip_Anual)~LATITUD+LONGITUD,data)
plot(Vario.Var)
# En el punto 9 se observa la asintota del variograma (sill).
#El valor de 9 se sustituye en los parámetros de sill y range para obtener los parámetros
sill=Vario.Var[11,2]
range=Vario.Var[11,2]
nugget=min(Vario.Var[,3])Figura 8. Semivariograma ajustado. Fuente: Propia
Se calcularon los errores de los diferentes modelos de variograma, para establecer cual representaba mejor los datos de precipitación de Antioquia (Tabla 3).
Tabla 3. Calculo error para cada modelo de variograma
#Se generan 7 modelos permisibles con los parámetros de inicio
Esférico=vgm(sill,"Sph",range,nugget,variogramModel=Vario.Var)
Exponencial=vgm(sill,"Exp",range,nugget,variogramModel=Vario.Var)
Gaussiano=vgm(sill,"Gau",range,nugget,variogramModel=Vario.Var)
Lineal=vgm(sill,"Lin",range,nugget,variogramModel=Vario.Var)
Matern=vgm(sill,"Mat",range,nugget,variogramModel=Vario.Var)
Bessel=vgm(sill,"Bes",range,nugget,variogramModel=Vario.Var)
Pentaesf=vgm(sill,"Pen",range,nugget,variogramModel=Vario.Var)
#Se ajustan los modelos permisibles
fit.Esférico=fit.variogram(Vario.Var,Esférico)
fit.Exponencial=fit.variogram(Vario.Var,Exponencial)
fit.Gaussiano=fit.variogram(Vario.Var,Gaussiano)
fit.Lineal=fit.variogram(Vario.Var,Lineal)
fit.Matern=fit.variogram(Vario.Var,Matern)
fit.Bessel=fit.variogram(Vario.Var,Bessel)
fit.Pentaesf=fit.variogram(Vario.Var,Pentaesf)
#Se evalúan los modelos y se selecciona el que tenga menor error, en este caso, fue Exponencial "
ErEs = attr(fit.Esférico,"SSErr")
ErEx=attr(fit.Exponencial,"SSErr")
ErGau=attr(fit.Gaussiano,"SSErr")
ErLin=attr(fit.Lineal,"SSErr")
ErMat=attr(fit.Matern,"SSErr")
ErBes=attr(fit.Bessel,"SSErr")
ErPen=attr(fit.Pentaesf,"SSErr")
Variogramas=c(ErEs,ErEx,ErGau,ErLin,ErMat,ErBes,ErPen)
ErrVar =matrix(Variogramas,nrow=7,ncol=1,byrow=TRUE,
dimnames = list(c("Esferico", "Exponencial","Gaussiano","Lineal", "Matern", "Bessel","Pentaesferico"), c("Error")))
ErrVar## Error
## Esferico 1.3438136
## Exponencial 0.7530954
## Gaussiano 2.7085439
## Lineal 2.5396841
## Matern 0.7530954
## Bessel 0.9176738
## Pentaesferico 1.1284034En la figura 9 se puede observar gráficamente, la representación de los modelos exponencial, gaussiano, esférico y Bessel. El de menor error es el exponencial, seguido del Bessel y el esférico. El modelo gaussiano y lineal son los que menos se ajustan a los datos.
#Se grafica el variograma ajustado que fue el exponencial
plot(Vario.Var, pl=F,model=fit.Exponencial,main="Variograma ajustado con modelo Exponencial",xlab="Distancia",ylab="Semivarianza")
#se Ajusta al ojo para definir parametros iniciales de los diferentes modelos
vExp <- fit.variogram(Vario.Var, vgm(model = "Exp"),fit.method = 2)
plot ( Vario.Var , vExp,main="Variograma ajustado con modelo Exponencial")
vGau <- fit.variogram(Vario.Var, vgm(model = "Gau"),fit.method = 2)
plot ( Vario.Var , vGau,main="Variograma ajustado con modelo Gaussiano")
vSph <- fit.variogram(Vario.Var, vgm(model = "Sph"),fit.method = 2)
plot ( Vario.Var , vSph, main="Variograma ajustado con modelo Esferico")
vMat <- fit.variogram(Vario.Var, vgm(model = "Mat", nugget = 1,kappa = 0.5),fit.method = 2)
plot ( Vario.Var , vMat, main="Variograma ajustado con modelo Matern")
vBes <- fit.variogram(Vario.Var,vgm("Bes"),fit.method = 2)
plot ( Vario.Var , vBes, main="Variograma ajustado con modelo Bessel")
vSte <- fit.variogram(Vario.Var,vgm("Ste"),fit.method = 2)
plot ( Vario.Var , vSte)
Figura 9. Modelos de semivariograma.Fuente: Propia
4.5 Kriging Simple
Empleando el método de Kriging Simple se obtuvo un mapa pronóstico (Figura 12) y un mapa de desviaciones (Figura 13) que refleja menor error estándar en la predicción en comparación con el kriging ordinario, sin embargo, el hecho de asumir una media conocida, hace que se obtenga una varianza menor en los datos (Rodríguez, 2011).
A caption
En el departamento se evidencia que la precipitación tiene fuertes variaciones relacionadas con la topografía y el uso del suelo, se observan sectores con valores bajos (azul) en la zona del Urabá y el centro del departamento. Por otro lado, las altas precipitaciones se encuentran demarcadas en los sectores de Bajo Cauca y parte de Oriente y Magdalena Medio.
4.6 Kriging Ordinario
Se aplicó el método de interpolación kriging ordinario, obteniendo dos mapas: el primer mapa de pronóstico (figura 10) que permite visualizar el comportamiento de la precipitación en el departamento y el segundo mapa de predicción del error estándar (figura 11) en el que se identifican las zonas en las cuales hay mayor incertidumbre en cuanto a la predicción obtenida con la interpolación. En la figura 10 se observa que las zonas de mayor precipitación se encuentran en el bajo Cauca y el Urabá antioqueño, siendo este último el que presenta mayor incertidumbre en la predicción dado a que en esta zona se encuentran pocas estaciones pluviométricas. También coincide con los ojos de toro que se formaron con el método de interpolación IDW. Aun así los valores concordaron con los mapas de precipitaciones que genera el IDEAM.
A caption
#crear grilla
x1<-min(Preci_Ant$LONGITUD)
x2<-max(Preci_Ant$LONGITUD)
y1<-min(Preci_Ant$LATITUD)
y2<-max(Preci_Ant$LATITUD)
r=1
mat<-matrix(0,nrow=((x2-x1)/1+1)*((y2-y1)/1+1),ncol=2)
for (i in seq(x1,x2,by=1))
{for (j in seq(y1,y2,by=1))
{mat[r,1]<-i
mat[r,2]<-j
r=r+1}}
write.csv(mat,file="grid.csv")Preci_Ant2 <- shapefile('D:/Dropbox/Asig_maestria_Geomatica/Geoestadistica/Proyecto_final/Antio_precip.shp')
proj4string(Preci_Ant2) <- CRS("+init=epsg:3116")
data.grid<-read.csv("grid.csv")
coordinates(data.grid) <- c("V1", "V2")
proj4string(data.grid) <- CRS("+init=epsg:3116")
k<-krige(Precip_Anu~1,locations=Preci_Ant2,newdata=data.grid,model=vExp )write.csv(k,file="OK_sph.csv")
#display maps of concentrations
spplot(k,"var1.pred",asp=1,col.regions=bpy.colors(64),xlim=c(x1,x2),ylim=c(y1,y2),main="Predicción kriging ordinario con modelo esferico")5. Conclusiones
El método de interpolación de distancia inversa ponderada no arrojó buenos resultados en el pronóstico de la precipitación del departamento de Antioquia, esto sucedió en parte, por la ausencia de información en algunos sectores y la distancia entre los mismos. Este método al depender únicamente de la distancia genera dificultades en zonas con puntos muestrales de valores extremos y en zonas con baja cobertura de muestras; sin embargo resulta útil como base preliminar para identificar el comportamiento de la variable medida. Se puede decir que los métodos de interpolación geoestadísticos, en el caso de la precipitación del departamento de Antioquia, son los métodos más adecuados para obtener un mapa de pronóstico, resaltando el kriging ordinario ya que asumir una media conocida, como se hace en el kriging simple, con la distribución espacial de los pluviómetros que se tiene actualmente, puede inducir a un error. La ventaja del kriging ordinario sobre el kriging simple (Tabla 4), radica en permitir que la media varie localmente y de esta forma se generan variaciones mayores en puntos con medias superiores o inferiores a la global, lo que es adecuado para el departamento ya que su precipitación tiene un rango amplio (entre 1000 y 6.450 mm).
#La validación cruzada se realiza empleando 10-fold para medir los errores de ajuste de kriging universal,ordinario y simple
#### Ordinario
#Esferico con logaritmo
OESL.VC=krige.cv(log(Precip_Anual)~1,data,fit.Esférico,nfold=10)
#Gausiano con logaritmo
OGAL.VC=krige.cv(log(Precip_Anual)~1,data,fit.Gaussiano,nfold=10)
#exponencial con logaritmo
OEXL.VC=krige.cv(log(Precip_Anual)~1,data,fit.Exponencial,nfold=10)
#Simple
SESL.VC=krige.cv(log(Precip_Anual)~1,data,fit.Esférico,nfold=10,beta=5)
SGAL.VC=krige.cv(log(Precip_Anual)~1,data,fit.Gaussiano,nfold=10,beta=5)
SEXL.VC=krige.cv(log(Precip_Anual)~1,data,fit.Exponencial,nfold=10,beta=5)
#Transformar a hoja de datos
#### Ordinario
OrEsL=as.data.frame(OESL.VC)$residual
OrGAL=as.data.frame(OGAL.VC)$residual
OrEXL=as.data.frame(OEXL.VC)$residual
#Simple
SrEsL=as.data.frame(SESL.VC)$residual
SrGAL=as.data.frame(SGAL.VC)$residual
SrEXL=as.data.frame(SEXL.VC)$residual
#Mean Error(ME)
#### Ordinario
OMEsL=mean(OrEsL)
OMEGAL=mean(OrGAL)
OMEXL=mean(OrEXL)
#Simple
SMEsL=mean(SrEsL)
SMEGAL=mean(SrGAL)
SMEXL=mean(SrEXL)
MEsL=c(OMEsL,SMEsL)
MEGAL=c(OMEGAL,SMEGAL)
MEEXL=c(OMEXL,SMEXL)
#Root Mean Squared Error (RMSE)
#### Ordinario
O.RMSE.ESL=sqrt(mean(OrEsL^2))
O.RMSE.GAL=sqrt(mean(OrGAL^2))
O.RMSE.EXL=sqrt(mean(OrEXL^2))
#Simple
S.RMSE.ESL=sqrt(mean(SrEsL^2))
S.RMSE.GAL=sqrt(mean(SrGAL^2))
S.RMSE.EXL=sqrt(mean(SrEXL^2))
RMSE.ESL=c(O.RMSE.ESL,S.RMSE.ESL)
RMSE.GAL=c(O.RMSE.GAL,S.RMSE.GAL)
RMSE.EXL=c(O.RMSE.EXL,S.RMSE.EXL)
#Mean Standardized Prediction Error (MSPE)
#### Ordinario
O.MSPE.ESL=mean(OrEsL^2)
O.MSPE.GAL=mean(OrGAL^2)
O.MSPE.EXL=mean(OrEXL^2)
#Simple
S.MSPE.ESL=mean(SrEsL^2)
S.MSPE.GAL=mean(SrGAL^2)
S.MSPE.EXL=mean(SrEXL^2)
MSPE.ESL=c(O.MSPE.ESL,S.MSPE.ESL)
MSPE.GAL=c(O.MSPE.GAL,S.MSPE.GAL)
MSPE.EXL=c(O.MSPE.EXL,S.MSPE.EXL)
datos.ESl=c(MEsL,RMSE.ESL,MSPE.ESL)
datos.GAL=c(MEGAL,RMSE.GAL,MSPE.GAL)
datos.EXL=c(MEEXL,RMSE.EXL,MSPE.EXL)
#Se imprime una tabla donde se muestran los errores obtenidos por la validación cruzada
# Evaluación krigging esferico logaritmo
Evaluación.ESL=matrix(datos.ESl,nrow=3,ncol=2,byrow=TRUE,
dimnames = list(c("ME", "RMSE","MSPE"), c("Ordinario", "Simple")))
Evaluación.ESL
# Evaluación krigging Gaussiano con logaritmo
Evaluación.GAL=matrix(datos.GAL,nrow=3,ncol=2,byrow=TRUE,
dimnames = list(c("ME", "RMSE","MSPE"), c("Ordinario", "Simple")))
Evaluación.GAL
# Evaluación krigging exponenciaL CON LOGARITMO
Evaluación.EXL=matrix(datos.EXL,nrow=3,ncol=2,byrow=TRUE,
dimnames = list(c("ME", "RMSE","MSPE"), c("Ordinario", "Simple")))
Evaluación.EXLTeniendo en cuenta que la precipitación es una variable que se encuentra directamente relacionada con la altura, sería útil aplicar métodos de interpolación como el cokriging y evaluar la exactitud temática del mismo. La implementación de estos mapas puede ser un valioso recurso para la toma de decisiones, como evaluar zonas adecuadas de localización de cultivos, políticas ambientales y futura disposición de pluviómetros en el departamento.
6. Bibliografia
Abedini, M. & Nasseri, M. (2009). Inverse Distance Weighting revisited. Dept. of Civil and Environmental Engineering, Shiraz University, Shiraz, Iran.
Achilleos, Georgios. (2011). The Inverse Distance Weighted interpolation method and error propagation mechanism – creating a DEM from an analogue topographical map. Journal of Spatial Science. 56. 283-304. 10.1080/14498596.2011.623348.
Ahrens, B. (2006). Distance in spatial interpolation of daily rain gauge data. Hydrology and Earth System Sciences, 10, 197–208.
Algarra, P., Sanchez, X. . (2006). Métodos Estadísticos Aplicados. España: Universitat de Barcelona.
Boletin Temático. (2011). Anuario Estadístico. Medio Ambiente. Consulta: sábado, 04 de junio de 2022, sitio web: http ://Antioquia .gov.co/ PDF2 /boletin_2011_ medio _ambiente. pdf
Burrough, P., Mcdonnell, R. (1998). Principles of geographical information systems. Regional Science. New York: Oxford University Press, 333.
ESRI. (2013). Consulta: sábado, 04 de junio de 2022, ArcGIS Geostatistical Analyst, Sitio web:http://www.esri.es/es/productos/arcgis/arcgis-for-desktop/extensiones-arcgis-for-desktop/arcgis-geostatistical-analyst/.
ESRI. (2022). ArcGIS (Versión 10.5) [Software de procesamiento digital de imágenes satelitales]. Los Ángeles, Estados Unidos: Environmental Systems Research Institute, Inc.
FONADE. (2022) Antioquia, Consulta: sábado, 04 de junio de 2022, sitio web: http ://www. fonade.gov.co/ GeoTec/ inventario1 / zonas/ Antioquia.php.
Kleijnen, Jack. (2017). Kriging: Methods and Applications. SSRN Electronic Journal. 10.2139/ssrn.3075151.
Lu, G. & Wong, D. (2008). An adaptive inverse-distance weighting spatial interpolation technique. Computers & Geosciences, Vol.34, 1044–1055.
Gómez, M. (2005). Inferencia Estadística. España: Ediciones Díaz de Santos.
Li, J., & Heap, A. D. (2008). A Review of Spatial Interpolation Methods for Environmental Scientists. Australian Geological Survey Organisation (Vol. GeoCat# 68). http://doi.org/http://www.ga.gov.au/image_cache/GA12526.pdf
Lichtenstern, A. (2013). Kriging methods in spatial statistics. Technische Universitat Müchen Department of Mathematics. Garching, August 13
Lixin, L., Weitian, T. & Reinhard, P. (2020). Chapter 7: Spatiotemporal interpolation methods for air pollution. Spatiotemporal Analysis of Air Pollution and Its Application in Public Health. https://doi.org/10.1016/B978-0-12-815822-7.00007-8 .
Marquínez, J., Lastra, J., García, P., (2003). Estimation models for precipitation in mountainous regions: the use of GIS and multivariate analysis. Journal of Hydrology. 270, 1-11.
Mateu, J. & Morell, I. (2003). Geoestadística y Modelos Matemáticos en Hidrogeología . España: Publicacions de la Universitat Jaume I.
Paramasivam, C.R. & Venkatramanan, S. 2019. Chapter 3: An Introduction to Various Spatial Analysis Techniques. GIS and Geostatistical Techniques for Groundwater Science. https://doi.org/10.1016/B978-0-12-815413-7.00003-1.
R Development Core Team. (2008). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Sitio web: http://www.R-project.org
Rodríguez, J., (2010). Estimación de la Distribución Espacial de la Precipitación en Zonas Montañosas Mediante Métodos Geoestadísticos. Escuela Técnica Superior de Ingenieros de Caminos, Canales y Puertos.
Rojas, E., Arce,B.,Peña, A., Boshell, F. y Ayarza,M.. (2010). Cuantificación e interpolación de tendencias locales de temperatura y precipitación en zonas alto andinas de Cundinamarca y Boyacá (Colombia). Revista Corpoica, 11, 173-182.
Sánchez, F. J. (2022). Hidrología Superficial y Subterránea. 2ª ed. Kindle Direct Publishing., 440 pp.
Sánchez, P. (2006). Métodos estadísticos aplicados, Universidad de Barcelona, España, p. 81- 85
Zimmerman, D., Pavlik, C., Ruggles, A., Armstrong, M. P., (1999). An Experimental Comparison of Ordinary and Universal Kriging and Inverse Distance Weighting. Mathematical Geology, Vol. 31, No. 4, 375 - 390.
Li, Z., Kuo,W., Hao, M. and Yaoxiang, Wu. (2018). An Adjusted Inverse Distance Weighted Spatial Interpolation Method. 3rd International Conference on Communications, Information Management and Network Security. Advances in Computer Science Research, volume 65