Interpolacion

paquetes <- c('terra', 'sf', 'sp', 'stars', 'gstat', 'automap', 'leaflet', 'leafem')
## install.packages(paquetes)

1. Introducción y Preparación del Entorno

La interpolación espacial es un procedimiento geoestadístico fundamental que nos permite estimar los valores de una variable continua (como las propiedades del suelo o el clima) en áreas donde no existen mediciones directas, partiendo de un conjunto de muestras puntuales. En este estudio, analizaremos el comportamiento espacial de la Capacidad de Intercambio Catiónico (CIC) en los suelos del departamento del Cauca.

Para lograr este objetivo, emplearemos un ecosistema de librerías especializadas: terra y sf para el manejo matricial y vectorial, gstat y automap para ejecutar los modelos matemáticos de interpolación, y el conjunto leaflet/leafem para la cartografía interactiva en web.

## llamamos las librerias: 
library(terra)
library(sf)
library(sp)
library(stars)
library(gstat)
library(automap)
library(RColorBrewer)
library(leaflet)
library(leafem)

list.files(path="data", pattern = "*.gpkg")

## [1] "Cauca.gpkg"          "Depto_cauca.gpkg"    "muestras_cec_2.gpkg"
## [4] "muestreo.gpkg"       "muestreo_cauca.gpkg" "munic_cau.gpkg"     
## [7] "otros_puntos.gpkg"   "puntos.gpkg"

## Muestreos generados en QGIS: 2000 de manera aleatoria 
samples <- sf::st_read("data/muestras_cec_2.gpkg")

## Reading layer `muestras_2' from data source 
##   `C:\Users\pc\OneDrive\Documentos\GB2\RSTUDIO\cuaderno7\data\muestras_cec_2.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 2000 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -77.89965 ymin: 0.950381 xmax: -75.90228 ymax: 2.947627
## Geodetic CRS:  WGS 84

## Revisar sistemas de coordenas y forzar a que todas esten el las mismas een este 
munic <- sf::st_read("data/munic_cau.gpkg")

## Reading layer `cortado' from data source 
##   `C:\Users\pc\OneDrive\Documentos\GB2\RSTUDIO\cuaderno7\data\munic_cau.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 42 features and 90 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -77.92834 ymin: 0.9580285 xmax: -75.74782 ymax: 3.328941
## Geodetic CRS:  MAGNA-SIRGAS

2. Análisis Exploratorio de Datos

Antes de ejecutar cualquier modelo espacial, recisaremos la distribución estadística de nuestra variable de interés. En este caso, la variable SAMPLE_1 corresponde a la Capacidad de Intercambio Catiónico (CIC), un indicador crucial de la fertilidad del suelo que mide la capacidad de retener y liberar iones positivos.

Procedemos a evaluar la normalidad de los datos mediante un histograma, redondeando los valores a dos cifras decimales para estandarizar los resultados.

## Revisamos el resumen de la variable samples
summary(samples)

##        id            SAMPLE_1                geom     
##  Min.   :   0.0   Min.   :  0.0   POINT        :2000  
##  1st Qu.: 499.8   1st Qu.:182.0   epsg:4326    :   0  
##  Median : 999.5   Median :248.0   +proj=long...:   0  
##  Mean   : 999.5   Mean   :247.8                       
##  3rd Qu.:1499.2   3rd Qu.:310.0                       
##  Max.   :1999.0   Max.   :571.0

## Realizaremos un histograma, nota SAMPLES_1 corresponde a CIC
hist(samples$SAMPLE_1)

## Aproximaremos a 2 cifras decimales para los valores de CIC
samples$SAMPLE_1 = round(samples$SAMPLE_1,2)

## Vizualizaremos y ajustaremos la peleta de colores de interes. 
pal <- colorNumeric(
  c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"),
  # colors depend on the count variable
  domain = samples$SAMPLE_1,
  )

## Plotearmos las muestras 

