Carga de librerías

# ========================
# 1. Instalación de paquetes (solo si es necesario)
# ========================
# install.packages("sp")
# install.packages("gstat")
# install.packages("raster")
# install.packages("geoR")

library(sp)
library(gstat)
library(raster)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(geoR)
## Warning in fun(libname, pkgname): couldn't connect to display ":0"
## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.9-4 (built on 2024-02-14) is now loaded
## --------------------------------------------------------------

Cargar datos

#setwd("/cloud/project")  # Cambiar si es necesario
df <- read.table("DatosCaso1.txt", header = TRUE)
df_espacial <- df %>% select(x, y, Cu)
coordinates(df_espacial) <- ~x + y
class(df_espacial)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

Análisis exploratorio

par(mfrow = c(3,2))
plot(df_espacial, asp = 1, pch = 19, main = "A", col = heat.colors(100)[cut(df_espacial$Cu, 100)])
plot(df$x, df$Cu, main = "B", pch = 19, xlab = "Coordenada x", ylab = "Cu")
plot(df$Cu, df$y, pch = 19, main = "C", xlab = "Cu", ylab = "Coordenada y")
plot(density(df$Cu), main = "D")
boxplot(df$Cu, main = "E")

Semivariograma empírico

df_espacial_geo <- as.geodata(df_espacial)
cloud1 <- variog(df_espacial_geo, option = "cloud", max.dist = 1)
## variog: computing omnidirectional variogram
plot(cloud1, xlab = "Distancia", ylab = "Semivarianza")

bin1 <- variog(df_espacial_geo, uvec = seq(0, 1, length.out = 11))
## variog: computing omnidirectional variogram
plot(bin1, pch = 19, type = "b", xlab = "Distancia", ylab = "Semivarianza")

variogram_empirical <- variogram(Cu ~ 1, df_espacial)
plot(variogram_empirical, pch = 19)

Ajuste de modelos teóricos

modelo_teorico <- vgm(psill = 350, model = "Sph", range = 1.7, nugget = 250)

smooth <- variog(df_espacial_geo, option = "smooth", max.dist = 1, n.points = 100, kernel = "normal", band = 0.2)
## variog: computing omnidirectional variogram
plot(variog(df_espacial_geo, max.dist = 1), main = "", xlab = "Distancia", ylab = "Semivarianza", pch = 19, ylim = c(0,800))
## variog: computing omnidirectional variogram
lines.variomodel(cov.model = "exp", cov.pars = c(350, 1.7),  nugget = 250, max.dist = 2.2, lwd = 3)
lines.variomodel(cov.model = "sph", cov.pars = c(350, 1.7), nugget = 250, max.dist = 1, lwd = 3, col = 2)
lines.variomodel(cov.model = "gau", cov.pars = c(350, 1.7), nugget = 250, max.dist = 1, lwd = 3, col = 3)
lines(smooth, type = "l", lty = 2)
legend("bottomright", c("Empírico", "Modelo Exponencial", "Suavizado", "Smooth"), lty = c(1,1,1,2), lwd = c(3,3,3), col = c(1,2,3))

Ajuste modelo exponencial (geoR)

variograma <- variog(df_espacial_geo, max.dist = 1)
## variog: computing omnidirectional variogram
plot(variograma, main = "Ajuste de semivariograma exponencial", xlab = "distance", ylab = "semivariance")
lines.variomodel(cov.model = "exp", cov.pars = c(350, 1.7), nugget = 250, max.dist = 2.2, lwd = 3)

Análisis direccional

vario72 <- variog(df_espacial_geo, max.dist = 1, direction = pi/2.5)
## variog: computing variogram for direction = 72 degrees (1.257 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
vario.4 <- variog4(df_espacial_geo, max.dist = 1)
## 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)
## variog: computing omnidirectional variogram
plot(vario.4, lwd = 2, xlab = "Distancia", ylab = "Semivarianza")

Kriging y predicción espacial

grid_pred <- expand.grid(x = seq(min(df$x), max(df$x), by = 0.05),
                         y = seq(min(df$y), max(df$y), by = 0.05))
coordinates(grid_pred) <- ~x + y
gridded(grid_pred) <- TRUE

variograma_empirico <- variogram(Cu ~ 1, df_espacial, max.dist = 1)
## Warning in variogram.default(y, locations, X, trend.beta = beta, grid = grid, :
## the following arguments are ignored: 1
model_exp <- vgm(psill = 350, model = "Exp", range = 1.7, nugget = 250)
model_sph <- vgm(psill = 350, model = "Sph", range = 1.7, nugget = 250)

fitted_exp <- fit.variogram(variograma_empirico, model_exp)
## Warning in fit.variogram(variograma_empirico, model_exp): No convergence after
## 200 iterations: try different initial values?
fitted_sph <- fit.variogram(variograma_empirico, model_sph)
## Warning in fit.variogram(variograma_empirico, model_sph): singular model in
## variogram fit
kriging_exp <- krige(Cu ~ 1, df_espacial, grid_pred, model = fitted_exp)
## [using ordinary kriging]
kriging_sph <- krige(Cu ~ 1, df_espacial, grid_pred, model = fitted_sph)
## [using ordinary kriging]
spplot(kriging_exp, "var1.pred", main = "Kriging Exponencial")

spplot(kriging_sph, "var1.pred", main = "Kriging Esférico")

Validación cruzada

df_espacial <- df %>% select(x, y, Cu)
coordinates(df_espacial) <- ~x + y

cv_exp <- krige.cv(Cu ~ 1, df_espacial, model = fitted_exp, nfold = nrow(df_espacial))
cv_sph <- krige.cv(Cu ~ 1, df_espacial, model = fitted_sph, nfold = nrow(df_espacial))

calcular_metricas <- function(cv_result) {
  errores <- cv_result$residual  
  ECM <- mean(errores^2)
  RECM <- sqrt(ECM)
  EAM <- mean(abs(errores))
  return(list(ECM = ECM, RECM = RECM, EAM = EAM))
}

metricas_exp <- calcular_metricas(cv_exp) 
metricas_sph <- calcular_metricas(cv_sph) 

cat("Modelo Exponencial - ECM:", metricas_exp$ECM, ", RECM:", metricas_exp$RECM, ", EAM:", metricas_exp$EAM, "\n")
## Modelo Exponencial - ECM: 391.8154 , RECM: 19.79433 , EAM: 11.63021
cat("Modelo Esférico - ECM:", metricas_sph$ECM, ", RECM:", metricas_sph$RECM, ", EAM:", metricas_sph$EAM, "\n")
## Modelo Esférico - ECM: 405.2641 , RECM: 20.13117 , EAM: 12.04876