1. Introducción

Este cuaderno presenta y explica dos métodos de interpolación espacial: Distancia Inversa Ponderada (IDW) y Kriging Ordinario (OK). La técnica IDW es un enfoque determinista que asigna valores a puntos desconocidos en función de la proximidad a puntos de muestreo conocidos, asumiendo que los más cercanos tienen mayor influencia. Por otro lado, OK es un método geoestadístico basado en modelos probabilísticos, que no solo considera la distancia entre puntos, sino también la estructura espacial de los datos a través de una función de variograma. Ambas técnicas se aplican en este análisis para generar una superficie continua del carbono orgánico del suelo (SOC) en la profundidad de 15-30 cm, utilizando datos de SoilGrids con una resolución espacial de 250 metros.

Este cuaderno se construye como forma de presentación del segundo informe del curso de Geomática Básica del período 2024-II, por lo cual, además de la interpretación de los datos de interpolación espacial, también se busca dar una explicación detallada sobre el código utilizado con el fin de demostrar los conocimientos adquiridos en el curso sobre R en el curso de geomática.

2. Instalación de programas

En primer lugar, limpiaremos la memoria, para optimizar el rendimiento del código

#rm(list=ls())

Para manejar datos geoespaciales, utilizaremos las bibliotecas terra para el procesamiento de datos ráster y sf para el manejo de datos vectoriales. Para realizar la interpolación espacial, emplearemos las bibliotecas gstat y automap, que facilitan la estimación y ajuste de modelos geoestadísticos. En cuanto a la visualización, recurriremos a ggplot2 para gráficos estáticos, leaflet para mapas interactivos y tmap para la representación cartográfica avanzada. Antes de comenzar, es fundamental asegurarse de que todas estas bibliotecas estén instaladas. Se recomienda hacerlo directamente desde la consola de R, en lugar de ejecutarlo dentro del cuaderno (notebook), para evitar posibles conflictos en la instalación y asegurar un mejor rendimiento del código.

En caso de que no se tenga ningún paquete instalado se recomienda usar el siguiente código: paquetes = c(‘terra’, ‘sf’, ‘sp’, ‘stars’, ‘gstat’, ‘automap’, ‘leaflet’, ‘leafem’) install.packages(paquetes)

En esta caso como ya tenemos la mayoría de los paquetes instalados, solo procederemos a instalar los paquetes que nos hagan falta de la siguiente manera: paquetes = c(‘gstat’, ‘automap’) install.packages(paquetes)

Ahora sí podemos importar todas las bibliotecas para asegurarnos que estén en funcionamiento

library(terra)
## terra 1.8.10
library(ggplot2)
library(sf)
## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE
library(sp)
library(stars)
## Cargando paquete requerido: abind
library(gstat)
library(automap)
library(RColorBrewer)
library(leaflet)
library(leafem)

3. Leer datos de entrada

Como se indicó en la introducción, nuestras muestras tienen valores de SOC a 15-30 cm de profundidad. Se obtuvieron de SoilGrids 250 m utilizando un cuaderno anterior donde se guardaron los datos en un archivo geopack.

En primer lugar, comprobemos que los datos están disponibles para su uso, para esto usaremos la función “list.files” que nos permitirá leer los archivos de una carpeta específica (En nuestro caso la carpeta se llama “DATOS”) que tengan la extensión “gpkg” puesto que los datos que obtuvimos en el otro cuaderno se guardaron en un geopack. Cabe aclarar que el Notebook y los datos deben estar en la misma carpeta de archivos, ya que de lo contrario la función no podrá encontrarlos.

list.files(path="DATOS", pattern = "*.gpkg")
## [1] "aguas_caquetá.gpkg"    "caquetá_box.gpkg"      "ciudades_Caquetá.gpkg"
## [4] "dpto.shp.gpkg"         "dpto_caquetá.gpkg"     "mun_caquetá.gpkg"     
## [7] "rios_Caquetá.gpkg"     "roads_caquetá.gpkg"    "soc_caq.gpkg"

El archivo donde se guardaron los datos es el de “soc_caq.gpkg”. Ahora se buscará la ruta específica de este archivo y se leerá con la función sf para posteriormente guardarse en el vector “samples”

