1. °| 🧩 Introducción |°

Este Notebook ilustra dos técnicas de interpolacion espacial: Distancia Inversa Ponderada (IDW, por sus siglas en ingles) y Kriging Ordinario (OK, por sus siglas en ingles). IDW es una técnica deterministica y OK es una técnica probabilistica.

Ambas técnicas son usadas para obtener una superficie continua SOC a 15 y 30 cm a partir de datos obtenidos desde SoilGrids 250m

2. °| 🛠️ Preparacion de la herramienta |°

Antes de iniciar vamos a asegurarnos que los paquetes y librerías que necesitamos están instaladas y cargadas, para darle un buen manejo a esto se recomienda limpiar la memoria antes de iniciar:

# Para limpiar la memoria:
rm(list=ls())

# Para cargar las librerías:
library(sp)
library(terra)
library(sf)
library(stars)
library(gstat)
library(leaflet)
library(leafem)
library(ggplot2)
library(dplyr)
library(curl)
library(automap)

Para realizar solicitudes HTTP más eficientes, se crea un objeto handle con “new_handle()” y se configura para usar el protocolo HTTP/2 mediante handle_setopt(h, http_version = 2). Esta configuración permite aprovechar las ventajas del protocolo HTTP/2, como una mayor velocidad de transferencia y compatibilidad con servidores modernos que requieren este estándar.

# h <- new_handle()
# handle_setopt(h, http_version = 2)

3. °| 📖 Leer los datos de entrada |°

Necesitamos leer los datos para simular la información del mundo real. Por lo tanto vamos a leer la capa SOC que descargamos desde ISRIC usando la librería “terra”

archivo <- ("C:/Users/juane/Downloads/Valentina/GB2/p4/carbono organico_valle.tif")
(soc <- rast(archivo))
## class       : SpatRaster 
## size        : 817, 814, 1  (nrow, ncol, nlyr)
## resolution  : 0.002259887, 0.002389078  (x, y)
## extent      : -77.54802, -75.70847, 3.093174, 5.045051  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : carbono organico_valle.tif 
## name        : carbono organico_valle 
## min value   :                    142 
## max value   :                   2352

Ahora convertimos estos datos a porcentaje:

soc.perc <-  soc/10

Nota: Los valores que se descarguen de SoilGrids para SOC están expresados en dg/kg. Por lo tanto, para convertirlos a la unidad estándar (g/kg): valor_g_per_kg = valor_soc_mapeado / 10.

Recordando que siempre es muy importante validar el CRS de nuestros datos:

crs(soc.perc)
## [1] "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        MEMBER[\"World Geodetic System 1984 (G2296)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

Como vemos que ya esta en el sistema de coordenadas en el que queremos trabajar (WGS84), nos es necesario hacer una transformación. De igual manera se dejara anotado como se realizaría esta transformación de ser necesario:

##geog ="+proj=longlat +datum=WGS84"
#(geog.soc = project(soc.perc, geog))

Ahora podemos convertir la capa SpatRaster en objetos “stars”:

stars.soc = st_as_stars(soc.perc)

Nota:Los objetos stars pertenecen al paquete stars, que está diseñado para representar y analizar datos espaciales tipo ráster y multidimensionales (como tiempo, profundidad, bandas, etc.). El nombre viene de “Spatiotemporal Arrays”.Se puede imaginar como una extensión moderna de un raster común, pero más flexible y orientado a operaciones multidimensionales y más integraciones con sf.

Y para ver esto de manera gráfica:

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      stars.soc,
      opacity = 0.8,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 8:130)
    ) 
## Para ver el mapa
m

4. °| 🌎 Representando el mundo |°

Vamos a obtener una muestra de aproximadamente 500 sitios a partir de los datos reales del mundo, utilizando un muestreo aleatorio de ubicaciones.

set.seed(123)

# Random sampling of 500 points
samples <- spatSample(soc.perc, 500, "random", as.points=TRUE)