leaflet() %>%
  addPolygons(
    data = munic,
    color = "gray",
    
    # set the opacity of the outline
    opacity = 1,
    
    # set the stroke width in pixels
    weight = 1,
    
    # set the fill opacity
    fillOpacity = 0.2) %>%
 addCircleMarkers(
    data = samples,
    radius= 1.5, 
    label = ~SAMPLE_1,
    color = ~pal(SAMPLE_1),
    fillOpacity = 1,
    stroke = F
  ) %>%
  addLegend(
    data = samples,
    pal = pal,
    values = ~SAMPLE_1,
    position = "bottomleft",
    title = "CIC:",
    opacity = 0.9) %>%
  addProviderTiles("OpenStreetMap")

## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'

3. Modelamiento Espacial: Inverso de la Distancia Ponderada (IDW)

El primer método que implementaremos es el IDW (Inverse Distance Weighting). Este modelo determinístico asume que los puntos más cercanos tienen mayor influencia sobre el punto a predecir que los puntos más alejados. Configuramos la fórmula SAMPLE_1 ~ 1, lo que indica que la interpolación dependerá únicamente de los valores de la variable en sí misma, sin variables explicativas externas (covariables).

Cuadrícula de Predicción

Para que el modelo matemático pueda “dibujar” la superficie continua, necesita un lienzo en blanco. Extraemos la caja delimitadora (BBox) de nuestras muestras y creamos un objeto ráster vacío (rrr) con 370 filas y 329 columnas. Estas dimensiones definen la resolución espacial del resultado; a mayor número de celdas, más detallado (y pesado) será el mapa resultante.

## Realizaremos la interpolacion correspondiente: asignando los datos y tipo de modelo mediante la libreria gstat. (IDW)

g1 = gstat(formula = SAMPLE_1 ~ 1, data = samples)

# Nota: definir cuadricula, pasarle al modelo los datos, creacion de raster de interes desde terra, aplicar funcion de prediccion.

# Creacion de raster de interes: 
caja <- sf::st_bbox(samples)
rrr <- terra::rast(
  xmin = caja["xmin"], xmax = caja["xmax"], 
  ymin = caja["ymin"], ymax = caja["ymax"], 
  nrows = 100, ncols = 100, vals = 1, crs = "epsg:4326"
)

## Convertiremos de SpatRaster a stars: 
stars.rrr <- stars::st_as_stars(rrr)

z1 = predict(g1, stars.rrr)

## [inverse distance weighted interpolation]

## Vizualiremos los datos de la variable z1
z1

## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                 Min.  1st Qu.   Median     Mean 3rd Qu.    Max.  NA's
## var1.pred  0.8626693 209.3914 255.8014 249.0214 294.457 513.254     0
## var1.var          NA       NA       NA      NaN      NA      NA 10000
## dimension(s):
##   from  to offset    delta refsys x/y
## x    1 100  -77.9  0.01997 WGS 84 [x]
## y    1 100  2.948 -0.01997 WGS 84 [y]

# Nota: Comentamos esta parte del codigo ya que el CIC puede tomar valores al rededor de cero por lo que esta linea de codigo disminuiria la vizualizacion de las muestras.
## NA se remplace por 0
#a <- which(is.na(z1[[1]]))
#z1[[1]][a] = 0.0

names(z1) <- c("CIC", "Varianza")

paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), domain = samples$SAMPLE_1, na.color = "transparent")

paleta <- colorNumeric(
  palette = c("orange", "yellow", "cyan", "green"), 
  domain = samples$SAMPLE_1, 
  na.color = "transparent"
)

# Extraemos las coordenadas de las 2000 muestras 
sitios <- as.data.frame(sf::st_coordinates(samples))
colnames(sitios) <- c("longit", "latit")
sitios$cic <- samples$SAMPLE_1 # Asignamos a la de CIC

