Taller: Análisis estructural
modelado del variograma empírico
Análisis estructural de variables regionalizadas: Semivariograma empírico
El semivariograma es una técnica de gran importancia como paso exploratorio y previo a la estimación de valores desconocidos sobre un área de estudio, este suele utilizarse para identificar el nivel de variación espacial de un valor \(Z\) en diferentes puntos sobre el espacio, estos están representados por \(\mathbf{s}_i\).
La construcción del semivariograma experimental está dada por la expresión:
\[\hat{\gamma}(h) = \frac{1}{2N(h)} \sum_{i=1}^{N(h)} \left[ Z(\mathbf{s}_i) - Z(\mathbf{s}_i + \mathbf{h}) \right]^2\] en donde:
\(\hat{\gamma}(h)\) es el valor del semivariograma experimental para una distancia dada.
\(N(h)\) es el número de pares de puntos separados por la distancia.
\(Z(\mathbf{s}_i)\) es el valor de la variable en el punto espacial.
\(Z(\mathbf{s}_i + h)\) es el valor de la variable en el punto espacial desplazado.
\(h\) es la distancia entre los puntos considerados.
A continuación se detallan algunos procedimientos que conducen a su construcción e interpretación. En primer lugar es necesario conocer y graficar los valores sobre un espacio coordenado.
Datos espaciales
Para el ejemplo se toma una matriz \(A\) de \(m=6 \times n=6\) con espaciado 1, en donde cada espacio constituye un valor muestral o realización de la variable por analizar.
\[A \in \mathbb{R}^{m \times n}\] Para este caso la matriz de valores será:
\[\begin{bmatrix} 26 & 22 & 19 & 14 & 16 & 19 \\ 23 & 20 & 17 & 20 & 14 & 23 \\ 22 & 17 & 18 & 19 & 18 & 25 \\ 21 & 15 & 20 & 18 & 20 & 20 \\ 19 & 18 & 15 & 15 & 18 & 23 \\ 18 & 16 & 10 & 16 & 14 & 18 \\ \end{bmatrix}\]
# Crear valores espaciados en matriz A de m*n
matriz <- matrix(c(
26, 22, 19, 14, 16, 19,
23, 20, 17, 20, 14, 23,
22, 17, 18, 19, 18, 25,
21, 15, 20, 18, 20, 20,
19, 18, 15, 15, 18, 23,
18, 16, 10, 16, 14, 18
), nrow = 6, byrow = T)El espacio coordenado de esta variable se asume como su número de fila y columna en la matriz de entrada, iniciando el conteo desde la parte superior izquierda. Para analizar estos valores como una variable regionalizada es necesario declarar como variables independientes las coordenadas \(X\) y \(Y\)
# Obtener las dimensiones de la matriz
num_filas_y <- nrow(matriz)
num_columnas_x <- ncol(matriz)
# Crear el DataFrame con las coordenadas y los valores
df <- data.frame(x = rep(1:num_columnas_x, each = num_filas_y), # X= Número de columna
y = rep(1:num_filas_y, times = num_columnas_x), #Y= Número de fila
valor = as.vector(matriz) # Aplana la matriz en un vector
)
# Imprimir las primeras 5 filas
knitr::kable(df[1:5,],
caption = "Valores y coordenadas para variable de análisis",
align=c('c','c','c'),col.names = c('Coordenada X','Coordenada Y','Valor puntual'))| Coordenada X | Coordenada Y | Valor puntual |
|---|---|---|
| 1 | 1 | 26 |
| 1 | 2 | 23 |
| 1 | 3 | 22 |
| 1 | 4 | 21 |
| 1 | 5 | 19 |
Dado que el valor coordenado inicial está en la parte superior izquierda de la matriz, esto generará una inversión del orden de los valores a la hora de espacializarlo, esto se hace evidente al comparara el gráfico que usa las coordenadas para espacializar la variable, frente al orden que presenta la matriz de entrada.
# Crear objeto espacial preservando coordenadas en columnas
df2 = st_as_sf(df,coords=c('x','y'))
df2 <- cbind(df2,st_coordinates(df2))
# Graficar los puntos con ggplot2
# Crear paleta de colores
fun_color_range <- colorRampPalette(c("yellow3", "red")) # Create color generating function
# Crear el gráfico en ggplot2
plot_gg <- ggplot(data = df2) +
geom_point(aes(x = X, y = Y, color=valor), size = 3) + # valor representa la columna de valores que quieres representar
scale_colour_gradientn(colors = fun_color_range(16))+
geom_text(aes(x = X, y = Y, label = valor),
nudge_y = .3, size = 3, color = "black") + # Muestra los valores como etiquetas sobre los puntos
labs(title = "Distribución de valores en espacio coordenado",
x = "Coordenada X", y = "Coordenada Y",
color = "Valor") +
theme_gray() +
theme(legend.position = "bottom")
# Convertir el gráfico de ggplot en interactivo con plotly
ggplotly(plot_gg)Una vez espacializada la variable, es posible indagar someramente sobre algunas características de los valores en el espacio, permitiendo explorar si existe o no una tendencia en los datos, fenómeno que puede influir en el cálculo del semivariograma y por ende, en el modelo a implementar a la hora de estimar valores desconocidos.
#
# Crear objeto espacial para calcular gamma (h)
#
# crear objeto geodata
df = as.geodata(obj = df,coords.col = 1:2,data.col = 3)
# Graficar valores valor y su relación con las coordenadas
plot(df,trend = 'cte',lowess = T)
Una exploración visual sencilla permite identificar la presencia de tendencias entre el valor y sus coordenadas, sin que esta sea lineal. En estos casos puede ser recomendable buscar estimar la tendencia para sustraerla y trabajar con los valores residuales, con el fin de evitar que esto se inmiscuya en las estimaciones y genere subestimaciones o sobrestimaciones de los valores de interés.
Sin embargo, en este caso el cálculo de la semivarianza se realiza con los datos en su estado original.
Semivariograma
Una vez espacializada la variable es posible indagar sobre algunos aspectos básicos a tener en cuenta en la construcción del semivariograma, tales como la distribución estadística de la variable numérica regionalizada, y el alcance máximo de la distancia para el área de trabajo. Esta información puede ser útil a la hora de interpretar los primeros resultados del semivariograma, o al momento de definir los intervalos de distancia para la búsqueda de pares en el cálculo de la semivarianza.
## Number of data points: 36
##
## Coordinates summary
## x y
## min 1 1
## max 6 6
##
## Distance summary
## min max
## 1.000000 7.071068
##
## Data summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.0 16.0 18.0 18.5 20.0 26.0Para este caso de han definido 4 \(h\) para la búsqueda, estos son:
- 1, 2, 3 y 4
de igual forma estos \(\hat{\gamma}(h)\) se calcularán en diferentes direcciones, usando para ello el los siguientes valores azimut en grados:
- 0° - 45° - 90° y 135°
# Calcular valores de gamma(h) para los casos en los que h vale: 1, 2, 3 y 4
# Dirección: Este Oeste - Az:90
dir_90=variog(df,
unit.angle = "degrees",
direction = 90,
uvec = seq(1, 4, by = 1))
# Dirección: Norte Sur - Az:0
dir_0=variog(df,
unit.angle = "degrees",
direction = 0,
uvec = seq(1, 4, by = 1))
# Dirección: Sur/Occidente- Nor/Oriente - Az:45
dir_45=variog(df,
unit.angle = "degrees",
direction = 45,
pairs.min = 2,
uvec = c(1.414214,2.828427,4.242641,5.656854))
# Dirección: Nor/Occidente - Sur/Oriente - Az:135
dir_135=variog(df,
unit.angle = "degrees",
direction = 135,
pairs.min = 2,
uvec = c(1.414214,2.828427,4.242641,5.656854))Una vez calculadas las matrices de semivarianza para cada \(h\) y dirección es posible resumir los resultados en una tabla de comparación
#
# Crear tabla con resultados
#
h = as.vector(cbind(dir_90$uvec,dir_0$uvec,dir_45$uvec,dir_135$uvec))
gamma = as.vector(cbind(dir_90$v,dir_0$v,dir_45$v,dir_135$v))
# Redondear valores a dos decimales
gamma = round(gamma,2)
directions = as.vector(cbind(rep("Az:90",4),rep("Az:0",4),rep("Az:45",4),rep("Az:135",4)))
n = as.vector(cbind(dir_90$n,dir_0$n,dir_45$n,dir_135$n))
# Crear Dataframe
df_variogram <- data.frame(
h = h, # Distancias (h)
gamma = gamma, # Valores de semivariograma (gamma)
directions = directions, # Direcciones (Azimut)
n_pairs = n # Número de pares asociados a una distancia h y dirección Azimut
)
# Igualaes valores de Gamma a los h teŕicos utilizados, sin repetir
df_variogram <- df_variogram %>%
group_by(h, directions,n_pairs) %>%
summarize(gamma = mean(gamma), .groups = 'drop')
# Convertir el dataframe a formato de matriz donde las columnas sean las direcciones
res_gamma <- df_variogram %>%
pivot_wider(names_from = directions, values_from = gamma)
# Eliminar repeticiones de n
res_gamma = res_gamma |>
group_by(h) |>
summarize(
n = sum(n_pairs, na.rm = TRUE),
gamma_Az_0 = mean(`Az:0`, na.rm = TRUE),
gamma_Az_90 = mean(`Az:90`, na.rm = TRUE),
gamma_Az_135 = mean(`Az:135`, na.rm = TRUE),
gamma_Az_45 = mean(`Az:45`, na.rm = TRUE))
# Eliminar 'NaN' por espacio vacio
res_gamma = res_gamma %>%
mutate(across(everything(), ~ ifelse(is.na(.), "", .)))
# Ver resultados en tabla
knitr::kable(res_gamma,
caption = "Cálculo de gamma(h) en diferentes direcciones",
align=c('c','c','c','c','c','c'),col.names = c('h','N','Az 0°','Az 90°','Az 135°','Az 45°'),
digits = 1)| h | N | Az 0° | Az 90° | Az 135° | Az 45° |
|---|---|---|---|---|---|
| 1.0 | 30 | 4.6 | 8.22 | ||
| 1.4 | 25 | 6.6 | 9.22 | ||
| 2.0 | 24 | 8.06 | 9.44 | ||
| 2.8 | 16 | 11.22 | 13.88 | ||
| 3.0 | 18 | 9.11 | 14.06 | ||
| 4.0 | 12 | 9.71 | 14.62 | ||
| 4.2 | 9 | 20.11 | 12.11 | ||
| 5.7 | 4 | 9.38 | 18.75 |
Los valores de \(h\) cambian necesariamente para los azimut de 45° y 135°, esto dado que la distancia mínima entre los puntos, en estas direcciones, corresponde al valor de la hipotenusa. En estos caso \(h_1\) será:
\[h_{1} = \sqrt{\Delta x^2 + \Delta y^2}\\ h1 = 1.414214\]
La tabla resumen separa estos valores para identificar los valores de \(\gamma\) y los \(N\) involucrados en cada uno de los cálculos con las distancias y direcciones dadas. El valor de \(n\) será de gran ayuda a la hora de estimar un alcance adecuado, pues un \(\gamma\) calculado a partir de pocas observaciones puede no ser necesariamente representativo de la varianza en el espacio de los datos analizados.
A partir de estos datos es posible graficar los semivariogramas empíricos para cada una de las direcciones analizadas.
# Graficar variogramas
# Extraer las distancias y los valores de cada variograma
data_45 <- data.frame(h = dir_45$u, gamma = dir_45$v, direction = "45 grados")
data_135 <- data.frame(h = dir_135$u, gamma = dir_135$v, direction = "135 g rados")
data_0 <- data.frame(h = dir_0$u, gamma = dir_0$v, direction = "0 grados")
data_90 <- data.frame(h = dir_90$u, gamma = dir_90$v, direction = "90 grados")
# Combinar los datos en un solo dataframe
df_variograms <- rbind(data_45,data_135, data_0, data_90)
# Crear el gráfico con ggplot2
vario_multiple = ggplot(df_variograms, aes(x = h, y = gamma, color = direction)) +
geom_line(size = 1,linetype = 'dashed') + # Añadir líneas para cada dirección
geom_point()+
theme_minimal() +
labs(title = 'Variogramas en diferentes direcciones',
x = 'h',
y = 'Gamma (h)',
color = 'Dirección') +
scale_color_manual(values = c('45 grados' = 'green','135 grados' = 'cyan', '0 grados' = 'yellow', '90 grados' = 'purple')) +
theme_get()+
theme(legend.position = "topright")
# Convertir el gráfico de ggplot en interactivo con plotly
ggplotly(vario_multiple)Los resultados obtenidos en cada uno de los variogramas permiten concluir que, en este caso el fenómeno se comporta de forma anisotropica sobre es espacio, presentando valores altos y diferenciales de semivarianza en función de la dirección Azimuth 135°. Sin embargo es necesario tener en cuenta el número de pares utilizados para su cálculo, pues además de ser sensible a datos atípicos, una semivarianza que considera pocos puntos puede no ser representativa de la variación del fenómeno sobre el espacio.