Lo que podemos ver en las características es: - Es un vector espacial “SpatVector” - Cada pieza en el objeto tiene geometría tipo punto lo que significa que tienen coordenadas X/Y - Tiene 500 geometrías, o sea 500 puntos diferentes. Y un tributo asociado a cada punto - Para la característica “extent” define las coordenadas espaciales que contienen todos los puntos - Y para la característica “coord. ref.” indica el sistema de coordenadas geográfico, en este caso es lon/lat WGS84

Ahora vamos a convertir este vector espacial “SpatVector” en un objeto de “simple feature”:

muestras <- sf::st_as_sf(samples)
head(muestras)
## Simple feature collection with 6 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -77.52429 ymin: 4.054778 xmax: -75.7096 ymax: 5.012799
## Geodetic CRS:  WGS 84
##   carbono organico_valle                   geometry
## 1                  197.1 POINT (-77.17853 4.054778)
## 2                     NA  POINT (-75.7096 4.618601)
## 3                     NA POINT (-76.61808 4.580375)
## 4                     NA POINT (-76.75593 4.764334)
## 5                     NA POINT (-77.52203 4.499147)
## 6                     NA POINT (-77.52429 5.012799)

Podemos eliminar las filas con valores NA o faltantes con el siguiente código:

nmuestras <- na.omit(muestras)

Esto evitara errores al aplicar algún modelo o calculo que no admita valores NA, ademas permite trabajar solo con observaciones completas.

Ahora, para visualizar las muestras:

longit <- st_coordinates(muestras)[,1]
latit <- st_coordinates(muestras)[,2]
soc1 <- muestras$`carbono organico_valle`

id <- seq(1,500,1) 

sitios <- data.frame(id, longit, latit, soc1)
head(sitios)
##   id    longit    latit  soc1
## 1  1 -77.17853 4.054778 197.1
## 2  2 -75.70960 4.618601    NA
## 3  3 -76.61808 4.580375    NA
## 4  4 -76.75593 4.764334    NA
## 5  5 -77.52203 4.499147    NA
## 6  6 -77.52429 5.012799    NA

En este caso al igual que antes podemos remover los valores NA:

sitios <- na.omit(sitios)

head(sitios)
##    id    longit    latit  soc1
## 1   1 -77.17853 4.054778 197.1
## 7   7 -77.06102 3.457509  37.0
## 9   9 -76.65650 3.352389  60.2
## 15 15 -76.84859 4.028498 180.9
## 19 19 -76.16158 4.155119  29.2
## 21 21 -76.01017 3.875597  50.5

Ahora veamos las muestras de manera gráfica:

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

# Print the map
m1

Nota: Este mapa muestra la distribución espacial del carbono orgánico del suelo en el departamento de interés que en nuestro caso es Valle de Cauca, para esto: 1⃣- Se utilizó un objeto tipo “star” (“stars.soc”) que contiene los valores de SOC como capa ráster espacial. 2⃣- Se aplicó una paleta de colores personalizada.Esto con la finalidad de visualizar fácilmente gradientes y transiciones de zonas más pobres a más ricas en materia orgánica. 3⃣- Además, se superpusieron 500 puntos de muestreo reales (“sitios”) con marcadores. Cada marcador contiene un popup con el valor de SOC extraído en ese sitio. Se agrupan dinámicamente con “markerClusterOptions()” para mantener el mapa limpio. 4⃣- El rango de valores (“domain = 8:130”) establece los límites visuales del ráster.

En resumen: Muy útil para estudios de fertilidad, salud del suelo o modelamiento ecológico.

5. °| 🔄 Interpolación |°

5.1. Creación del objeto gstat 📝

Para realizar una interpolacion espacial, primero debemos crear un objeto “gstat”, utilizando una función con el mismo nombre. Un objeto de gstat tiene toda la información para realizar una interpolación espacial, siendo: - La definición del modelo - Los datos de calibración

Ademas, según los argumentos, la función gstat “entiende” que tipo de modelo de interpolacion vamos a usar: - Si no se tiene variograma del modelo → IDW - Si se tiene variograma, sin covariantes → OK

