#Geoestadistica
#Análisis variográfico
##############################################
##############################
# Instalar los paquetes necesarios
[Link]("readxl")
[Link]("ggplot2")
[Link]("sp")
[Link]("gstat")
[Link]("geoR")
[Link]("spdep")
[Link]("ncf")
[Link]("e1071")
##############################################
##############################
# Cargar las librerías necesarias
library(readxl)
library(ggplot2)
library(sp)
library(gstat)
library(geoR)
library(spdep)
library(ncf)
library(e1071)
##############################################
#############################
#Preparación de datos
# Función para la transformación de puntuaciones normales
nscore <- function(x) {
nscore <- qqnorm(x, [Link] = FALSE)$x # puntuación normal
[Link] <- [Link](x = sort(x), nscore = sort(nscore))
return(list(nscore = nscore, [Link] = [Link]))
}
# Leer los datos desde el archivo Excel
data <- read_excel("geo_data.xlsx")
# Convertir a data frame
data <- [Link](data)
plot(data) # muestra como están distribuidos sus datos en el espacio
##############################################
############
# Cargar librería dplyr
library(dplyr)
# Crear un resumen estadístico
stats_summary <- geo_data %>%
summarise(
Promedio = mean(copper, [Link] = TRUE), # Cambia 'copper' si
es necesario
Mínimo = min(copper, [Link] = TRUE), # Cambia 'copper' si es
necesario
Máximo = max(copper, [Link] = TRUE), # Cambia 'copper' si
es necesario
Desviación_Estandar = sd(copper, [Link] = TRUE) # Cambia 'copper'
si es necesario
)
# Imprimir el resumen estadístico
print(stats_summary)
##############################################
#######################
# Análisis exploratorio básico de datos
# Aplicar la transformación de puntuaciones normales de los datos de cobre
# Esta transformación hace que los datos tengan una distribución normal
[Link] <- nscore(data$copper) # Transformación de puntuaciones
normales
data[["Ncopper"]] <- [Link]$nscore # Añadir la transformación de
puntuaciones normales
# Configurar la visualización de gráficos en una matriz de 2x2
par(mfrow = c(2, 2), mar = c(4, 4, 2, 1)) # Ajustar márgenes para gráficos
# Histograma de cobre
hist(data$copper, main = "Cobre", xlab = "Cobre", col = "lightblue", breaks
= 15)
# Gráfico CDF de cobre
plot(ecdf(data$copper), main = "CDF de Cobre", xlab = "Cobre", ylab =
"Probabilidad Acumulada")
# Histograma de puntuaciones normales de cobre
hist(data$Ncopper, main = "N[Cobre]", xlab = "N[Cobre]", col = "lightblue",
breaks = 15)
# Gráfico CDF de puntuaciones normales de cobre
plot(ecdf(data$Ncopper), main = "CDF de N[Cobre]", xlab = "N[Cobre]", ylab
= "Probabilidad Acumulada")
# Verificar la asimetría
# Calcular skewness para cobre
skewness_copper <- skewness(data$copper, [Link] = TRUE)
cat("Skewness de Cobre: ", skewness_copper, "\n")
# Calcular skewness para Ncopper
skewness_ncopper <- skewness(data$Ncopper, [Link] = TRUE)
cat("Skewness de N[Cobre]: ", skewness_ncopper, "\n")
##############################################
#############
# Autocorrelación espacial
# Índices de Moran y Geary
# Convertir a data frame si es necesario
data <- [Link](data)
# Crear una matriz de coordenadas
coords_matrix <- cbind(data$x, data$y)
# Crear un objeto de vecindad (neighbors) y una matriz de pesos (weights)
# Usamos una distancia máxima de la mitad de la distancia máxima entre
puntos
nb <- spdep::dnearneigh(coords_matrix, 0, max(dist(coords_matrix)) / 2)
lw <- spdep::nb2listw(nb, style = "W")
# Calcular el índice de Moran
moran <- spdep::[Link](data$Ncopper, lw)
# Calcular el índice de Geary
geary <- spdep::[Link](data$Ncopper, lw)
# Imprimir los resultados
# Índice de Moran: El Índice de Moran mide la autocorrelación espacial
global, lo que significa que evalúa la relación entre valores en todo el
espacio.
# Moran I statistic: este es el valor que nos interesa, ya que es el resultado
del índice de Morán
# P-value: El valor p asociado con el índice de Moran indica la significancia
estadística. Un valor p menor que 0.05 sugiere que la autocorrelación
espacial observada es significativa.
# Un valor positivo de Moran's I indica autocorrelación positiva, mientras
que un valor negativo indica autocorrelación negativa. Un valor cercano a
cero indica que no hay autocorrelación espacial.
cat("Índice de Moran:\n")
print(moran)
# Índice de Geary: evalúa la autocorrelación espacial local o la variabilidad
entre vecinos cercanos, es más sensible
# Geary C statistic: este es el valor que nos interesa, ya que es el resultado
del índice de Geary
# P-value: El valor p asociado con el índice de Geary indica la significancia
estadística. Un valor p menor que 0.05 sugiere que la autocorrelación
espacial observada es significativa.
# Un valor de Geary menor que 1 indica autocorrelación positiva, mientras
que un valor mayor que 1 indica autocorrelación negativa.
cat("\nÍndice de Geary:\n")
print(geary)
##############################################
###########################################
# Distribución espacial
# Seleccionar la columna de interés: 'Ncopper'
XCoord <- data$x
YCoord <- data$y
Ncopper <- data$Ncopper
# Función para calcular la distribución espacial
DEspacial <- function(x, y, z, n_bins = 10, xlab = "X", ylab = "Y", zlab = "Z",
main = "Distribución Espacial") {
df <- [Link](x = x, y = y, z = z)
ggplot(df, aes(x = x, y = y, color = z)) +
geom_point() +
scale_color_viridis_c() +
labs(x = xlab, y = ylab, color = zlab, title = main) +
theme_minimal() +
theme([Link] = "bottom")
}
# Calcular y mostrar la distribución espacial
DEspacial(XCoord, YCoord, Ncopper, n_bins = 10,
xlab = 'X (m)', ylab = 'Y (m)', zlab = 'Ncopper', main = 'Distribución
Espacial de Ncopper')
##############################################
###############
# Análisis de Tendencia (Estacionariedad)
# Convertir los datos a un objeto SpatialPointsDataFrame
data_sp <- [Link](X = XCoord, Y = YCoord, ncopper = Ncopper)
coordinates(data_sp) <- ~ X + Y
# Extraer las coordenadas y datos
coords <- coordinates(data_sp) # Extrae las coordenadas
ncopper <- data_sp$ncopper # Extrae el dato de cobre transformado
# Convertir a un objeto geoR
geo_data <- [Link]([Link](coords, ncopper), [Link] = c("X",
"Y"), [Link] = "ncopper")
# Graficar los datos usando geoR
plot(geo_data)
title(main = "Distribución Espacial de Ncopper")
##############################################
####################
# Visualizar la tendencia superficial
# Estacionariedad - Interpretación
# No-estacionariedad: si el valor de p es menor que 0.05 implica no-
estacionariedad
# Estacionariedad: si el valor de p es mayor que 0.05 implica
estacionariedad
# Crear una malla para la visualización
trend_surface <- [Link](x = seq(min(data$x), max(data$x), [Link]
= 100),
y = seq(min(data$y), max(data$y), [Link] = 100))
# Ajustar el modelo lineal de tendencia
trend_model <- lm(Ncopper ~ x + y, data = data)
summary(trend_model)
# Predecir los valores en la superficie de tendencia
trend_surface$Ncopper <- predict(trend_model, newdata = trend_surface)
# Crear el gráfico de la superficie de tendencia
ggplot(trend_surface, aes(x = x, y = y, fill = Ncopper)) +
geom_tile() +
scale_fill_viridis_c() +
labs(title = "Superficie de Tendencia", x = "x", y = "y", fill = "Ncopper") +
theme_minimal()
ggplot() +
geom_tile(data = trend_surface, aes(x = x, y = y, fill = Ncopper)) +
geom_contour(data = trend_surface, aes(x = x, y = y, z = Ncopper), color
= "white") +
geom_point(data = data, aes(x = x, y = y), color = "black", size = 2) +
scale_fill_viridis_c() +
labs(title = "Superficie de Tendencia con Datos Originales", x = "X (m)", y
= "Y (m)", fill = "Ncopper") +
theme_minimal()
# Convertir los datos a un objeto SpatialPointsDataFrame
coordinates(data) <- ~x + y
proj4string(data) <- CRS("+proj=utm +zone=17 +south +datum=WGS84")
# Actualiza el CRS según sea necesario
# Crear el gráfico de burbujas
bubble(data, "Ncopper",
col = c("#00ff0088", "#00ff0088"), # Color de las burbujas
main = "Concentraciones de Cobre (ppm)") # Título del gráfico
##############################################
###########################
# Detrending
# El detrending se está haciendo con un modelo polinómico de 3er orden,
porque las tendencias de los datos son no lineales
# El detrending se puede hacer con modelos polinómicos de 1er orden
(lineal), de 2do orden, 3er orden, y órdenes superiores
# Preparación de datos
data <- [Link](data)
# Ajustar un modelo polinómico de tercer grado
poly_model_3rd <- lm(Ncopper ~ poly(x, 3) + poly(y, 3), data = data)
summary(poly_model_3rd)
# Esta parte es para que sepan como se usan los modelos de detrendado,
pero vamos a trabajar con los datos Ncopper
##############################################
#############################
# Análisis variográfico
# Nube del variograma
# Convertir los datos a un objeto geoR
geo_data <- [Link]([Link](X = data$x, Y = data$y, Z =
data$Ncopper),
[Link] = c("X", "Y"), [Link] = "Z")
# Calcular la nube del variograma para los datos
variog_result <- variog(geo_data, option = "cloud")
# Graficar la nube del variograma
plot(variog_result, main = "Variograma Cloud")
# Calcular el variograma cloud con estimación modulus
# Explicación: La estimación modulus es una técnica que puede ser más
robusta frente a outliers. En lugar de usar la diferencia directa entre los
pares de puntos, utiliza el módulo, que suaviza el efecto de los outliers.
p_cloud.m <- variog(geo_data, option="cloud", est="modulus")
plot(p_cloud.m, main = "Variograma Cloud con Estimación Modulus")
# Calcular el variograma cloud con lambda = 1 (transformación simple)
# Explicación: La transformación lambda = 1 es una transformación simple
que puede ayudar a estabilizar la variabilidad en los datos. Aunque no es
tan común como la transformación logarítmica (lambda = 0), puede ser útil
cuando se tiene una mezcla de valores positivos y negativos. En este
ejemplo, los datos detrendados tienen valores positivos y negativos.
p_cloud_lambda1 <- variog(geo_data, option="cloud", lambda=1)
plot(p_cloud_lambda1, main = "Variograma Cloud con Transformación
Lambda=1")
# Sustentar el uso de la transformación lambda=1
# Verificar valores no positivos en los datos
sum(data$Ncopper <= 0) # Contar valores no positivos
sum(data$data <= 0) # Contar valores no positivos
##############################################
##############################
##############################################
##############################
# Ajustar modelos de variograma con parámetros iniciales
model_spherical <- vgm(model = "Sph", range = range_value, nugget =
nugget_value, psill = sill_partial)
model_exponential <- vgm(model = "Exp", range = range_value, nugget =
nugget_value, psill = sill_partial)
model_gaussian <- vgm(model = "Gau", range = range_value, nugget =
nugget_value, psill = sill_partial)
# Ajustar los modelos a los datos empíricos
fit_spherical <- [Link](v_empirical, model_spherical)
fit_exponential <- [Link](v_empirical, model_exponential)
fit_gaussian <- [Link](v_empirical, model_gaussian)
# Imprimir resultados de modelos ajustados
cat("Modelo Esférico Ajustado:\n")
print(fit_spherical)
cat("Modelo Exponencial Ajustado:\n")
print(fit_exponential)
cat("Modelo Gaussiano Ajustado:\n")
print(fit_gaussian)
# Tremendo código para 1 gráfico!!!!
# Graficar el variograma con los modelos ajustados
# Definir funciones de variograma teórico con parámetros fijos
spherical_variogram <- function(dist, range, nugget, psill) {
semivariance <- ifelse(dist <= range,
nugget + psill * (1.5 * (dist / range) - 0.5 * (dist / range)^3),
nugget + psill)
return(semivariance)
}
exponential_variogram <- function(dist, range, nugget, psill) {
semivariance <- nugget + psill * (1 - exp(-dist / range))
return(semivariance)
}
gaussian_variogram <- function(dist, range, nugget, psill) {
semivariance <- nugget + psill * (1 - exp(- (dist^2) / (range^2)))
return(semivariance)
}
# Parámetros ajustados para cada modelo
# Modelo Esférico
range_spherical <- 784.1718 # Reemplazado por el valor ajustado
nugget_spherical <- 0.1584315 # Reemplazado por el valor ajustado
psill_spherical <- 0.9805530 # Reemplazado por el valor ajustado
# Modelo Exponencial
range_exponential <- 354.5204 # Reemplazado por el valor ajustado
nugget_exponential <- 0.04278656 # Reemplazado por el valor ajustado
psill_exponential <- 1.18343083 # Reemplazado por el valor ajustado
# Modelo Gaussiano
range_gaussian <- 332.6261 # Reemplazado por el valor ajustado
nugget_gaussian <- 0.2665510 # Reemplazado por el valor ajustado
psill_gaussian <- 0.8231943 # Reemplazado por el valor ajustado
# Cambiar el metal de cadmio a cobre para el análisis
# Obtener los límites de los ejes del variograma empírico
xlims <- range(v_empirical$dist)
ylims <- range(v_empirical$gamma)
# Crear una secuencia de distancias para calcular la semivarianza teórica
dist_seq <- seq(from = min(xlims), to = max(xlims), [Link] = 100)
# Calcular la semivarianza teórica para cada modelo
spherical_values <- spherical_variogram(dist_seq, range_spherical,
nugget_spherical, psill_spherical)
exponential_values <- exponential_variogram(dist_seq, range_exponential,
nugget_exponential, psill_exponential)
gaussian_values <- gaussian_variogram(dist_seq, range_gaussian,
nugget_gaussian, psill_gaussian)
# Verificar longitudes (deberían ser 100 cada uno)
print(length(dist_seq)) # Debería ser 100
print(length(spherical_values)) # Debería ser 100
print(length(exponential_values)) # Debería ser 100
print(length(gaussian_values)) # Debería ser 100
# Ajustar los parámetros gráficos
par(mfrow = c(1, 1)) # Configura para un solo gráfico
par(mar = c(5, 5, 4, 2) + 0.1) # Ajusta los márgenes (c(bottom, left, top,
right))
# Crear el gráfico del variograma empírico
plot(v_empirical$dist, v_empirical$gamma, main = "Variograma Empírico
para Cobre con Modelos Ajustados",
xlab = "Distancia", ylab = "Semivarianza", pch = 20, col = "black", cex =
1.2,
xlim = xlims, ylim = ylims)
# Añadir los modelos ajustados al gráfico
lines(dist_seq, spherical_values, col = "red", lwd = 2, lty = 1)
lines(dist_seq, exponential_values, col = "green", lwd = 2, lty = 2)
lines(dist_seq, gaussian_values, col = "purple", lwd = 2, lty = 3)
# Añadir una leyenda para distinguir los modelos
legend("bottomright", legend = c("Empírico", "Esférico", "Exponencial",
"Gaussiano"),
col = c("black", "red", "green", "purple"), lwd = 2, lty = c(NA, 1, 2, 3),
pch = c(20, NA, NA, NA), bg = "white", cex = 0.8,
bty = "n",
[Link] = strwidth("Empírico") * 1.5)
##############################################
#########################
# Selección del mejor modelo para el metal cobre
# Calcular la semivarianza teórica para cada modelo en las distancias del
variograma empírico
spherical_values_empirical <- spherical_variogram(v_empirical$dist,
range_spherical, nugget_spherical, psill_spherical)
exponential_values_empirical <- exponential_variogram(v_empirical$dist,
range_exponential, nugget_exponential, psill_exponential)
gaussian_values_empirical <- gaussian_variogram(v_empirical$dist,
range_gaussian, nugget_gaussian, psill_gaussian)
# Calcular el SSE para cada modelo
sse_spherical <- sum((v_empirical$gamma - spherical_values_empirical)^2,
[Link] = TRUE)
sse_exponential <- sum((v_empirical$gamma -
exponential_values_empirical)^2, [Link] = TRUE)
sse_gaussian <- sum((v_empirical$gamma - gaussian_values_empirical)^2,
[Link] = TRUE)
# Determinar el modelo con el menor SSE
sse_values <- c(Spherical = sse_spherical, Exponential = sse_exponential,
Gaussian = sse_gaussian)
# Mostrar el SSE para cada modelo
print(sse_values)
# Determinar el modelo con el menor SSE
best_model <- names([Link](sse_values))
# Imprimir el modelo seleccionado
cat("El mejor modelo basado en el método de mínimos cuadrados para
cobre es:", best_model, "\n")
#Geoestadistica
#Determinar áreas de predicción
##############################################
###############################
# Cargar las librerías necesarias
library(readxl)
library(ggplot2)
##############################################
########
#Convex hull
# Esto nos permite identificar el área donde haremos la interpolación
# Leer el archivo de Excel
geo_data <- read_excel("geo_data.xlsx")
# Calcular el convex hull
hull_indices <- chull(geo_data$x, geo_data$y)
hull_coords <- geo_data[hull_indices, ]
# Plotear los puntos y el convex hull
ggplot(geo_data, aes(x = x, y = y)) +
geom_point() +
geom_polygon(data = hull_coords, aes(x = x, y = y), fill = NA, color =
"red") +
labs(title = "Convex Hull de las Coordenadas Proporcionadas", x = "x", y =
"y") +
theme_minimal()
# Imprimir las coordenadas del convex hull
print("Coordenadas que forman el polígono (Convex Hull):")
print(hull_coords)