1) Introducción

En este libro se van a utilizar dos tecnicas de interpolación espacial para obtener una superficie continua d la disponibilidad de nitrogeno en el suelo del departamento del Magdalena a una profundidad de 15 a 30 cm a artir de muestras obtenidas de SoilGrids 2.0.

Las tecnicas de interpolación espacial son: * La interpolación por distancia ponderada inversa (IDW) * Kriging Ordinario (OK)

Introducción a la interpolación por distancia ponderada inversa (IDW)

Cuando se trabaja con datos espaciales, como precipitaciones medidas en estaciones metereologicas o niveles de contaminación en ciertos puntos de una ciudad, muchas veces nos enfrentamos a una realidad incómoda: no tenemos datos para todos los lugares. ¿Qué hacemos entonces si queremos conocer un valor en una zona donde no se tomó ninguna medición?

Una de las soluciones más usadas es la interpolación espacial. Entre los métodos más sencillos y populares se encuentra la Interpolación por Distancia Ponderada Inversa

IDW parte de una idea muy simple pero poderosa, “Lo que está cerca, se parece más que lo que está lejos.”, Esto significa que si queremos estimar un valor en un punto donde no hay datos, podemos mirar a los puntos vecinos y calcular un promedio, pero no cualquier promedio: los puntos más cercanos tendrán más peso, y los más lejanos influirán menos.

Por ejemplo, si quieres saber cuánto llovió en un lugar donde no hay estación climática, puedes usar la información de estaciones cercanas, pero confiarás más en la que está a solo 5 km de distancia que en otra que está a 50 km.

*¿Como se hace el calculo?

Matemáticamente, la estimación del valor en un punto se hace así:

sus coeficientes significan: * Z(X0): Valor estimado en el punto donde no tenemos datos * Z(Xi): Es el valor conocido en el punto i, que puede ser por ejemplo la temperatura medida. * d(X0, Xi): distancia entre el punto desconocido y el punto i. *p: Es el parametro de potecia, que decide cuanto peso se le da a los puuntos cercanos. n: Es el número de puntos usados para estimar.

  • ¿Que hay que tener en cuenta? IDW depende mucho del parametro de potencia p, por lo que si lo hacemos muy alto, solo importaaran los puntos mas cercanos, lo que puede hacer que el mapa tenga picos bruscos, en cambio si lo hacemos muy bajo, todos los puntos influyen casi de la misma manera, lo que puede suavizar demasiado los datos. Tampoco detecta patrones especiales complejos, a diferencia de kriging ordinari, IDW no toma en cuenta correlaciones espaciales ni variablidad. Simplemente promedia según la distancia. Idw puede ser engañoso si los datos estan mal distribuidos, por lo que si tenemos muchos puntos juntos en una zona y pocos en otra, los resultados pueden ser sesgados.

Introducción a Ordinari kriging (OK)

Cuando nos enfrentamos a la tarea de estimar valores en ubicaciones donde no hay datos, queremos hacerlo con la mayor precisión posible. Para ello, existen muchos métodos de interpolación, pero uno de los más sofisticados, precisos y poderosos es el Kriging Ordinario (Ordinary Kriging).

A diferencia de métodos más simples como el IDW, el kriging no se basa solo en la distancia entre puntos, sino que también considera cómo se comporta la variable en el espacio, es decir, la estructura espacial o correlación geográfica de los datos. En palabras simples: aprende del patrón del fenómeno que estamos midiendo.

El Kriging Ordinario es un método geoestadístico de interpolación que asume que la media del fenómeno es constante pero desconocida en el área de estudio. Tambien uiliza los valores conocidos de una variable como niveles de pH, altura o temperatura. Puede modelar como es que estos valores se relacionan en función de la distancia y dirección, usando un semivariograma.

*¿Como se calcula?

La fórmula general del kriging es una combinación lineal ponderada de los valores conocidos:

En donde sus coeficientes significan: Z(X0): Es el valor estimado en el punto sin datos. Z(Xi): Es el valor obtenido en el punto i. λi: Es el peso aignado a cada punto, calculando con base en la estructura espacial de la variable por medio del semiovariograma. ∑λi = 1: Los pesos suman 1 para que la estimación sea insesgada.

Pero a diferencia de IDW, los pesos no se asignan simplemente según la distancia, sino también según la relación espacial estadística.

*¿Que es un semivariograma?

Cuando usamos Kriging Ordinario, no solo queremos interpolar, sino hacerlo de forma óptima, basándonos en cómo se relacionan los datos en el espacio. Y aquí es donde entra el semivariograma que nos permite cuantificar y modelar la dependencia espacial entre las observaciones.