samples <- sf::st_read("C:/Users/Karol Ruiz/Documents/GEOMATICA 2024-2/PROYECTO4/DATOS/soc_caq.gpkg")
## Reading layer `SOC_reducido' from data source 
##   `C:\Users\Karol Ruiz\Documents\GEOMATICA 2024-2\PROYECTO4\DATOS\soc_caq.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 2000 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -76.30642 ymin: -0.7034065 xmax: -71.22078 ymax: 2.9375
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

En este archivo podemos ver las coordenadas de latitud y altitud del marco del departamento en formato CRS WGS84 y la cantidad de muestras tomadas en el departamento.

Adicionalmente traeremos los datos del departamento del Caquetá que es nuestro departamento de estudio.

munic <- sf::st_read("C:/Users/Karol Ruiz/Documents/GEOMATICA 2024-2/PROYECTO4/DATOS/mun_caquetá.gpkg")
## Multiple layers are present in data source C:\Users\Karol Ruiz\Documents\GEOMATICA 2024-2\PROYECTO4\DATOS\mun_caquetá.gpkg, reading layer `municipios_colombia'.
## Use `st_layers' to list all layer names and their type in a data source.
## Set the `layer' argument in `st_read' to read a particular layer.
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## automatically selected the first layer in a data source containing more than
## one.
## Reading layer `municipios_colombia' from data source 
##   `C:\Users\Karol Ruiz\Documents\GEOMATICA 2024-2\PROYECTO4\DATOS\mun_caquetá.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 16 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.30622 ymin: -0.70584 xmax: -71.25385 ymax: 2.938334
## Geodetic CRS:  MAGNA-SIRGAS

Como podemos ver, el archivo de los municipios esta en formato MAGNA SIRGAS así que vamos a transformar su formato para que este en coordenadas de latitud y longitud al igual que las muestras de SOC.

(munic_wgs84 <- st_transform(munic, crs = 4326))
## Simple feature collection with 16 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.30622 ymin: -0.70584 xmax: -71.25385 ymax: 2.938334
## Geodetic CRS:  WGS 84
## First 10 features:
##    dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr             mpio_cnmbr
## 1          18        001      18001    CAQUETÁ              FLORENCIA
## 2          18        029      18029    CAQUETÁ                ALBANIA
## 3          18        094      18094    CAQUETÁ BELÉN DE LOS ANDAQUÍES
## 4          18        150      18150    CAQUETÁ   CARTAGENA DEL CHAIRÁ
## 5          18        205      18205    CAQUETÁ                CURILLO
## 6          18        247      18247    CAQUETÁ            EL DONCELLO
## 7          18        256      18256    CAQUETÁ              EL PAUJÍL
## 8          18        410      18410    CAQUETÁ           LA MONTAÑITA
## 9          18        460      18460    CAQUETÁ                  MILÁN
## 10         18        479      18479    CAQUETÁ                MORELIA
##                              geom
## 1  MULTIPOLYGON (((-75.47547 1...
## 2  MULTIPOLYGON (((-75.89531 1...
## 3  MULTIPOLYGON (((-75.78705 1...
## 4  MULTIPOLYGON (((-74.78722 1...
## 5  MULTIPOLYGON (((-76.09066 1...
## 6  MULTIPOLYGON (((-75.36167 2...
## 7  MULTIPOLYGON (((-75.36638 2...
## 8  MULTIPOLYGON (((-75.47712 1...
## 9  MULTIPOLYGON (((-75.32574 1...
## 10 MULTIPOLYGON (((-75.77185 1...

4. Exploración de datos de entrada

Ahora leeremos un resumen de los datos descargados con la función summary

summary(samples)
##  soc_igh_15_30                geom     
##  Min.   :  91.32   POINT        :2000  
##  1st Qu.: 129.22   epsg:NA      :   0  
##  Median : 176.37   +proj=long...:   0  
##  Mean   : 244.29                       
##  3rd Qu.: 308.34                       
##  Max.   :1193.53                       
##  NA's   :47

Como podemos observar, en la base de datos encontramos un total de 2000 muestras siendo el dato promedio de SOC en todas las muestras 244.29, el valor mínimo de 91.32 y el valor máximo de 1193.53. También podemos observar que hay 47 valores faltantes que probablemente no tiene muestras de datos.