Para esto vamos a usar tres parámetros en la función gstat: - Formula: la predicción “formula” especifica las variables dependientes e independientes (conocidas también como covariantes) - Datos/Data: Los datos de calibración o datos de entrenamiento - Modelo: El modelo del variograma.

5.2. Interpolación IDW 🗺️

Para interpolar nuestro SOC usando IDW podemos usar la función “predict()” para estimar los valores. Esta función acepta: - Un raster → Un objeto stars, como un DEM - Un modelo → Un objeto gstat, como g1

El raster nos sirve para dos funciones: - Especificar la ubicación donde queremos hacer las predicciones (en todos los métodos) - Especificar las covariantes (solo en Kriging Universal (UK, por sus siglas en ingles))

De este modo vamos a crear un objeto raster con valores de celda iguales a 1:

rrr = aggregate(soc.perc, 4)

Definimos los valores nuevos:

values(rrr) <- 1

Definimos los nombres nuevos:

names(rrr) <- "valor"

Validamos que tenemos en “rrr” ahora:

rrr
## class       : SpatRaster 
## size        : 205, 204, 1  (nrow, ncol, nlyr)
## resolution  : 0.009039548, 0.009556314  (x, y)
## extent      : -77.54802, -75.70395, 3.086007, 5.045051  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : valor 
## min value   :     1 
## max value   :     1

Con esto creamos un objeto stars:

stars.rrr = st_as_stars(rrr)

Ahora, podemos crear el siguiente objeto gstat, especificando la formula y los datos:

g1 <- gstat(formula = carbono.organico_valle ~ 1, data = nmuestras)

De esta manera, con el siguiente código interpolamos el SOC de acuerdo con el modelo definido en g1 y la plantilla del raster que definimos:

## Interpolacion por distancia inversa - IDW
z1 = predict(g1, stars.rrr)
## [inverse distance weighted interpolation]
z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean  3rd Qu.     Max.  NA's
## var1.pred  16.00197 43.65463 53.83906 67.24908 86.96741 195.0015     0
## var1.var         NA       NA       NA      NaN       NA       NA 41820
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 204 -77.55   0.00904 WGS 84 [x]
## y    1 205  5.045 -0.009556 WGS 84 [y]

Nota: La interpolación por distancia inversa (IDW) es un método determinista que estima valores en ubicaciones que no se tengan los datos usando una media de los puntos cercanos. La idea central es que los puntos más cercanos tienen mayor influencia en el valor estimado que los más lejanos. Este método es útil cuando no contamos con un modelo de variograma (como en kriging, en el siguiente numeral se realiza con este método) y queremos hacer una estimación basada puramente en la proximidad espacial.

La función “predict()” aplica la interpolación a cada celda del raster, usando los puntos de muestreo definidos en g1 y devolviendo un nuevo raster “z1” con los valores estimados.

Vamos a extraer el primer atributo de z1 y renombrarlo:

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

Ahora vamos a empezar a graficar, necesitamos un paleta de colores:

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

Organizamos el código para ver el SOC que interpolamos:

m2 <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 11:55)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal=paleta, values= z1$soc,
    title = "IDW SOC interpolation [%]"
    )
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
# Print the map
m2

De esta manera, podemos ver un mapa interactivo que nos permite visualizar los datos que interpolamos de manera mas sencilla, se superpone los datos que interpolamos de carbono orgánico (SOC). Ademas se muestran los sitios reales de muestreo con marcadores.

5.3. Interpolación OK 🗺️

Los métodos de Kriging requieren un modelo de variograma. El variograma es una herramienta fundamental para cuantificar de manera objetiva el patrón de autocorrelación en los datos y permite asignar pesos adecuados al momento de realizar predicciones.

La función requiere dos argumentos: - “formula”: especifica la variable dependiente y las covariables, de forma similar a como se hace en la función gstat. - “data”: es la capa de puntos (en nuestro caso seria el objeto “muestras”) que contiene los valores de la variable dependiente (el SOC) y, si es el caso, las covariables, almacenadas como atributos.