En otras palabras: el semivariograma nos dice qué tan similares son los valores en función de su separación espacial, y eso define los pesos estadísticamente óptimos que usará el kriging para estimar puntos desconocidos.

Mientras que en IDW los pesos se asignan según la distancia, en Kriging los pesos dependen de la estructura de correlación espacial modelada por el semivariograma, tambien busca el estimador lineal insesgado con varianza mínima. Para cumplir esa promesa, necesita conocer cómo varían los datos en el espacio, y eso lo determina el modelo del semivariograma, por lo que una buena elección del modelo del mismo mejora notablemente la precisión de las estimaciones del kriging y la confiabilidad de los mapas.

  • ¿Como se usa en la practica?

Se calcula a partir de los pares de puntos disponibles, midiendo qué tanto se diferencian los valores a distintas distancias.

Sobre el gráfico experimental, se ajusta una función matemática como el modelo esférico, exponencial, gaussiano, etc. para describir la variación espacial. Esto permite estimar:

Nugget (Co): variabilidad a escala muy pequeña o error de medición. Sill (Co + C): Varianza total. *Range (a): distancia hasta la cual los valores estan correlacionados.

2) Configuración

Para empezar, lo primero que debemos hacer es descargar y llamar a las siguientes librerias:

library(sp)
library(terra)
## terra 1.8.54
library(sf)
## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(stars)
## Cargando paquete requerido: abind
library(gstat)
library(automap)
library(leaflet)
library(leafem)
library(ggplot2)
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:terra':
## 
##     intersect, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(curl)
## Using libcurl 8.10.1 with Schannel

Cada una de las librerias sirve para lo siguiente en relación con el manejo espacial y raster.

3) leer los datos de entrada

Para empezar a realizar nuestro mapa de interpolación espacial debemos de empezar cargando la capa raster .tif que descargamos anteriormente la cual tiene datos espaciales del suelo y lo imprime en pantalla para verificar que fue cargado de la manera correcta.

archivo <- ("C:\\Users\\brand\\OneDrive\\Escritorio\\GB2\\P9\\MAGDALENA_INT\\nitrogen_igh_15_30.tif")
(nitrogeno <- rast(archivo))
## class       : SpatRaster 
## size        : 1074, 594, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8386750, -8238250, 995000, 1263500  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine 
## source      : nitrogen_igh_15_30.tif 
## name        : nitrogen_igh_15_30

Podemos ver que Este raster representa el contenido de nitrógeno en el suelo a 15-30 cm de profundidad, en celdas de 250 x 250 m, con una sola capa, usando una proyección global.

Una vez cargado el rasster cramos una copia de la variable de estudio que en este caso es el nitrogeno y le asignamos una nueva variable llamada nitrogeno.perc.

nitrogeno.perc <- nitrogeno

Luego de crear la copia, realizamos una reproyección raster al sistema de coordenadas geograficas WGS 84, este paso es importante para combinar ste raster con otros datos espaciales que están en WGS84 como shapefiles de municipios, capas de puntos GPS, etc. Inicialmente, el raster estaba en la proyección Interrupted Goode Homolosine, adecuada para análisis globales pero poco compatible con coordenadas geográficas. La reproyección a WGS84 va a permitir integrar estos datos con otros insumos espaciales de uso local o global.

## Por medio de la función project(nitrogeno.perc, geog) reproyectamos el raster a wgs84
geog ="+proj=longlat +datum=WGS84"
(geog.nitro = project(nitrogeno.perc, geog))
## class       : SpatRaster 
## size        : 1106, 710, 1  (nrow, ncol, nlyr)
## resolution  : 0.002181281, 0.002181281  (x, y)
## extent      : -75.03631, -73.4876, 8.937717, 11.35021  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : nitrogen_igh_15_30 
## min value   :           78.13547 
## max value   :         1163.01575

La salida del comando permite evidenciar que el raster ha sido reproyectado exitosamente a coordenadas WGS84, manteniendo sus valores originales pero adaptado para trabajo con sistemas globales de coordenadas.

una vez que reproyectamos el raster procedemos a convertirlo en un spatraster del paquete terra, de clase stars, el cual es util cuando vamos a trabajar con funciones propias de esta biblioteca, integrarse de manera facil con sf y visualizarse por medio de ggplot2.

stars.nitro = st_as_stars(geog.nitro)