Para tener más clara la información vamos a visualizar los datos en un histograma que nos permitirá interpretar los datos de forma gráfica

hist(samples$soc)

En el histograma se corroboran los datos presentados anteriormente, sobre los valores máximos y mínimos de muestra, los cuartiles, la media y la mediana.

Ahora redondearemos los datos a 2 cifras para su mejor manejo. Esto se hará con la siguiente función:

samples$soc = round(samples$soc,2)

Ahora visualizaremos los datos en un mapa interactivo con ayuda de la función leaflet. El primer paso es escoger una paleta de colores adecuada para la visualización de las muestras. En mi caso decidí escoger los colores ““#66CD00” (Verde oscuro), “#BCEE68” (Verde claro), “#F9D423” (Amarillo), “#FC913A” (Naranja) y “#FF4E50” (Rojo), en ese orden respectivamente.

pal <- colorNumeric(
  c("#66CD00", "#BCEE68", "#F9D423", "#FC913A", "#FF4E50"),
  domain = samples$soc,
  )

Ahora si usaremos la función leaflet para generar el mapa, en primer lugar, vamos a agregar los polígonos de la capa “munic” que es la que contiene los municipios, con algunas especificaciones como por ejemplo tener un contorno gris, una opacidad total para los bordes y un grosor de los bordes de 1 píxel. Después se agregarán los puntos de las coordenadas de la capa “samples” que es la que contiene los datos de las muestras SOC. Estos puntos los vamos a específicar con un tamaño de radio de 1.5 y con los colores respectivos de la paleta que hemos construido anteriormente, les daremos también una opacidad total del 100% e indicaremos que no tengan borde alrededor. Finalmente agregaremos una leyenda que nos permita interpretar los colores de la muestra y la ubicaremos en la esquina inferior izquierda.

