Variograma Experimental en R

---

En este tutorial revisaremos de forma breve y resumida, algunas de las propiedades vectoriales del plano y funciones trigonometricas usadas en la creación de la banda del variograma (Figura1). Esta revision será realizada textualmente y graficamente, quedando para el lector la verificación de las mismas. Finalmente, un código del computo del variograma experimental sera proporcionado en el lenguaje de programación R project.

Figura 1. Banda del Variograma

El variograma experimental de una variable regionalizada quantifica la dependencia espacial en una dirección especifica en función de la distancia. Existiendo una disimulitud espacial a mayores distancias entre las muestras, denominada amplitud. En la práctica el variograma experimental para datos irregularmente espaciados es calculado asumiendo algun grado de tolerancia (Figura1), como:

A continuación describiremos en tres, principales pasos, como el variograma experimental (Ecuación anterior) es calculado, usando las propiedades entre vectores y funciones trigonometricas.

El primer paso en el cálculo del variograma experimental es la definición de una o variás direcciones. Esta dirección es representada de forma vectorial por el par ordenado (x,y), donde cada componente del vector es calculado usando las coordenadas polares y el sistema acimutal (Figura 2).

Figura 2. Representacion del vector u, en coordenadas polares.

El segundo paso en el cálculo del variograma experimental, de forma computacional, es calcular la altura del cono (coneh) usando la función trigonométrica tangente del ángulo (tol), como (Figura 3):

Figura 3. Altura del cono (coneh).

El tercer y último paso es computar el vector componente entre los puntos p1 y p2, lo cual resulta en el vector v. Este nuevo vector es proyectado sobre el vector unitario, u, usando el producto punto escalar, con el objetivo de verificar si el número real obtenido, es menor o igual o mayor que la altura del cono (coneh) (Figura 4). Además de lo anterior, la diferencia entre las magnitudes de los vectores u y v son calculadas, con el objetivo de identificar si esta diferencia no excede la largura de banda (bandh) (Figura 4).

Figura 4. Aceptación o rechazo de la diferencia y proyección entre pares de puntos.

A seguir el código en R, del cálculo de variograma experimental. Así como el link de descarga del archivo de prueba.

library(gridExtra) library(ggplot2) data = read.table(file = "data.csv", header = TRUE, sep = ";") # Convirtiendo data.frame to matrix data = unname(as.matrix(data, dimnames = NULL)) # Parametros de tolerancias lagtol = 1.25 lag = 2.5 numlags = 12 # Parametros del vector direccional ndir = 2 azimuth = c(45, 135)*pi/180 tolerance = rep(45, ndir)*pi/180 bandwidth = rep(500, ndir) # Variable de interese y sus coordenadas coord = data[,c(1,2)] variable = data[,c(3)] # Definiendo tamaño de los vectores para el calculo del variograma # Experimental npairs = matrix(rep(0, numlags*ndir), nrow = ndir) mcorr = matrix(rep(0, numlags*ndir), nrow = ndir) lags = lag + seq(0,numlags-1,1)*lag for(dir in 1:ndir){ if(tolerance[dir]>=pi/2){ coneh = -1 bandwidthsq = -1 }else{ coneh = bandwidth[dir]/tan(tolerance[dir]) bandwidthsq = coneh*tan(tolerance[dir]) bandwidthsq = bandwidthsq*bandwidthsq } # Sistema de coordenadas polares y el sistema cimutal dirx = cos(pi/2 - azimuth[dir]) diry = sin(pi/2 - azimuth[dir]) # Calculando la correlacion espacial entre cada punto dentro del area for(i in 1:(length(variable)-1)){ for(j in (i+1):length(variable)){ v = coord[i,]-coord[j,] if(tolerance[dir]>=pi/2){ bool = TRUE }else{ h = dirx*v[1] + diry*v[2] d2 = sum(v^2) - h*h if(h>coneh){ bool = d2<bandwidthsq }else{ dmax = h*tan(tolerance[dir]) bool = d2<(dmax*dmax) } } if(!bool){ next } d = sqrt(sum(v^2)) lagid = -1 for(k in 1:numlags){ if(((lags[k]-lagtol)<=d) && (d<(lags[k]+lagtol))){ lagid = k } } if(lagid<0){ next } hpropi = variable[i] hpropj = variable[j] tpropi = variable[i] tpropj = variable[j] mcorr[dir,lagid] = mcorr[dir,lagid] + 0.5*(hpropi-hpropj)*(tpropi-tpropj) npairs[dir,lagid] = npairs[dir,lagid] + 1 } } mcorr[dir,] = mcorr[dir,]/npairs[dir,] } g1 = ggplot() + geom_point(aes(x=lags, y=mcorr[1,]), color = "red") + xlim(c(0,31)) + ylim(c(0,65)) + xlab("Distancia") + ylab("Variograma") g2 = ggplot() + geom_point(aes(x=lags, y=mcorr[2,]), color = "blue") + xlim(c(0,31)) + ylim(c(0,65)) + xlab("Distancia") + ylab("Variograma") grid.arrange(g1, g2, ncol = ndir)