Con esta información transformada, procedemos a realizar un mapa a partir de la capa raster que se descargo anteriormente de soilGrids, permitiendo visualizar el contenido de nitrogeno disponible en el suelo por medio de una escala de colores personalizada que deja identificar las zonas con diferentes concentraciones, facilitando la interpretación geografica del fenomeno.

(N <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      stars.nitro,
      opacity = 0.8,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain =70:1170)
    ))

Se logró visualizar la distribución espacial del nitrógeno en el departamento del Magdalena mediante un mapa interactivo. Las zonas con menor concentración se representan en naranja, mientras que las áreas más ricas en nitrógeno destacan en verde y cian, lo cual permite identificar fácilmente patrones regionales que pueden orientar decisiones agrícolas o de manejo del suelo.

4) Muestreando el Magdalena

Por medio del siguiente comando se van a tomar 3000 muestras al azar en el departamento del magdalena con los datos reproyectados del contenido de nitrogeno en el suelo.

set.seed(123456) ## garantiza que el muestreo sea reproducible, es decir, siempre genera los mismos puntos aleatorios.

(samples <- spatSample(geog.nitro, 3000, "random", as.points=TRUE)) ## extrae 3.000 coordenadas espaciales (como puntos) distribuidas aleatoriamente dentro del área del raster.
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 3000, 1  (geometries, attributes)
##  extent      : -75.03522, -73.48869, 8.938808, 11.34912  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : nitrogen_igh_15_30
##  type        :              <num>
##  values      :              147.6
##                                NA
##                             140.3

Al extraer valores del raster de nitrógeno, se obtuvieron datos numéricos válidos en la mayoría de los puntos, aunque también se identificaron valores NA. Esto confirma que la cobertura no es completa en toda el área muestreada y será necesario filtrar los NA para asegurar análisis válidos.

## Por medio del paquete sf convertimos el objeto samples en un objeto de sf que es el formato indicado para manejar datos espaciales vectoriales en R.
(muestras <- sf::st_as_sf(samples))

Por medi de la tabla se puede observar que los 3.000 puntos muestreados fueron convertidos correctamente a formato sf, manteniendo su ubicación en coordenadas geográficas (WGS84) y su valor asociado de nitrógeno en el suelo. La extensión espacial cubre adecuadamente el departamento del magdalena y los atributos están correctamente estructurados. La presencia de valores NA en algunos puntos refleja zonas sin datos en el raster original, por lo que será importante filtrarlos antes de cualquier análisis estadístico o modelado espacial. Esta estructura ya permite integrar los datos con herramientas de análisis vectorial, visualización y exportación.

## Por medio de esta linea eliminamos todas las filas con valores NA en el objeto muestras, filtrando los datos muestreados que no tienen datos validos de nitrogeno, este resultado se guardara en el objeto n_muestras.
n_muestras <- na.omit(muestras)

Se depuró la muestra eliminando los puntos sin datos válidos de nitrógeno, lo que dejó un conjunto limpio n_muestras listo para análisis. Este filtrado garantiza que los resultados reflejen únicamente zonas con información real y evita distorsiones causadas por valores ausentes.

Ahora visualizamos las muestras.

longit <- st_coordinates(muestras)[,1]  ## Extrae la longitud de cada punto espacial del objeto sf muestras.
latit <- st_coordinates(muestras)[,2] ## Extrae la latitud de cada punto.
soc <- muestras$nitrogen_igh_15_30 ## Extrae los valores del contenido de nitrógeno

Se separaron las coordenadas geográficas latitud y longitud y los valores de nitrógeno asociados a cada punto de muestreo, preparándolos para análisis estadísticos, modelado espacial o visualización. Esta estructura plana facilita el uso de herramientas como modelos de regresión o interpolación.

id <- seq(1, 300, 1)

Se creó una secuencia de 300 identificadores numéricos únicos, útil para asignar etiquetas a un subconjunto de puntos muestreados. Esta numeración puede facilitar tareas como la selección, organización, trazabilidad o visualización de datos específicos dentro del análisis espacial.

Luego se crea un data frame que agrupe el identificador unico del 1 al 300, las coordenadas de longitud, coordenadas de latitud y los valores del contenido de nitrogeno en esos puntos.

(sitios <- data.frame(id, longit, latit, soc))

Se construyó una tabla con los primeros 300 puntos muestreados, cada uno identificado con su ID, coordenadas geográficas y valor de nitrógeno.

Ahora eliminamos las filas del dataframe cque contenga valores NA en alguna de sus columnas. el resultado va a ser una tabla depurada con solo los puntos que tienen información completa.

sitios <- na.omit(sitios)