Con la siguiente expresión podemos calcular el variograma empírico de “muestras”, sin las covariables:

v_emp_ok = variogram(carbono.organico_valle ~ 1, data=nmuestras)

Lo graficamos:

plot(v_emp_ok, main = "Variograma Empirico")

En la gráfica del variograma empírico se observa que la semivarianza inicia en un valor aproximado de 350–400 incluso a distancias muy cortas, lo cual corresponde al nugget, indicando que existe cierta variabilidad no explicada por la distancia (posiblemente asociada a errores de medición o variaciones a escalas más pequeñas que el muestreo). A medida que la distancia entre puntos aumenta, también lo hace la semivarianza, lo que sugiere que los datos cercanos son más similares entre sí. Esta tendencia parece estabilizarse alrededor de un valor de semivarianza entre 2500 y 2600, lo que indica el sill, y ocurre a una distancia cercana a los 80 metros, que sería el rango. Por lo que los puntos que tengan mas de 80 unidades se considerarían independientes.

Por lo tanto, parece que existe una autocorrelacion espacial con los datos.

Hay varias maneras de ajustar un modelo de variograma a un variograma empírico. En este caso usaremos la mas sencilla, ajustando automáticamente con la función “autofitVariogram” de la librería automap. De la siguiente manera:

v_mod_ok = autofitVariogram(carbono.organico_valle ~ 1, as(nmuestras, "Spatial"))

Ahora graficamos otra vez:

plot(v_mod_ok)

En este gráfico podemos ver que hay autocorrelacion espacial, por lo que igual que en el anterior, la semivarianza aumenta con la distancia, el nugget de 317 surgiere que hay una variabilidad considerable, se ve una pequeña problemática principal por la discrepancia entre el sill y el range, por lo que probablemente signifique que no se muestreo lo suficientemente lejos para captar por completo la estructura de autocorrelacion de los datos.

El objeto resultante es en realidad una lista de varios componentes, incluyendo el variograma empírico y el ajustado que vimos en los códigos pasados, el componentes “$var_model” del objeto “v_mod_ok” contiene el modelo como tal. Para verlo en lista:

v_mod_ok$var_model
##   model       psill    range kappa
## 1   Nug    316.7849    0.000   0.0
## 2   Ste 114288.3201 2589.703   0.7
  • “Nug” se refiere al nugget, es decir, la variabilidad a distancia cero.
  • “Ste” es el modelo esférico (Spherical) usado para ajustar el variograma.

| psill | Sill parcial o contribución de varianza de cada componente del modelo:

  • El valor para “Nug” es 316.78, lo cual representa el nugget.
  • El valor para “Ste” es 114288.32, y representa la varianza estructurada que se estabiliza con la distancia.

| range | Rango: indica a partir de qué distancia las observaciones se vuelven independientes.

  • Para “Nug” es 0 porque no tiene estructura espacial.
  • Para “Ste” es 2589.70, lo que significa que los puntos separados por más de ~2590 unidades ya no están correlacionados espacialmente.

| kappa | Parámetro adicional usado por algunos modelos como “Matérn” o “Ste” (aunque Ste no siempre requiere kappa, algunos paquetes lo incluyen). En este caso, “Ste” tiene un kappa de 0.7.

Ahora, el modelo del variograma puede ser pasado a una función “gstat” y podemos realizar el OK:

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

De igual manera que en IDW, vamos a extraer el primer atributo y renombrarlo:

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

Las predicciones creadas por OK las podemos visualizar de la siguiente manera:

m3 <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 11:55)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal = paleta, values= z2$soc,
    title = "OK SOC interpolation [%]"
    )
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
# Print the map
m3

De esta manera, podemos ver un mapa interactivo que nos permite visualizar los datos que interpolamos de manera mas sencilla, se superpone los datos que interpolamos de carbono orgánico (SOC). Ademas se muestran los sitios reales de muestreo con marcadores.