# Ploteo
mapa_interpolado <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1["CIC"], # Lee la capa ya calculada 
      opacity = 0.7,                
      colorOptions = colorOptions(
        palette = c("orange", "yellow", "cyan", "green"), 
        domain = range(z1$CIC, na.rm = TRUE) # Rango automático del raster
      )
    ) %>%
  # Incorporamos los marcadores 
  addMarkers(
    lng = sitios$longit, 
    lat = sitios$latit, 
    popup = paste("<b>Valor CIC:</b>", sitios$cic), 
    clusterOptions = markerClusterOptions() # Evita que se sature el mapa
  ) %>%
  addLegend(
    "bottomright", pal = paleta, values = as.vector(z1$CIC),
    title = "IDW CIC interpolation" 
  )

mapa_interpolado

4. Modelamiento Espacial: Kriging Ordinario (KO)

A diferencia del IDW, el Kriging Ordinario es un método geoestadístico estocástico. No solo considera la distancia entre los puntos, sino también el grado de autocorrelación espacial.

Para entender esta correlación, primero generamos un variograma empírico. Luego, mediante la función autofitVariogram(), ajustamos automáticamente el modelo matemático que mejor describe el comportamiento de los puntos.

• El Nugget: Representa el error de medición o la microvariabilidad a distancias menores de las que pudimos muestrear.
• Extracción del Modelo (v_mod_ok$var_model): Esta línea de código es crítica porque aísla únicamente la ecuación matemática teórica (Esférica, Exponencial, etc.) obtenida del ajuste, la cual es el insumo obligatorio para que el algoritmo Kriging pueda predecir datos.

Ejecución del Modelo Predictivo

Posteriormente, construimos el modelo predictivo, donde: * g2: Representa el objeto de modelo Kriging empaquetado (con los datos, la fórmula y el modelo de variograma ajustado). * z2: Es el resultado final de la función predict(). Contiene la proyección del modelo g2 sobre el ráster vacío,

## El siguiente codigo calcula el variograma empirico de las muestras sin covariables. 
v_emp_ok = variogram(SAMPLE_1 ~ 1, data=samples) #<- aqui va la variable de interes (SAMPLES_1 = CIC)
plot(v_emp_ok)

## autofitvariagrma permite encontrar el mejor funcion para los puntos.
v_mod_ok = automap::autofitVariogram(SAMPLE_1 ~ 1, as(samples, "Spatial"))#<- aqui va la variable de interes(SAMPLES_1 = CIC)
plot(v_mod_ok) ## Revisar nugget

# Vizualizaremos los datos del variagrama 
v_mod_ok$var_model

##   model     psill    range kappa
## 1   Nug  359.6286  0.00000   0.0
## 2   Ste 7300.9342 61.03846   0.4

## z2 =  Predicción Final
g2 = gstat(formula = SAMPLE_1 ~ 1, model = v_mod_ok$var_model, data = samples) #<- aqui va la variable de interesc(SAMPLES_1 = CIC)
z2= predict(g2, stars.rrr)

## [using ordinary kriging]

#a <- which(is.na(z2[[1]]))
#z2[[1]][a] = 0.0
names(z2) = c("CIC", "Varianza")

sitios <- as.data.frame(sf::st_coordinates(samples))
colnames(sitios) <- c("longit", "latit")

paleta_cic <- colorNumeric(
  palette = c("orange", "yellow", "cyan", "green"), 
  domain = samples$SAMPLE_1, 
  na.color = "transparent"
)


m <- leaflet() %>%
  addTiles() %>%  
  

  leafem:::addGeoRaster(
      z2[1], 
      opacity = 0.7,                
      colorOptions = colorOptions(
        palette = c("orange", "yellow", "cyan", "green"), 
        domain = range(samples$SAMPLE_1, na.rm = TRUE) 
      )
    ) %>%
  

  addMarkers(
    lng = sitios$longit, 
    lat = sitios$latit, 
    popup = paste("<b>Valor CIC:</b>", sitios$cic), 
    clusterOptions = markerClusterOptions() 
  ) %>%
  

  addLegend(
    "bottomright", 
    pal = paleta_cic, 
    values = as.vector(z2[[1]]), 
    title = "OK CIC interpolation"
  )

## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA

m

colores <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), domain = samples$SAMPLE_1, na.color = "transparent")

# Configuración de color para los mapas continuos (Rasters)
colores_cic <- colorOptions(
  palette = c("orange", "yellow", "cyan", "green"), 
  domain = range(samples$SAMPLE_1, na.rm = TRUE) # Ajuste automático al rango a CIC
)

# Paleta 
paleta_leyenda <- colorNumeric(
  palette = c("orange", "yellow", "cyan", "green"), 
  domain = samples$SAMPLE_1,
  na.color = "transparent"
)

m <- leaflet() %>%
  addTiles() %>%  
  
    addCircleMarkers(
    data = samples, 
    radius = 1.5, 
    label = ~SAMPLE_1,
    color = ~paleta_leyenda(SAMPLE_1), 
    fillOpacity = 0.9, 
    stroke = FALSE,
    group = "RealWorld" 
  ) %>%
  
  # Capa "IDW": Superficie continua generada por el primer modelo
  leafem:::addGeoRaster(
    z1["CIC"], 
    colorOptions = colores_cic, 
    opacity = 0.8, 
    group = "IDW" 
  )  %>%
  
  # Capa "OK": Superficie continua generada por el Kriging Ordinario
  leafem:::addGeoRaster(
    z2["CIC"], 
    colorOptions = colores_cic, 
    opacity = 0.8, 
    group = "OK" 
  )  %>%
  
  
  addLayersControl(
    overlayGroups = c("RealWorld", "IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE) 
  ) %>% 
  
 
  addLegend(
    "bottomright", 
    pal = paleta_leyenda, 
    values = as.vector(z1[[1]]), 
    title = "CIC interpolation [meq/100g]" 
  )


m

5. Validación Cruzada (Cross-Validation) y Rendimiento del Modelo

Para determinar cuál de los dos métodos (IDW o Kriging Ordinario) es más preciso, aplicamos la técnica de Validación Cruzada (Leave-One-Out Cross-Validation). Este método consiste en retirar temporalmente un punto muestral del conjunto de datos, usar los puntos restantes para predecir el valor en esa ubicación vacía, y luego calcular la diferencia (el residuo) entre el valor real observado y el predicho por el modelo.

El gráfico de burbujas nos permite visualizar espacialmente estos residuos, identificando zonas donde el modelo subestima o sobreestima significativamente la variable.

Para cuantificar el error de forma global, calculamos el Error Cuadrático Medio (RMSE - Root Mean Square Error) utilizando la siguiente formulación matemática:

\[RMSE = \sqrt{\frac{\sum_{i=1}^{n} (Predicho_i - Observado_i)^2}{n}}\]

6. Conclusiones: Evaluación Cuantitativa de los Modelos

Tras ejecutar la validación cruzada para ambos métodos de interpolación, procedimos a evaluar su rendimiento global mediante el cálculo del Error Cuadrático Medio (RMSE). Los resultados arrojaron un RMSE de 44.83 para el modelo determinístico (IDW) y un RMSE de 39.92 para el modelo geoestadístico (Kriging Ordinario).

Interpretación de Resultados

La reducción del error en el segundo modelo confirma estadísticamente que el Kriging Ordinario es el estimador espacial superior para modelar la Capacidad de Intercambio Catiónico (CIC) en esta región.

Oportunidades de Mejora

A pesar de que el Kriging Ordinario mejoró la precisión, un error promedio de aproximadamente 40 unidades sigue siendo considerable. Como evidenciamos en el análisis de residuales la CIC no depende exclusivamente de la distancia. Para futuros estudios, la precisión de este mapa podría incrementarse sustancialmente evolucionando hacia un modelo de Regresión Kriging (Kriging Universal), integrando covariables.