Una vez que se confirma que el conjunto sitios no contiene valores faltantes, asegurando todos los puntos que tienen coordenadas y valores de nitrogeno completos, manteniendo la integridad del dataset para los analisis espaciales o estadisticos siguientes.

head(sitios)

Se verificó visualmente la estructura del conjunto sitios. Los primeros registros confirman que cada punto de muestreo cuenta con su ID, coordenadas geográficas y valor numérico de nitrógeno, en un formato limpio y listo para análisis. Esta revisión inicial facilita la detección temprana de errores de formato o valores inesperados.

Ahora procedemos a realizar nuestro mapa interactivo de la cantidad de nitrogeno en el suelo, usando una paleta de colores degradada y con los puntos de muestreo como marcadores en donde al dar click se va mostrar el valor del nitrogeno.

(N <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      stars.nitro,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 70:1170)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()))

El mapa integra exitosamente la visualización del contenido de nitrógeno en el suelo con la ubicación de 300 puntos muestreados. La combinación de la capa raster coloreada y los marcadores con agrupación mejora significativamente la interpretación espacial. Esta representación facilita identificar regiones con mayor densidad de muestreo, superposición con zonas críticas y variabilidad de los valores registrados, lo que apoya la toma de decisiones basadas en patrones espaciales observados.

5) Interpolación)

5.1) Creación del objeto gstat

La interpolación espacial es una técnica clave para estimar valores en lugares no muestreados (por ejemplo, contenido de nitrógeno en suelo) a partir de datos recogidos en puntos específicos. En R, el paquete gstat permite realizar este proceso mediante dos métodos principales: determinísticos (como IDW) y geostadísticos (como Kriging).

Para configurar la interpolación en R se crea un objeto gstat() especificando:

  • formula: la variable dependiente (y posibles covariables).

  • data: los puntos de calibración.

  • model: el modelo de variograma (opcional).

Según los argumentos utilizados, gstat() decide automáticamente:

  • Sin variograma → IDW.

  • Con modelo de variograma sin covariables → Ordinary Kriging.

  • Con variograma y covariables → Universal Kriging

El flujo típico de trabajo al trabajr con interpolación espacial es el siguiente: 1) Explorar datos: analizar la distribución espacial y realizar análisis descriptivo.

  1. Calcular variograma: definir la estructura semivariográfica de los datos.

  2. Ajuste: aplicar un modelo (Esférico, Exponencial, etc.) al variograma empírico.

  3. Crear gstat: definir el sistema de interpolación con los parámetros ajustados.

  4. Interpolar: aplicar predict() para calcular valores en un grid espacial.

Este enfoque facilita una toma de decisiones basada en la variabilidad geográfica, permitiendo detectar patrones y zonas con incertidumbre definida.

5.2) interpolación IDW

Para poder interpolar el contenido de nitrogeno disponible en el suelo, es necesario utilizar el metodo IDW, por lo cual vamos a crear el siguiente objeto gstat, en donde vamos a especificar solo la fórmula y los datos:

n1 <- gstat(formula = nitrogen_igh_15_30 ~ 1, data = n_muestras)

Se creó el objeto gstat con los datos de nitrógeno y sin definir un modelo de variograma, lo cual configura automáticamente una interpolación espacial del tipo IDW. Este método asigna mayor peso a los puntos cercanos, permitiendo generar una superficie continua que refleja la variación espacial del nitrógeno con una suposición simple y sin necesidad de modelar estructura espacial.

rrr = aggregate(geog.nitro, 2)

Se redujo la resolución espacial del raster de nitrógeno a la mitad, mediante agregación por bloques de 2x2 celdas. Esta operación suaviza los datos y mejora el rendimiento computacional para tareas como la interpolación o visualización, manteniendo patrones generales sin perder excesiva información.

rrr
## class       : SpatRaster 
## size        : 553, 355, 1  (nrow, ncol, nlyr)
## resolution  : 0.004362561, 0.004362561  (x, y)
## extent      : -75.03631, -73.4876, 8.937717, 11.35021  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : nitrogen_igh_15_30 
## min value   :           80.32337 
## max value   :         1082.81303

El objeto rrr es una versión agregada del raster original de nitrógeno, con menor resolución espacial. Esta transformación simplifica el raster, reduce el peso computacional y es útil para pruebas rápidas de interpolación o visualizaciones más generales.

Su salida representa una versión suavizada y más ligera del raster original de nitrógeno, obtenida al agrupar bloques de 2×2 celdas. Aunque mantiene la misma cobertura geográfica y sistema de coordenadas, ahora cuenta con menos celdas y valores promedio, lo que facilita análisis más rápidos y reduce la carga computacional, sin perder la estructura espacial general de la variable.