6. °| 📊 Evaluacion de resultados |°

6.1. Analisis cualitativo 🧐

Vamos a ver de manera conjunta las dos interpolaciones y el mundo real, para tener una visión mas clara de lo realizado:

Primero creamos un paleta de colores acorde:

colores <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), domain = 10:100, na.color = "transparent")

Unimos en una misma gráfica:

m4 <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(stars.soc, opacity = 0.8, colorOptions = colores, 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$soc,
    title = "Soil organic carbon [%]"
)
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
# Print the map
m4

En este mapa se muestran los resultados de dos métodos de interpolación espacial aplicados sobre la variable carbono orgánico del suelo (%) en el Valle del Cauca

La leyenda nos permite agregar y quitar los valores de las técnicas que queramos ver, haciendo mas fácil comparar y revisar las diferencias entre los dos métodos usados (IDW y OK)

Método Descripción general Ventajas Desventajas Observaciones visuales
RealWorld (1° imagen) Muestra los valores reales observados en puntos de muestreo. Representa fielmente los datos medidos. No permite conocer valores entre puntos. Se visualiza una cobertura dispersa y discontinua. La información espacial está limitada a donde se tomaron muestras.
IDW (2° imagen) Interpolación basada en la cercanía a los puntos muestreados. Valores más cercanos tienen mayor influencia. Método simple, rápido, sin necesidad de un modelo de variograma. No tiene en cuenta la estructura espacial; puede generar artefactos o bordes duros. Se observan zonas con transiciones abruptas y patrones algo artificiales, especialmente en áreas con pocos datos.
OK - Kriging ordinario (3° imagen) Método geoestadístico que utiliza un modelo de variograma para ponderar las distancias y correlaciones espaciales. Ofrece estimaciones más suaves y confiables al modelar la dependencia espacial. Requiere ajuste del variograma; es más complejo y computacionalmente intensivo. El mapa muestra una superficie más continua y suavizada, con patrones espaciales que parecen más naturales y generalizados.

En la visualización se comparan tres representaciones del carbono orgánico del suelo: los puntos reales de muestreo, la interpolación por distancia inversa (IDW) y el kriging ordinario (OK). Mientras que IDW ofrece una interpolación rápida pero sensible a la distribución de los puntos, el método de kriging produce una superficie más coherente al incorporar la estructura espacial mediante el variograma. Esta comparación evidencia cómo la elección del método influye en la representación espacial del fenómeno.

6.2. Validacion cruzada 🧠

Hemos estimado superficies climáticas utilizando dos métodos distintos de interpolación: IDW y OK. Aunque es útil comparar los resultados de forma visual, también es necesario contar con un criterio objetivo para evaluar la precisión de las interpolaciones. Esto nos permite seleccionar de manera fundamentada el método más preciso entre las opciones disponibles.

Una forma común de validación es la Validación Cruzada Leave-One-Out Cross Validation, que consiste en: 1. Excluir un punto del conjunto de datos de calibración. 2. Realizar una predicción para ese punto usando los datos restantes. 3. Repetir el proceso para todos los puntos.

En R, podemos realizar esta validación usando la función gstat.cv(), que trabaja directamente con un objeto gstat:

#### Este código se corre desde la consola y linea por linea
##-1 cv1 = gstat.cv(g1)
#-2  cv2 = gstat.cv(g2)
cv1 = gstat.cv(g1)
cv2 = gstat.cv(g2)

El objeto resultante de “gstat.cv” incluye varias columnas que nos permiten analizar detalladamente el desempeño de los métodos de interpolación aplicados. Se describen los principales atributos:

Atributo Descripción
var1.pred Valor predicho.
var1.var Varianza (solo disponible en métodos de Kriging).
observed Valor observado.
residual Diferencia entre el valor observado y el predicho: observed - var1.pred.
zscore Puntaje Z, útil para evaluar si los errores de predicción están dentro del rango esperado. Solo en Kriging.
fold Identificador del pliegue de validación cruzada. 
cv1 = na.omit(cv1)