leaflet() %>%
  addPolygons(
    data = munic,
    color = "gray",
    opacity = 1,
    weight = 1,
    fillOpacity = 0.2) %>%
 addCircleMarkers(
    data = samples,
    radius= 1.5, 
    label = ~soc,
    color = ~pal(soc),
    fillOpacity = 1,
    stroke = F
  ) %>%
  addLegend(
    data = samples,
    pal = pal,
    values = ~soc,
    position = "bottomleft",
    title = "SOC:",
    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'

En el mapa podemos ver un total de 2000 muestras con 47 valores NA, es decir que se trabaja plenamente con 1953 muestras, de las cuales la mayoría se encuentran en la gama de colores verde, lo cual indica que en la mayoría de muestras presentan bajos niveles de SOC en el suelo a 15-30cm de profundidad. Con estos datos procederemos a hacer las interpolaciones correspondientes.

5. Interpolación de datos

5.1. Creación del objeto gstat

Para hacer la interpolación de datos primero deberemos crear un objeto de clase gstat usando una función con el mismo nombre. En R esta función permite construir modelos de interpolación espacial almacenando toda la información necesaria para realizar predicciones en ubicaciones donde no se tienen mediciones directas. Es una herramienta clave en geoestadística, utilizada para aplicar métodos de interpolación como Kriging y Distancia Inversa Ponderada (IDW).

Su principal función es definir un modelo de interpolación con base en datos geoespaciales seleccionando el método de interpolación, ya sea Kriging ordinario (OK), Kriging universal (UK) o IDW, aunque también se utiliza para ajustar modelos de semivarianza para poder modelar la relación espacial entre puntos y predecir valores en áreas no muestreadas.

Ahora hablaremos un poco de los modelos que puede usar gstat. Por un lado, tenemos la distancia inversa ponderada (IDW) la cual se aplica cuando no se define un modelo de variograma y estima los valores en función de la distancia de los puntos conocidos. Y por otro lado, tenemos el modelo kriging ordinario (OK) que se usa cuando si se específica el modelo de variograma, este modelo supone que la media del proceso es constante en toda la región de estudio.

Para usar esta función tendremos en cuenta 3 parámetros:

  1. La fórmula → Especifica la relación entre la variable a predecir y las covariables. • Para IDW: formula = soc ~ 1 • Para Kriging con covariables: formula = soc ~ elev + slope

  2. datos → Contiene los datos de entrenamiento, que pueden estar en formato data.frame, sf o sp.

  3. modelo → Define el modelo de variograma. Recordemos que: • Si no se especifica, se aplica IDW. • Si se define, se utiliza Kriging para la interpolación.

5.2 Interpolación IDW

Para hacer la interpolación del SOC mediante el método IDW tenemos que definir el espacio gstat usando la fórmula descrita anteriormente. Sin embargo, como observamos en el mapa interactivo de la muestra, teníamos algunos valores NA, por lo cual el primer paso es eliminar esos valores para que no nos generen problemas más adelante

sum(is.na(samples))
## [1] 94
g1_data <- na.omit(samples)

Ahora si podemos definir el espacio con la fórmula para el método IDW

g1 = gstat(formula = soc ~ 1, data = g1_data, set = list(idp = 2))

Ahora que nuestro modelo de interpolación g1 está definido, podemos usar la función de predicción para estimar los valores de SOC. Esta función acepta tanto objetos ráster como objetos de modelo gstat.

El uso de un ráster en el análisis espacial cumple dos funciones principales: - Definir las ubicaciones donde se realizarán las predicciones, independientemente del método de interpolación utilizado. - Asignar valores a las covariables en el caso de aplicar kriging universal (UK) permitiendo mejorar la precisión de las estimaciones.

Para ello, utilizaremos la biblioteca terra para generar un objeto ráster cuyas celdas tengan un valor uniforme de 1. La extensión del ráster debe abarcar completamente el área de estudio para garantizar que las predicciones cubran toda la región de interés.

Nota: Puede utilizar los metadatos de cualquier ráster que cubra su departamento (por ejemplo, un modelo digital de elevación) para definir los parámetros de creación del ráster:

(rrr <- terra::rast(xmin=-76.31313, xmax=-71.21854, ymin=-0.7056429, ymax=2.939737  , nrows=815, ncols = 1139, vals = 1, crs = "epsg:4326"))
## class       : SpatRaster 
## dimensions  : 815, 1139, 1  (nrow, ncol, nlyr)
## resolution  : 0.004472862, 0.004472859  (x, y)
## extent      : -76.31313, -71.21854, -0.7056429, 2.939737  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : lyr.1 
## min value   :     1 
## max value   :     1

Vamos a reducir la resolución del ráster agregando celdas en bloques de 2x2 con el fin de ahorrar un poco de memoria.

rrr_reducido <- terra::aggregate(rrr, fact = 2, fun = mean)
terra::res(rrr_reducido)
## [1] 0.008945724 0.008945718

Ahora, vamos a convertir este SpatRaster reducido en un objeto stars con el siguiente código:

stars.rrr <- stars::st_as_stars(rrr_reducido)

La razón por la cual convertimos el objeto es porque los objetos stars nos permiten manejar datos ráster y vectoriales en un solo objeto. Además, nos facilita la manipulación de datos especiales con sf y nos permite manejar múltiples dimensiones. Ahora si podemos proceder a realizzar la respectiva predicción. Esta se hará con la función predict() que realiza una interpolación espacial sobre una cuadrícula (stars.rrr) utilizando un modelo previamente definido (g1)

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

En el vector z1 nos quedará guardada la predicción o las estimaciones del modelo IDW con los datos de SOC (Cantidad de carbono orgánico en el suelo) en el departamento del Caquetá, esto según los datos conocidos.

Para entender un poco mejor, vamos a ver que tiene este objeto específicamente

z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean 3rd Qu.     Max.   NA's
## var1.pred  92.06756 149.7301 195.6765 244.8358 323.215 1181.687    977
## var1.var         NA       NA       NA      NaN      NA       NA 232560
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 570 -76.31  0.008946 WGS 84 [x]
## y    1 408   2.94 -0.008946 WGS 84 [y]

z1 tiene valores mínimos de 92,06756 y valores máximos de 203,0913. También podemos ver que su promedio esta sobre 244.8358. También podemos evidenciar que el modelo no ofrece una varianza conocida.

Ahora vamos a crear un subconjunto solo del primer atributo que son los valores estimados de SOC asignándole a los valores de NA un valor de 0 y renombrándolo.

a <- which(is.na(z1[[1]]))
z1[[1]][a] = 0.0
z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##            Min. 1st Qu.   Median     Mean 3rd Qu.     Max.   NA's
## var1.pred     0 149.236 195.2021 243.8072 322.672 1181.687      0
## var1.var     NA      NA       NA      NaN      NA       NA 232560
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 570 -76.31  0.008946 WGS 84 [x]
## y    1 408   2.94 -0.008946 WGS 84 [y]
names(z1) = "soc"

Ahora vamos a elegir una paleta de colores adecuada para visualizar la interpolación de los datos.

paleta <- colorNumeric(palette =  c("cyan", "green", "yellow", "orange", "red"), domain = 0:110, na.color = "transparent")
color = ~paleta(soc)

Y ahora visualizaremos la interpolación en un mapa interactivo con la función leaflet el cual guardaremos en el vector “m”

m <- leaflet() %>%
  addTiles() %>%
  leafem:::addGeoRaster(
    z1,
    opacity = 0.7,                
    colorOptions = colorOptions(palette = c("cyan", "green", "yellow", "orange", "red"), 
                                domain = range(as.vector(z1[[1]]), na.rm = TRUE))  
  ) %>%
   addCircleMarkers(
    data = samples,
    radius= 1.5, 
    label = ~soc,
    color = ~paleta(soc),
    fillOpacity = 1,
    stroke = F
  ) %>%
  addLegend(
    "bottomright", 
    pal = paleta,
    values = as.vector(z1[[1]]),
    title = "IDW SOC interpolation [%]",
    opacity = 0.9
  )
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
## Warning in paleta(soc): Some values were outside the color scale and will be
## treated as NA
## Warning in paleta(soc): Some values were outside the color scale and will be
## treated as NA
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m

Este mapa muestra la interpolación de carbono orgánico del suelo (SOC) utilizando el método de distancia inversa ponderada (IDW). Los valores más altos de SOC parecen concentrarse en ciertas áreas del mapa representadas con colores cálidos (amarillo, naranja, rojo) y las áreas con SOC más bajo están representadas en tonos azules y verdes. Los puntos rojos representan los lugares donde se tomaron muestras de SOC por lo cual las zonas donde hay mayor densidad de puntos indican mayor cantidad de datos, por lo que pueden tener una interpolación más precisa. En este mapa se evidencia un problema en la escala de la leyenda puesto que muestra un rango mucho mayor al necesario y no se evidencian claramente que indican los colores en la leyenda, pero por más de que se intentó eliminar valores nulos y ajustar los rangos de la escala, no se encontró forma de corregir la leyenda sin afectar drásticamente la visualización de la interpolación.

5.3 Interpolación OK

Para el kriging ordinario (OK) se requiere un modelo de variograma para cuantificar la autocorrelación espacial de los datos. El siguiente código calcula el variograma empírico de muestras, sin covariables:

samples <- na.omit(samples)
v_emp_ok = variogram(soc ~ 1, data=samples)

Ahora lo vamos a visualizarlo gráficamente

plot(v_emp_ok)

Vamos a seguir ajustando el modelo del variograma a un variograma empírico, para esto existen varias formas, pero aquí usaremos la más sencilla: El ajuste automático mediante la función autofitVariogram del paquete automap, el cual permite evaluar varios modelos de variogramas sean esférico, exponencial, gaussiano, etc. Después de hacer la respectiva evaluación se elegirá el modelo que menor diferencia presente entre el valor empírico y el valor teórico asegurándonos que los parámetros encontrados sean los óptimos para realizar la interpolación de Kriging.

v_mod_ok = automap::autofitVariogram(soc ~ 1, as(samples, "Spatial"))

Hay que tener en cuenta que la función autofitVariogram no funciona en objetos sf y dado que sample es un raster de tipo sf, es necesario hacer la conversión, por lo que convertimos el objeto samples en un SpatialPointsDataFrame (paquete sp) Para conocer a mayor detalle los tipos de modelos que utiliza autofitVariagram, vamos a utilizar la funcion show.vgms() obteniendo la gráfica de 17 modelos que se van a comparar con los valores experimentales para identificar cuales se aproximan más.

show.vgms()

Ahora si vamos a representar gráficamente el modelo que más se aproxima al valor empírico

plot(v_mod_ok)

En el gráfico se puede ver el modelo de variograma ajustado a los datos espaciales. En este caso el modelo elegido es el Ste, que es un tipo de modelo geoestadístico. Tiene un Nugget de 2119 lo cual indica la variabilidad a escalas muy pequeñas o el error de medición, en este caso, el valor tan alto nos sugiere una mayor variabilidad. El sill de 88376, que es el umbral del variograma e indica la varianza total de los datos cuando las distancias entre puntos no son lo suficientemente grandes como para que la correlación espacial desaparezca. El rango de alcance que es 1869, es decir, que a partir de ese valor los valores dejan de estar correlacionados espacialmente y el Kappa de 0.4 que indica que la transición en la variabilidad espacial es más brusca.

la función autofitVariogram() ajusta un modelo de variograma a los datos espaciales y devuelve una lista con varios componentes. Uno de esos componentes es var_model, que contiene una tabla con los parámetros del modelo ajustado.

v_mod_ok$var_model
##   model     psill    range kappa
## 1   Nug  2118.897    0.000   0.0
## 2   Ste 86256.732 1869.145   0.4

Ahora, el modelo de variograma se puede pasar a la función gstat y podemos continuar con la interpolación de Kriging ordinario

stars.rrr_resampled <- st_downsample(stars.rrr, n = 2)
g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)
z2= predict(g2, stars.rrr_resampled)
## [using ordinary kriging]