## este codigo asigna el valor 1 a todas las celdas del raster rrr, sobreescribiendo completamente los valores anteriores del nitrogeno.
values(rrr) <-1

Se sobrescribieron todos los valores del raster rrr con 1, generando una capa completamente homogénea. Esto es útil para crear máscaras, delimitar áreas de predicción, o establecer una rejilla base sobre la cual aplicar métodos como interpolación, extracción o muestreo espacial sin depender de los datos originales

names(rrr) <- "valor"

Se renombró la capa del raster rrr como “valor” para reflejar que contiene una grilla uniforme con valor constante. Este cambio mejora la claridad del objeto y facilita su uso como capa base o de soporte para interpolaciones, predicciones o combinaciones espaciales posteriores.

Actualmente, rrr es un raster auxiliar con valor constante 1 en toda su extensión y una capa renombrada como “valor”. Esta grilla homogénea es útil como base de predicción en la interpolacion IDW o, ya que define el espacio geográfico y resolución donde se proyectarán los valores estimados.

rrr
## class       : SpatRaster 
## size        : 553, 355, 1  (nrow, ncol, nlyr)
## resolution  : 0.004362561, 0.004362561  (x, y)
## extent      : -75.03631, -73.4876, 8.937717, 11.35021  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : valor 
## min value   :     1 
## max value   :     1

El raster rrr ha sido transformado en una grilla base homogénea, donde todas las celdas tienen el valor 1 y cubren el área de estudio del Magdalena. Este tipo de raster es comúnmente utilizado como rejilla de predicción para interpolaciones espaciales, ya que define dónde se estimarán los valores a partir de los datos de muestreo. Su resolución reducida optimiza el tiempo de cómputo sin perder la cobertura geográfica.

stars.rrr = st_as_stars(rrr)

Se convirtió el raster uniforme rrr en un objeto de clase stars, con el fin de utilizarlo como grilla de predicción en funciones espaciales compatibles con este formato. Esta transformación facilita la integración del raster en mapas interactivos o en procesos de interpolación espacial con gstat.

z1 = predict(n1, stars.rrr)
## [inverse distance weighted interpolation]

Se realizó una interpolación espacial por IDW para estimar el contenido de nitrógeno en toda el área de estudio. Utilizando los puntos muestreados como entrada (n1) y una rejilla homogénea como soporte (stars.rrr), se generó una superficie continua (z1) que representa la variación espacial del nitrógeno en el Magdalena. Esta superficie podrá visualizarse o analizarse posteriormente para apoyar decisiones territoriales o agronómicas.

z1 = z1["var1.pred",,]
names(z1) = "nitrogeno"

Se refinó el resultado de la interpolación IDW extrayendo solo la capa de predicción y renombrándola como “nitrogeno”. Esto facilita su interpretación, visualización y uso posterior como mapa continuo de concentración de nitrógeno en el suelo del Magdalena.

Ahora creamos una paleta de colores personalizada para este mapa.

paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), domain = 87:833, na.color = "transparent")

Una vez que tenemos la paleta de colores para interpretar el mapa y se realizo la interpolación mediante IDW, procedemos a graficarlo.

(m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 87:833)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal=paleta, values= z1$nitrogeno,
    title = "Interpolacion IDW de nitrogeno [cg/kg]"
    ))

5.3) Interpolación Ok

Para aplicar Kriging, es necesario modelar la estructura espacial de la variable de interés. El primer paso consiste en calcular el variograma empírico, el cual resume cómo varía el contenido de nitrógeno según la distancia entre puntos. Este variograma luego se ajusta con un modelo teórico que será la base para realizar predicciones más realistas mediante Kriging.

v_emp_ok = variogram(nitrogen_igh_15_30 ~ 1, data=n_muestras)

Se calculó el variograma empírico de la variable nitrógeno a partir de los puntos de muestreo disponibles. Este variograma describe cómo cambia la similitud entre valores con la distancia, y es un paso esencial para ajustar un modelo teórico de variograma que será utilizado en la interpolación por Kriging. El siguiente paso será visualizar este resultado y ajustar un modelo como el esférico, exponencial o gaussiano

plot(v_emp_ok)

El variograma empírico evidencia una relación espacial clara: a medida que aumenta la distancia entre puntos, la semivarianza también crece, lo que indica que los valores cercanos son más similares entre sí. Aunque no se identifica un rango claramente definido, la tendencia ascendente justifica aplicar un modelo teórico para avanzar con la interpolación por kriging.