cv1
## class       : SpatialPointsDataFrame 
## features    : 229 
## extent      : -77.42486, -75.77966, 3.101536, 4.953072  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 6
## names       :        var1.pred, var1.var, observed,          residual, zscore, fold 
## min values  : 19.2026712256364,       NA,     15.5, -83.3466075240574,     NA,    1 
## max values  : 163.767048489252,       NA,    197.1,  80.6998540681505,     NA,  229

Convertimos el objeto cv1:

cv1 = st_as_sf(cv1)

Ahora graficamos los residuales:

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

Se grafican los residuales, es decir, las diferencias entre los valores observados y los valores predichos en cada punto del muestreo. Los colores nos indican que:

  • Color verde: Residuo positivo, indica que el modelo subestimó el valor observado, osea, el modelo predijo menos de lo que realmente hay.
  • Color fucsia: Residuo negativo, el modelo sobreestimó el valor observado, el modelo predijo más de lo que había.

El tamaño de cada burbuja representa la magnitud absoluta del residuo, por lo tanto:

  • Burbuja grande = error alto
  • Burbuja pequeña = error bajo

Conclusión del gráfico: Hay una mezcla de errores positivos y negativos, lo cual es esperable. Existen algunos residuos grandes, como uno negativo de -83.3 y uno positivo de 80.7, son posiblemente outliners o zonas mal representadas por el modelo. El modelo parece tener un desempeño razonable en la mayoría de las áreas.

Ahora vamos a calcular los indices de precisión de las predicciones, como el Error Cuadrático Medio o RMSE por sus siglas en ingles:

# Para calcular el RMSE para la interpolacion IDW 
sqrt(sum((cv1$var1.pred - cv1$observed)^2) / nrow(cv1))
## [1] 23.2409

Repetimos el proceso para la interpolacion OK:

Convertimos el objeto “cv2”:

cv2 = st_as_sf(cv2)

Aplicamos el RMSE a los resultados del OK:

# Para calcular el RMSE para la interpolacion OK
sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))
## [1] 20.83259

Con base en los valores del Error Cuadrático Medio (RMSE) obtenidos mediante validación cruzada Leave-One-Out, se observa que la técnica de interpolación Kriging Ordinario (OK) presenta mayor precisión que la ponderación por distancia inversa (IDW). El RMSE para IDW fue de 23.24, mientras que para OK fue de 20.83, lo que indica que las predicciones realizadas con Kriging fueron, en promedio, más cercanas a los valores observados. Dado que un RMSE más bajo indica mejor desempeño del modelo, se puede concluir que Kriging Ordinario es el método de interpolación más preciso en este caso. Este resultado es coherente con la ventaja conocida del Kriging, que considera la autocorrelación espacial, mejorando así la calidad de las predicciones.

°💻 ENTORNO INFORMATICO°

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

sessionInfo()
## R version 4.5.1 (2025-06-13 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## 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] automap_1.1-20 curl_6.4.0     dplyr_1.1.4    ggplot2_3.5.2  leafem_0.2.4  
##  [6] leaflet_2.2.2  gstat_2.1-4    stars_0.6-8    abind_1.4-8    sf_1.0-21     
## [11] terra_1.8-60   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-7     digest_0.6.37      magrittr_2.0.3     RColorBrewer_1.1-3
##  [9] evaluate_1.0.4     grid_4.5.1         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] scales_1.4.0       crosstalk_1.2.1    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        tools_4.5.1        raster_3.6-32      parallel_4.5.1    
## [33] spacetime_1.3-3    vctrs_0.6.5        R6_2.6.1           png_0.1-8         
## [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.11.0      bslib_0.9.0       
## [45] gtable_0.3.6       glue_1.8.0         Rcpp_1.1.0         tidyselect_1.2.1  
## [49] tibble_3.3.0       xfun_0.52          knitr_1.50         farver_2.1.2      
## [53] htmltools_0.5.8.1  rmarkdown_2.29     xts_0.14.1         compiler_4.5.1