En el código anterior se está reduciendo la resolución del objeto para agilizar el proceso, seguido a esto se define el nuevo objeto (g2) para realizar la interpolación espacial mediante Kriging Ordinario (OK) junto con la fórmula apropiada que nos permitirá indicar los valores sin covariables. En el tercer renglón aplicamos la interpolación Kriging Ordinario (g2) a un objeto ráster (stars.rrr), lo que genera un nuevo conjunto de datos interpolados (z2).

Nuevamente, vamos a crear un subconjunto del atributo de valores predichos y lo renombraremos soc, además eliminamos los valores nulos:

a <- which(is.na(z2[[1]]))
z2[[1]][a] = 0.0
names(z2) = "soc"

Ahora revisaremos el rango del valor de la predicción z2 y los datos de entrada:

range(z2$soc, na.rm = TRUE)
## [1]   98.19627 1005.15943
range(samples$soc, na.rm = TRUE)
## [1]   91.32 1193.53

Ahora vamos a visualizar los datos por medio de un histograma que nos va a permitir comparar el comportamiento de los datos y observar si tienen la misma distribución y tendencia

valores_z2 <- as.vector(z2$soc)
valores_samples <- samples$soc
df_hist <- data.frame(
  Tipo = rep(c("Z2", "Samples"), c(length(valores_z2), length(valores_samples))),
  Valores = c(valores_z2, valores_samples)
)
ggplot(df_hist, aes(x = Valores, fill = Tipo)) +
  geom_histogram(position = "identity", alpha = 0.5, bins = 15) +
  facet_wrap(~Tipo, scales = "free") +
  theme_minimal() +
  labs(title = "Distribución de Valores",
       x = "Carbono Orgánico Disponible - SOC",
       y = "Cantidad de valores")