v_mod_ok = autofitVariogram(nitrogen_igh_15_30 ~ 1, as(n_muestras, "Spatial"))

Se utilizó la función autofitVariogram() para ajustar automáticamente un modelo teórico al variograma empírico, tomando como variable dependiente la concentración de nitrógeno sin covariables. Esta herramienta selecciona el modelo más adecuado (entre opciones como esférico, exponencial o gaussiano) y estima de forma óptima los parámetros nugget, sill y range. El resultado proporciona una base sólida para aplicar kriging con un modelo que representa adecuadamente la estructura espacial de los datos.

plot(v_mod_ok)

El gráfico muestra un ajuste del modelo teórico tipo Ste (modelo de tipo “stable”) al variograma empírico, utilizando la función autofitVariogram(). El modelo presenta un nugget de 0, un sill de 15,245 y un range de 168, lo que indica que no hay variación aleatoria a distancias muy cortas, y que la correlación espacial se mantiene hasta aproximadamente 168 unidades de distancia. La curva ajustada sigue de manera coherente la tendencia de los puntos empíricos, lo que confirma que la estructura espacial ha sido bien modelada y está lista para ser usada en la interpolación por kriging.

v_mod_ok$var_model

El modelo ajustado consta de dos componentes: un nugget igual a 0, lo que indica que no se detecta variabilidad aleatoria no explicada a distancias muy pequeñas, y un componente Ste (stable) con un sill de 15,245.37, un range de 168.45 y un kappa de 0.2, que define la suavidad de la curva. Estos parámetros reflejan una estructura espacial continua y bien definida en los datos de nitrógeno, lo cual es favorable para aplicar kriging con un nivel de confianza aceptable en las predicciones dentro del área de estudio.

g2 = gstat(formula = nitrogen_igh_15_30 ~ 1, model = v_mod_ok$var_model, data = n_muestras)
z2= predict(g2, stars.rrr)
## [using ordinary kriging]

Se implementó el kriging ordinario con la función gstat(), utilizando el modelo de variograma previamente ajustado. La predicción se aplicó sobre la grilla de referencia (stars.rrr), generando un nuevo objeto (z2) que representa la distribución espacial estimada del nitrógeno en el área de estudio. Al incorporar la estructura de dependencia espacial identificada en los datos, esta interpolación proporciona una superficie continua que refleja variaciones locales con mayor precisión que métodos determinísticos, como IDW, y está preparada para ser visualizada o analizada en detalle.

z2 = z2["var1.pred",,]
names(z2) = "nitrogeno"

Se ha extraído la capa de predicción principal (var1.pred) del resultado del kriging y se ha renombrado como “nitrogeno” para una mejor interpretación y organización. Esta etapa permite trabajar de forma más clara con los datos interpolados, ya sea para su visualización, análisis espacial o exportación a formatos como .tif. Al simplificar el nombre de la variable, se mejora la trazabilidad del flujo de trabajo en el proyecto.

g2
## data:
## var1 : formula = nitrogen_igh_15_30`~`1 ; data dim = 2169 x 1
## variograms:
##         model    psill    range kappa
## var1[1]   Nug     0.00   0.0000   0.0
## var1[2]   Ste 15245.37 168.4509   0.2

El objeto g2 contiene la especificación del modelo geoestadístico utilizado para el kriging ordinario. Incluye un componente nugget nulo y un modelo estable (Ste) con sill de 15,245.37, rango de 168.45 y kappa de 0.2, lo cual representa una estructura espacial suave y continua. Esta configuración asegura que el modelo utilizado en la predicción respeta la dependencia espacial observada, permitiendo una interpolación robusta y coherente con los datos de campo.

Una vez que se realizo la interpolación por kriging procedemos a realizar su respectivo mapa para visualizar sus datos.

(m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 87:832)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$nitrogeno, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal = paleta, values= z2$nitrogeno,
    title = "Nitrogeno disponible [cg/kg]"
    ))

6) Analisis de resultados

6.1) Evaluación cualitativa

colores <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), na.color = "transparent")
(m <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(stars.nitro, colorOptions = colores, opacity = 0.8, group="RealWorld") %>%
  addGeoRaster(z1, colorOptions = colores, opacity = 0.8, group= "IDW")  %>%
  addGeoRaster(z2, colorOptions = colores, opacity = 0.8, group= "OK")  %>%
  # Add layers controls
  addLayersControl(
    overlayGroups = c("RealWorld", "IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>% 
    addLegend("bottomright", pal = paleta, values= z1$nitrogeno,
    title = "Nitrogeno disponible [cg/kg]"
))

El visor desarrollado permite comparar de forma dinámica las tres capas principales: los valores reales de nitrógeno (RealWorld), la interpolación por IDW (IDW) y la obtenida mediante kriging (OK). Al alternar entre las capas, se observa que el método IDW tiende a generar transiciones más abruptas entre zonas de alta y baja concentración, mientras que el kriging produce una superficie más suavizada y continua, respetando la estructura espacial definida en el modelo de variograma. Esta diferencia es especialmente notoria en áreas intermedias, donde el kriging ofrece una representación más realista de las transiciones, lo que lo posiciona como una opción más robusta para modelar fenómenos naturales como la distribución del nitrógeno en suelos agrícolas.

6.2) Validación cruzada.

Se aplicó validación cruzada tipo Leave-One-Out (LOOCV) utilizando la función gstat.cv() sobre el modelo de kriging ajustado. Este procedimiento permite evaluar de manera objetiva la precisión de la interpolación, comparando los valores observados y predichos para cada punto de muestreo, excluyéndolos uno a uno del modelo. La tabla resultante proporciona la base para calcular métricas de error como el RMSE y el bias, lo cual es esencial para determinar qué método de interpolación ofrece mayor confiabilidad en el contexto del estudio.

#cv1 = gstat.cv(g1)
#cv2 = gstat.cv(g2)

Se configuró la validación cruzada Leave-One-Out para los modelos de interpolación por IDW (cv1) y Kriging Ordinario (cv2) mediante la función gstat.cv(). Esta técnica proporciona un conjunto de errores de predicción individuales, permitiendo evaluar objetivamente el desempeño de cada método en base a sus diferencias entre valores observados y predichos. Esta comparación será clave para determinar cuál de los dos métodos ofrece mayor precisión y estabilidad al modelar la distribución espacial del nitrógeno disponible en el suelo.

#cv1 = na.omit(cv1)

Se aplicó na.omit() al objeto cv1 para eliminar predicciones con valores faltantes en la validación cruzada del modelo IDW. Esta limpieza asegura que las métricas de evaluación (como RMSE o MAE) se calculen únicamente con pares válidos de valores observados y predichos, garantizando así una comparación justa y estadísticamente sólida con el modelo de Kriging.

#cv1 = st_as_sf(cv1)

Se convirtió el objeto cv1 a clase sf mediante st_as_sf(), permitiendo así el análisis espacial directo de los errores de predicción generados en la validación cruzada del modelo IDW. Esta conversión facilita la visualización de los residuales en el espacio geográfico, lo cual puede revelar patrones sistemáticos de sobre o subestimación asociados a la distribución espacial de los puntos de muestreo.

#sp::bubble(as(cv1[, "residual"], "Spatial"))

El gráfico de burbujas representa los residuos espaciales de la validación cruzada del modelo IDW. Cada punto indica la diferencia entre el valor observado y el predicho, donde los valores negativos (en fucsia) representan sobreestimaciones y los positivos (en verde) subestimaciones. La variabilidad espacial de los residuos sugiere que el método IDW tiende a cometer errores sistemáticos en ciertas zonas, especialmente en el norte del área de estudio. Esta representación visual permite identificar posibles sesgos regionales en la interpolación y refuerza la necesidad de compararla con métodos como el kriging, que integran de forma más explícita la estructura espacial.

#sqrt(sum((cv1$var1.pred - cv1$observed)^2) / nrow(cv1))
55.73338

El cálculo del RMSE para el modelo de interpolación por IDW arrojó un valor de 55.73 cg/kg, lo que indica que, en promedio, las predicciones del modelo difieren en esa magnitud respecto a los valores observados durante la validación cruzada. Este nivel de error cuantifica la precisión del método y establece una base objetiva para compararlo con el modelo de kriging. Una vez obtenido el RMSE del segundo método, será posible determinar con mayor claridad cuál ofrece un mejor ajuste a los datos reales. ______________

#cv2 = na.omit(cv2)
#cv2 = st_as_sf(cv2)
#sp::bubble(as(cv2[, "residual"], "Spatial"))

El mapa de residuos generado para el modelo de kriging muestra una distribución espacial de errores similar a la observada en el modelo IDW, con algunas diferencias clave. Aunque persisten errores tanto de sobreestimación (fucsia) como de subestimación (verde), se aprecia una mayor dispersión y menor concentración de errores extremos. Esto sugiere que el kriging, al incorporar la estructura espacial a través del variograma, logra suavizar las predicciones y reducir errores localizados. Esta representación gráfica respalda el uso del kriging como una alternativa más precisa, especialmente en regiones donde el IDW presenta acumulación de residuos elevados.

#sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))
49.86213

Tras aplicar validación cruzada Leave-One-Out, se obtuvo un RMSE de 55.73 cg/kg para el modelo IDW y 49.86 cg/kg para el modelo de kriging. Esta diferencia indica que, en promedio, el kriging comete errores de predicción ligeramente menores, lo que refleja una mayor precisión y mejor capacidad de ajuste a los datos observados. Al combinar este resultado con el análisis visual de los residuos, se concluye que el kriging ofrece un desempeño superior, especialmente al capturar patrones espaciales que el IDW no logra representar completamente.

con 3000 muestras, se observó que el método de Kriging ordinario (RMSE = 49.86) obtuvo un error menor que el método IDW (RMSE = 55.73). Esto sugiere que el Kriging presenta una mejor capacidad predictiva en la estimación espacial del nitrógeno disponible para este conjunto de datos, con un 10.5% menos error que IDW.

7) referencias