En el histograma podemos ver como ambos objetos tienden a la misma distribución de datos Ahora vamos a seleccionar la paleta de colores con la cual graficaremos la interpolación OK

paleta1 <- colorNumeric(palette =  c("#66CD00", "#BCEE68", "#FFD700", "#FF4500", "#FF0000"), samples$soc, na.color = "transparent")

Y con la función leaflet visualizaremos el mapa interactivo de la interpolación OK

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,  # Raster sin puntos
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("#66CD00", "#BCEE68", "#FFD700", "#FF4500", "#FF0000"), 
                                  domain = range(samples$soc, na.rm = TRUE)) # Ajusta el dominio a los valores de z2$soc
    ) %>%
   addCircleMarkers(  # Aquí sí mostramos los puntos
    data = samples,
    radius = 1.5, 
    label = ~soc,
    color = ~paleta1(soc),
    fillOpacity = 1,
    stroke = FALSE
  ) %>%
    addLegend("bottomright", pal = paleta, values = samples$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
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m

Las zonas verdes del mapa indican valores más bajos de SOC, mientras que las zonas amarillas a rojas muestran concentraciones más altas de SOC, siendo lo lugares con mayor concentración la zona sur y suroeste del departamento.

6. Evaluación de resultados

6.1 Evaluación cualitativa

Vamos a seleccionar otra paleta de colores para este análisis

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

Y ahora visualizaremos el mapa interactivo

m <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(z1, colorOptions = colores, opacity = 0.8, group= "IDW")  %>%
  addGeoRaster(z2, colorOptions = colores, opacity = 0.8, group= "OK")  %>%
  # Add layers controls
  addLayersControl(
    overlayGroups = c("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
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m

6.2 Validación cruzada

Hemos generado las superficies de SOC aplicando dos enfoques distintos: IDW y Kriging ordinario. Si bien la comparación visual de los resultados es útil, es fundamental contar con un método objetivo para evaluar la precisión de la interpolación. Esto nos permite determinar con mayor certeza cuál de los métodos disponibles ofrece mejores resultados.

Para ello, empleamos la validación cruzada, una técnica estadística utilizada para medir el desempeño de un modelo predictivo. Su propósito principal es evaluar la capacidad del modelo para generalizar a datos no vistos, reduciendo el riesgo de sobreajuste (overfitting).

El proceso consiste en dividir los datos en subconjuntos (folds), entrenando el modelo con una parte y validándolo con la otra. Este procedimiento se repite varias veces con diferentes particiones, y los resultados se promedian para obtener una métrica más representativa de la precisión del modelo.

Entre las variantes más comunes de validación cruzada se encuentran: • Validación Simple (Hold-Out) • Validación Cruzada k-Fold • Validación Leave-One-Out (LOO)

En particular, en la validación Leave-One-Out: 1. Se excluye un punto de la base de datos de calibración. 2. Se realiza una predicción para ese punto. 3. El proceso se repite para todos los puntos.

Al finalizar, obtenemos una tabla con los valores observados y sus correspondientes valores predichos.

Para este caso, aplicaremos la validación cruzada Leave-One-Out mediante la función gstat.cv, que trabaja con objetos de tipo gstat.

Al ejecutar el siguiente código, se ocultarán tanto los mensajes como los resultados

## this code should be run from the console not from here
# cv2 = gstat.cv(g2)

El código anterior permite almacenar los resultados en cv2, que contiene las predicciones en los puntos de muestreo y los errores de predicción. Ahora bien, esta variable contiene las siguientes columnas o atributos: • var1.pred: Valores predichos por el modelo de Kriging (g2) en esos mismos puntos. • var1.var: varianza (solo para Kriging) • observed: Valores reales de soc en los puntos de muestreo (samples). • residual: Diferencia entre el valor observado y el predicho • zscore: Residual normalizado (se usa para detectar errores sistemáticos en la interpolación). • fold: ID de validación cruzada • geometry: Punto espacial que fue extraido para hacer la validación cruzada Ahora vamos a convertir el objeto cv2

## run from the console
# cv2 = st_as_sf(cv2)

Ahora haremos un gráfico que nos representara el error de la estimación

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

Los circulos representan la diferencia o el error que hay entre los datos reales y las estimaciones generadas en la interpolación utilizando kriging ordinario (OK)

Los círculos verdes representan los “residuales positivos”, es decir, las sobreestimaciones del modelo y los morados representan las subestimaciones del modelo o también llamadas “residuales negativos”

El tamaño de los círculos indica el tamaño del error, entre más grande el circulo señala errores más grandes de predicción y entre más pequeños indica menor error.

Teniendo en cuenta que los círculos están distribuidos sobre el área del departamento del Caquetá podemos inferir que la zona donde se genero mayor error en las estimaciones fue la zona de la cordillera y la frontera con la selva amazónica, esto se puede dar a causa de los usos agrícolas del suelo en esas zonas o las fallas o variaciones geológicas que se presentan en dichas zonas.

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

En la línea anterior usamos un código para determinar el RMSE, que es un valor que nos mide la exactitud de las predicciones realizadas por los valores. En este caso nuestro RMSE es de 64,44095.

Vamos a repetir el proceso, pero ahora con los resultados de la interpolación IDW que están en el objeto g1

## this code should be run from the console not from here
## cv1 = gstat.cv(g1)
## run from the console
# cv1 = st_as_sf(cv1)

Vamos a graficar los residuos y a determinar el valor RMSE:

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

En el gráfico podemos ver que los circulos están distribuidos por la misma zona que en el modelo anterior, sin embargo, el tamaño en este caso es mayor, lo cual se confirma con el valor RMSE que en este caso es de 73.27025

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

El modelo OK, al tener un valor más pequeño de RMSE, tiene un mejor desempeño en la predicción de los valores espaciales en comparación con el módelo IDW. Esto se debe a que un valor menor de RMSE tiene un menor error promedio y, por lo tanto, una mejor capacidad de estimar valores en ubicaciones no muestreadas.

En este orden de ideas si buscamos minimizar el error en la interpolación, el modelo OK sería la opción más adecuada, sin embargo, algo que cabe resaltar es que, el modelo IDW es más fácil y rápido de calcular, por lo que sigue siendo útil si la precisión no es la prioridad o si vamos a trabajar con grandes conjuntos de datos.

7. Guardar las salidas de interpolación

vamos a usar la función write_stars() para guardar el objeto ráster (el resultado de la interpolación) en un archivo, en este caso, un GeoTIFF (.tif)

write_stars(
  z1, dsn = "C:/Users/Karol Ruiz/Documents/GEOMATICA 2024-2/PROYECTO4/DATOS/IDW_soc_caq.tif", layer = 1
)
write_stars(
  z2, dsn ="C:/Users/Karol Ruiz/Documents/GEOMATICA 2024-2/PROYECTO4/DATOS/OK_soc_caq.tif", layer = 1
)

8. Conclusiones

La interpolación espacial es una herramienta fundamental en el campo de la geomática, ya que nos permite estimar valores de variables en ubicaciones no muestreadas a partir de datos conocidos, facilitando la lectura de fenómenos geoespaciales y permitiendo la creación de modelos predictivos que pueden ayudar en la toma de decisiones frente a la gestión territorial.

En los países en vía de desarrollo, como Colombia, el análisis de variables geomorfológicas enfrenta múltiples desafíos, principalmente debido a la limitada disponibilidad de datos espaciales de alta calidad. Esto se debe a diversos factores, entre ellos la escasez de recursos económicos y tecnológicos, la dificultad de acceso a ciertas zonas geográficas y la falta de programas sistemáticos de monitoreo ambiental.

El análisis realizado en el departamento del Caquetá demostró que la presencia de cordilleras o zonas selváticas pueden generar fluctuaciones en los valores estimados a la hora de realizar análisis con variables relacionadas a la composición de los suelos. También se expuso las zonas del departamento donde se han realizado una mayor cantidad de muestras, las cuales son al norte y nororiente del departamento. Los valores más bajos de SOC predominan en la zona centro-norte del departamento y los valores más altos se encuentran en la zona sur, esto último se puede deber a las diferentes vocaciones del suelo del departamento respeto al centro y el sur de este, la vegetación o las condiciones edafológicas.

Por último, hay que recordar la importancia de la selección del método de interpolación, pues de esta elección dependerá la precisión de las estimaciones que luego serán cruciales para el desarrollo de estrategias de manejo ambiental y planificación territorial adaptadas a las condiciones específicas de la región. En el contexto del Caquetá, se ha evidenciado que técnicas como el Kriging ofrecen estimaciones más precisas en comparación con métodos determinísticos como el IDW, debido a su capacidad para modelar la estructura de dependencia espacial de los datos.

En resumen, la correcta estimación de variables geoespaciales mediante técnicas de interpolación adecuadas es vital para comprender y gestionar eficazmente los recursos y fenómenos en regiones como el departamento del Caquetá. La elección del método de interpolación debe basarse en las características específicas de los datos y los objetivos del estudio para garantizar resultados precisos y útiles.

9. Referencias