Lizarazo, I., 2023. Spatial interpolation of soil organic carbon. Available at https://rpubs.com/ials2un/soc_interp.

Bivand, R. S., Pebesma, E., & Gómez-Rubio, V. (2013). Applied spatial data analysis with R (2nd ed.). Springer. https://doi.org/10.1007/978-1-4614-7618-4

Pebesma, E. J. (2004). Multivariable geostatistics in S: The gstat package. Computers & Geosciences, 30(7), 683–691. https://doi.org/10.1016/j.cageo.2004.03.012

Pebesma, E. (2023). gstat: Spatial and spatio-temporal geostatistical modelling, prediction and simulation. R package version 2.1-1. https://cran.r-project.org/package=gstat

Hijmans, R. J. (2023). terra: Spatial data analysis. R package version 1.7-65. https://cran.r-project.org/package=terra

Pebesma, E. (2023). stars: Spatiotemporal arrays, raster and vector data cubes. R package version 0.6-2. https://r-spatial.github.io/stars/

Cheng, J., Karambelkar, B., & Xie, Y. (2023). leaflet: Create interactive web maps with the JavaScript ‘leaflet’ library. R package version 2.1.2. https://cran.r-project.org/package=leaflet

Cressie, N. (1993). Statistics for spatial data (Revised ed.). Wiley.

sessionInfo()
## R version 4.4.3 (2025-02-28 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8  LC_CTYPE=Spanish_Colombia.utf8   
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.utf8    
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] curl_6.2.2     dplyr_1.1.4    ggplot2_3.5.2  leafem_0.2.4   leaflet_2.2.2 
##  [6] automap_1.1-16 gstat_2.1-3    stars_0.6-8    abind_1.4-8    sf_1.0-21     
## [11] terra_1.8-54   sp_2.2-0      
## 
## loaded via a namespace (and not attached):
##  [1] generics_0.1.4     sass_0.4.10        class_7.3-23       KernSmooth_2.23-26
##  [5] lattice_0.22-6     digest_0.6.37      magrittr_2.0.3     evaluate_1.0.3    
##  [9] grid_4.4.3         RColorBrewer_1.1-3 fastmap_1.2.0      plyr_1.8.9        
## [13] jsonlite_2.0.0     e1071_1.7-16       reshape_0.8.10     DBI_1.2.3         
## [17] crosstalk_1.2.1    scales_1.4.0       codetools_0.2-20   jquerylib_0.1.4   
## [21] cli_3.6.5          rlang_1.1.6        units_0.8-7        intervals_0.15.5  
## [25] withr_3.0.2        base64enc_0.1-3    cachem_1.1.0       yaml_2.3.10       
## [29] FNN_1.1.4.1        raster_3.6-32      tools_4.4.3        parallel_4.4.3    
## [33] spacetime_1.3-3    png_0.1-8          vctrs_0.6.5        R6_2.6.1          
## [37] zoo_1.8-14         proxy_0.4-27       lifecycle_1.0.4    classInt_0.4-11   
## [41] htmlwidgets_1.6.4  pkgconfig_2.0.3    pillar_1.10.2      bslib_0.9.0       
## [45] gtable_0.3.6       glue_1.8.0         Rcpp_1.0.14        tidyselect_1.2.1  
## [49] tibble_3.2.1       xfun_0.52          rstudioapi_0.17.1  knitr_1.50        
## [53] farver_2.1.2       htmltools_0.5.8.1  rmarkdown_2.29     xts_0.14.1        
## [57] compiler_4.4.3