1 Introducción

En este cuaderno se presentan dos métodos de interpolación espacial ampliamente utilizados en análisis geoespacial: Inverse Distance Weighted (IDW) y Kriging Ordinario (OK). IDW corresponde a un procedimiento determinístico que estima valores en función de la distancia a los puntos conocidos, mientras que OK es un método geoestadístico que modela la estructura espacial del fenómeno y genera predicciones con base probabilística.

Ambas técnicas se aplican para generar una superficie continua del carbono orgánico del suelo (SOC) a una profundidad entre 15 y 30 cm, a partir de muestras extraídas de la base de datos SoilGrids, la cual posee una resolución espacial de 250 metros.

2 Instalación

Primero limpiamos el puesto de trabajo

rm(list=ls())

Para manejar la información espacial en este trabajo utilizaremos diferentes paquetes especializados. El paquete terra nos permitirá cargar, procesar y manipular datos en formato raster, mientras que sf será empleado para la lectura y gestión de capas vectoriales. Los procedimientos de interpolación espacial se realizarán con ayuda de gstat y automap, dos librerías que implementan métodos clásicos como IDW y Kriging Ordinario. Finalmente, para la elaboración de mapas y representaciones gráficas recurriremos a herramientas como leaflet, tmap y ggplot2, que facilitan tanto la visualización interactiva como la generación de productos cartográficos estáticos.

Antes de comenzar con el análisis, es importante asegurarse de que todos estos paquetes estén instalados en el entorno de R. La instalación debe hacerse desde la consola, y no dentro del cuaderno RMarkdown, para evitar errores y tiempos de carga innecesarios.

# NO INSTALES PAQUETES QUE YA HAYAS INSTALADO
#paquetes = c('terra', 'sf', 'sp', 'stars', 'gstat', 'automap', 'leaflet', 'leafem')
#install.packages(paquetes)
#

Cargamos las librerias

library(terra)
## terra 1.8.70
library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library(sp)
library(stars)
## Warning: package 'stars' was built under R version 4.5.2
## Cargando paquete requerido: abind
## Warning: package 'abind' was built under R version 4.5.2
library(gstat)
## Warning: package 'gstat' was built under R version 4.5.2
library(automap)
## Warning: package 'automap' was built under R version 4.5.2
library(RColorBrewer)
## Warning: package 'RColorBrewer' was built under R version 4.5.2
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.5.2
library(leafem)
## Warning: package 'leafem' was built under R version 4.5.2

3 Lectura de los datos de entrada

En esta etapa cargaremos la información que utilizaremos para la interpolación. Como se mencionó anteriormente, los puntos muestreados corresponden a valores de carbono orgánico del suelo (SOC) a una profundidad de 15–30 cm, derivados previamente de SoilGrids 250 m utilizando un cuaderno de procesamiento independiente.

Estos datos fueron almacenados en formato GeoPackage (.gpkg) dentro de la carpeta data/, siguiendo la convención de nombres soc_{departamento}.gpkg.

Antes de proceder con el análisis, es importante asegurarnos de que el archivo existe y está accesible desde el directorio donde estamos trabajando. Para ello, verificamos su presencia mediante una simple consulta al sistema de archivos.

list.files(path="C:/Users/José Aguirre/Documents/GB2/R-studio/InformeFinal", pattern = "ValleDelCaucaDANE.gpkg")
## [1] "ValleDelCaucaDANE.gpkg"

Después de confirmar que el archivo está disponible, procedemos a cargar los puntos de muestreo. Para ello empleamos la librería sf, que permite trabajar con datos espaciales en formato vectorial de manera eficiente. Esta lectura convierte automáticamente el GeoPackage en un objeto espacial listo para ser utilizado en los procesos de interpolación.

samples <- sf::st_read("soc_valle.gpkg")
## Reading layer `soc_valle' from data source 
##   `C:\Users\José Aguirre\Documents\GB2\R-studio\InformeFinal\soc_valle.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1748 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -77.5471 ymin: 3.090434 xmax: -75.70757 ymax: 5.047424
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

Como el análisis se desarrolla dentro del departamento del Valle del Cauca, también es necesario cargar el archivo que contiene los límites municipales del área de estudio. Este GeoPackage nos permitirá visualizar la interpolación dentro del contexto administrativo correcto y realizar análisis espaciales posteriores.

munic <- sf::st_read("wgs84valle.gpkg")
## Reading layer `wgs84valle' from data source 
##   `C:\Users\José Aguirre\Documents\GB2\R-studio\InformeFinal\wgs84valle.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 42 features and 53 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -77.54977 ymin: 3.091239 xmax: -75.70724 ymax: 5.047394
## Geodetic CRS:  WGS 84

4 Exploración de los datos de entrada

Antes de realizar cualquier interpolación, es importante revisar la distribución general de los valores muestreados. Para ello, generamos un resumen estadístico de las muestras, lo cual nos permite identificar rangos, valores atípicos y posibles inconsistencias en los datos.

summary(samples)
##       soc                    geom     
##  Min.   : 14.73   POINT        :1748  
##  1st Qu.: 36.48   epsg:NA      :   0  
##  Median : 61.50   +proj=long...:   0  
##  Mean   : 92.53                       
##  3rd Qu.:161.57                       
##  Max.   :242.38
hist(samples$soc)

Ahora procedemos a redondear los valores de SOC a dos decimales, con el fin de simplificar su lectura y mejorar la presentación de los resultados.

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

Ahora realizaremos una visualización inicial de los puntos muestreados para observar cómo se distribuyen espacialmente en el territorio. Para comprender mejor el funcionamiento de leaflet y sus opciones de visualización interactiva, te recomiendo consultar la documentación disponible en su página oficial.

pal <- colorNumeric(
  c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"),
  # los colores dependen de la variable de conteo
  domain = samples$soc,
  )
# Esta es una visualización simple
# Cambia el código para adaptarlo a tus preferencias
leaflet() %>%
  addPolygons(
    data = munic,
    color = "gray",
    # establecer la opacidad del contorno
    opacity = 1,
    # establecer el grosor del trazo en píxeles
    weight = 1,
    # establecer la opacidad del relleno
    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")

5 Interpolación

5.1 Creación del objeto gstat

Para realizar cualquier tipo de interpolación espacial, el primer paso consiste en construir un objeto de clase gstat, utilizando la función del mismo nombre. Este objeto actúa como un contenedor que reúne toda la información necesaria para ejecutar el proceso de interpolación, incluyendo:

La definición del modelo que se va a emplear.

Los datos de entrenamiento, es decir, las muestras con las que se calibrará el modelo.

La función gstat() interpreta automáticamente el tipo de interpolación que queremos aplicar según los argumentos que le proporcionemos:

Si no se especifica un modelo de variograma → se aplica IDW (Inverse Distance Weighting).

Si se proporciona un variograma y no hay covariables → se ejecuta Kriging Ordinario (OK).

En este cuaderno trabajaremos principalmente con tres parámetros de gstat():

formula: indica la relación entre la variable objetivo (SOC) y posibles covariables.

data: define los datos base o muestras que alimentan el modelo.

model: corresponde al variograma (solo necesario para Kriging).

5.2 Interpolación mediante IDW

Para generar una superficie continua de SOC mediante el método IDW, basta con construir un objeto gstat especificando únicamente la fórmula de interpolación y los datos de entrada. Al no incluirse un modelo de variograma, la función asumirá automáticamente que deseamos aplicar el método de distancia inversa.

Si deseas que también prepare la versión para Kriging, la explicación del variograma o una interpretación de resultados, dímelo y te lo elaboro.

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

Ahora que ya definimos nuestro modelo de interpolación g1, podemos emplear la función predict() para generar la superficie interpolada, es decir, para estimar los valores de SOC en toda la zona de estudio.

La función predict() requiere dos elementos principales:

Un raster (o un objeto de tipo stars), que indica dónde queremos obtener las predicciones.

Un modelo de interpolación, en este caso nuestro objeto g1 creado con gstat.

El raster de entrada cumple dos funciones importantes:

Define los puntos del espacio donde se calcularán las predicciones (para IDW, Kriging u otros métodos).

Proporciona valores de covariables, aunque esto solo aplica cuando se utiliza Kriging Universal.

Para continuar, crearemos un raster base usando la librería terra, asignando a todas sus celdas un valor de 1. La extensión espacial (extent) de este raster debe abarcar completamente nuestro departamento o área de interés.

Nota: Una forma práctica de definir las dimensiones, resolución o extensión del raster es tomar como referencia cualquier otro raster que cubra el departamento (por ejemplo, un DEM o un raster descargado anteriormente).

(rrr <- terra::rast(xmin= -77.55, xmax= -75.70, ymin= 3.09, ymax= 5.05, nrows=370, ncols = 329, vals = 1, crs = "wgs84"))
## class       : SpatRaster 
## size        : 370, 329, 1  (nrow, ncol, nlyr)
## resolution  : 0.0056231, 0.005297297  (x, y)
## extent      : -77.55, -75.7, 3.09, 5.05  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : lyr.1 
## min value   :     1 
## max value   :     1

Nota: Si tu computador no cuenta con mucha memoria o capacidad de procesamiento, es recomendable trabajar con una resolución espacial más baja (por ejemplo, 0.008 o 0.01 grados). Esto reduce la cantidad de celdas del raster y hace que la interpolación sea más ligera.

Una vez creado el SpatRaster que servirá como base para realizar las predicciones, debemos convertirlo en un objeto de clase stars, ya que la función predict() del paquete gstat trabaja directamente con este tipo de estructuras.

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

La siguiente instrucción realiza la interpolación de los valores de SOC utilizando:

el modelo de interpolación almacenado en g1, y

el raster base previamente convertido a objeto stars (stars.rrr).

Es decir, esta línea genera las estimaciones continuas de carbono orgánico del suelo en cada celda del raster definido para nuestra área de estudio.

## [interpolación por distancia inversa ponderada]
## ejecutar este código toma varios minutos
z1 = predict(g1, stars.rrr)
## [inverse distance weighted interpolation]

¿Qué es z1?

El objeto z1 corresponde al resultado de la interpolación. Es decir, es un raster continuo donde cada celda contiene un valor estimado de carbono orgánico del suelo (SOC), generado a partir de los puntos de muestreo originales mediante el método IDW.

En otras palabras:

z1 ya no son los puntos individuales de SOC,

sino una superficie completa que predice el valor del suelo en todas las ubicaciones dentro del área de estudio,

usando como referencia el modelo definido en g1.

Por eso z1 es el raster interpolado que luego podemos visualizar en mapas, reclasificar o usar para análisis posteriores.

z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean  3rd Qu.     Max.   NA's
## var1.pred  15.57863 50.20229 74.95066 96.55337 145.7212 240.7343      0
## var1.var         NA       NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -77.55  0.005623 WGS 84 [x]
## y    1 370   5.05 -0.005297 WGS 84 [y]

Debemos observar que el objeto z1 contiene dos atributos internos (dos capas o columnas de información). Para trabajar únicamente con el atributo que nos interesa —la capa interpolada de SOC— es necesario extraer esa primera capa y además asignarle un nombre adecuado.

En práctica:

Subsetting nos permite seleccionar solo la primera capa del objeto z1.

Luego, renombramos esa capa para evitar confusiones y darle un nombre descriptivo que podamos usar en los mapas y análisis posteriores.

# tratando los valores NA
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  15.57863 50.20229 74.95066 96.55337 145.7212 240.7343      0
## var1.var         NA       NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -77.55  0.005623 WGS 84 [x]
## y    1 370   5.05 -0.005297 WGS 84 [y]
names(z1) = "soc_valle"

Para visualizar la superficie interpolada es necesario definir previamente una paleta de colores. Esta paleta asignará un color específico a cada rango de valores del raster, lo cual facilita interpretar la variación espacial del SOC en el departamento. Elegir una buena paleta es clave para resaltar patrones y contrastes en la superficie resultante.

# cambia este fragmento para que se ajuste a tus preferencias
paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), domain = 10:120, na.color = "transparent")

Ahora procedemos a representar gráficamente el resultado de la interpolación. Con este mapa podremos observar cómo varían espacialmente los valores estimados de SOC en todo el departamento, facilitando la interpretación visual del patrón generado por el método de interpolación aplicado.

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 10:120)
    ) %>%
   addCircleMarkers(
    data = samples,
    radius= 1.5, 
    label = ~soc,
    color = ~paleta(soc),
    fillOpacity = 1,
    stroke = F
  ) %>%
    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
## 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  # Imprimir el mapa

5.3 Interpolación mediante Kriging Ordinario (OK)

A diferencia del método IDW, los procedimientos de Kriging requieren la construcción previa de un modelo de variograma. Este modelo describe cómo cambia la similitud entre puntos a medida que aumenta la distancia entre ellos, permitiendo asignar pesos óptimos cuando se realizan predicciones espaciales.

El variograma es, por tanto, una herramienta fundamental para cuantificar de forma objetiva el patrón de autocorrelación espacial presente en los datos.

Como primer paso, es necesario calcular y visualizar el variograma empírico, el cual muestra la relación entre la distancia y la disimilitud promedio entre los valores muestreados.

Para ello utilizamos la función variogram(), que requiere dos elementos:

formula: indica la variable dependiente que queremos interpolar, así como los posibles predictores (si existieran).

data: la capa de puntos que contiene los valores observados y sus atributos.

La siguiente instrucción calcula el variograma empírico de las muestras, sin incluir covariables:

v_emp_ok = variogram(soc ~ 1, data=samples)

Ahora vamos a graficar el variograma empírico para visualizar cómo varía la dispersión de los valores de SOC en función de la distancia entre los puntos muestreados. Esta gráfica nos permitirá identificar tendencias de autocorrelación espacial que luego utilizaremos para ajustar un modelo teórico de variograma.

plot(v_emp_ok)

Existen distintos métodos para ajustar un modelo teórico al variograma empírico. En este cuaderno utilizaremos la alternativa más sencilla: el ajuste automático mediante la función autofitVariogram del paquete automap. Esta herramienta selecciona de forma eficiente los parámetros del modelo (nugget, sill y rango), permitiendo obtener un variograma teórico adecuado sin necesidad de hacerlo manualmente.

# asegúrate de entender los parámetros necesarios para ejecutar esta función de ajuste
v_mod_ok = automap::autofitVariogram(soc ~ 1, as(samples, "Spatial"))

La función selecciona automáticamente el tipo de modelo que mejor se ajusta al variograma empírico y ajusta sus parámetros (como nugget, sill y rango). Si quieres explorar los distintos modelos de variograma disponibles, puedes usar show.vgms(), que muestra las formas comunes utilizadas en geoestadística.

Es importante tener en cuenta que autofitVariogram no opera directamente sobre objetos de tipo sf. Por esta razón, previamente convertimos el objeto samples a un SpatialPointsDataFrame del paquete sp.

Una vez ajustado el modelo, es posible visualizarlo utilizando la función plot().

plot(v_mod_ok)

El objeto generado por la función no es un único modelo, sino una lista que agrupa varios elementos, entre ellos el variograma empírico y el variograma ajustado. El componente llamado $var_model dentro de esa lista es el que contiene el modelo de variograma final, es decir, el que será utilizado posteriormente para el proceso de kriging.

v_mod_ok$var_model
##   model        psill    range kappa
## 1   Nug 2.022884e+02      0.0   0.0
## 2   Ste 3.045839e+07 248639.4   0.6
## [usando kriging ordinario]
## Esto toma varios minutos 
## (o incluso una eternidad dependiendo del tamaño de celda del ráster y de los recursos de tu computador)
g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)
z2= predict(g2, stars.rrr)
## [using ordinary kriging]

Nuevamente, vamos a extraer únicamente la capa que contiene los valores predichos y asignarle un nombre más apropiado para su interpretación.

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

Las predicciones generadas mediante Kriging Ordinario (OK) se visualizan en el siguiente mapa.

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 10:120)
    ) %>%
   addCircleMarkers(
    data = samples,
    radius= 1.5, 
    label = ~soc,
    color = ~paleta(soc),
    fillOpacity = 1,
    stroke = F
  ) %>%
    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
## 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  # Imprimir el mapa

6 Evaluación de los resultados

6.1 Evaluación cualitativa

colores <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), domain = 10:120, na.color = "transparent")
m <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(z1, colorOptions = colores, opacity = 0.8, group= "IDW")  %>%
  addGeoRaster(z2, colorOptions = colores, opacity = 0.8, group= "OK")  %>%
  # Añadir controles de capas
  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  # Imprime el mapa

7 Validación cruzada

Hasta este punto hemos generado superficies continuas de SOC empleando dos enfoques distintos: IDW y Kriging Ordinario. Aunque la inspección visual permite tener una idea general del comportamiento espacial de cada método, esta revisión no es suficiente para determinar cuál de ellos ofrece mejores estimaciones. Para lograr quietud científica, es necesario aplicar una evaluación cuantitativa que nos permita medir de manera objetiva la precisión de los resultados.

La técnica más utilizada para este propósito es la validación cruzada Leave-One-Out (LOOCV). Este procedimiento consiste en:

Excluir temporalmente un punto del conjunto de datos.

Predecir el valor de ese punto usando únicamente los demás puntos disponibles.

Comparar la predicción con el valor real.

Repetir el proceso para cada punto de la muestra.

Al finalizar, obtenemos una tabla donde para cada ubicación existen dos valores:

el SOC observado,

y el SOC estimado por el modelo.

Esta información permite calcular métricas de error y comparar rigurosamente los dos métodos de interpolación.

En R, este procedimiento puede ejecutarse fácilmente mediante la función gstat.cv(), que recibe como argumento un objeto creado con gstat.

Cuando escribas el bloque de código correspondiente en tu notebook, recuerda ocultar los mensajes y los resultados para mantener el documento ordenado.

### 
## Este codigo debe ser ejecutado desde la consola, no desde aqui
#cv2 = gstat.cv(g2)

Los resultados generados por gstat.cv contienen varias columnas que permiten evaluar en detalle el desempeño del método de interpolación. Entre las más importantes se encuentran:

var1.pred: el valor estimado por el modelo para cada punto.

var1.var: la varianza asociada a la predicción (solo disponible cuando se usa Kriging).

observed: el valor real medido en campo o tomado de la fuente original.

residual: la diferencia entre lo observado y lo predicho (observado – predicho).

zscore: puntaje estandarizado que evalúa si la predicción está dentro de lo esperado según la varianza (solo en Kriging).

fold: identificador del proceso de validación cruzada para cada punto.

Para poder trabajar con estos resultados de manera más cómoda, es necesario transformar el objeto cv2 en un formato tabular o espacial según nuestras necesidades. A continuación se muestra cómo realizar dicha conversión:

# Ejecuta desde la consola
#cv2 = st_as_sf(cv2)

Ahora procedamos a visualizar los residuos obtenidos en la validación cruzada. Esta gráfica permite identificar si el modelo tiende a sobreestimar o subestimar los valores de SOC en distintos puntos del territorio, lo cual es clave para evaluar su desempeño.

# Ejecuta desde la consola
#sp::bubble(as(cv2[, "OK residuals"], "Spatial"))
# Este es el valor de RMSE para la interpolación OK
# ejecútalo desde la consola
#sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))

Ahora registra el valor del RMSE obtenido para los resultados de la interpolación mediante Kriging Ordinario (OK).

Luego, realiza el mismo procedimiento para el método IDW y anota su RMSE correspondiente.

Una vez tengas ambos valores, compáralos y analiza en qué se diferencian.

El método de interpolación más confiable es aquel que presenta el menor valor de RMSE (Root Mean Square Error) durante la validación cruzada. El RMSE mide cuánto se desvían, en promedio, las predicciones del modelo respecto a los valores observados. Por lo tanto, cuanto más bajo sea el RMSE, mayor precisión y mejor desempeño tendrá el método.

Si el RMSE obtenido con Kriging Ordinario (OK) es menor que el obtenido con IDW, entonces OK es el método más confiable, ya que logra capturar mejor la variabilidad espacial mediante un modelo de variograma. En cambio, si el RMSE de IDW es menor, significa que la estructura espacial del fenómeno no está bien representada por el variograma y que un enfoque determinístico funciona mejor en este caso.

8 Guardar los resultados de la interpolación

write_stars(
  z1, dsn = "soc_valle.tif", layer = 1
)
write_stars(
  z2, dsn = "soc_valle.tif", layer = 1
)

9 Conclusiones

El desarrollo de este cuaderno permite comprender de manera práctica cómo funcionan los procesos de interpolación espacial y por qué son tan importantes en geomática. A lo largo del ejercicio aprendí a preparar, explorar y visualizar datos espaciales en R, así como a aplicar dos métodos de interpolación distintos: IDW y Kriging Ordinario. Descubrí que ambos generan superficies continuas a partir de puntos muestreados, pero utilizan enfoques conceptualmente diferentes: mientras IDW se basa únicamente en la distancia, el Kriging incorpora la estructura espacial de los datos mediante el variograma. Además, el uso de la validación cruzada y del RMSE me enseñó la importancia de evaluar los modelos cuantitativamente y no depender solo de la apariencia visual de los mapas. En general, este cuaderno reforzó mi capacidad para trabajar con datos espaciales en R, interpretar resultados y entender las ventajas y limitaciones de cada técnica, lo que me ayuda a tomar decisiones informadas al analizar información geográfica.

10 Referencia

Lizarazo, I. (2025). Interpolación espacial. Disponible en: https://rpubs.com/ials2un/Spatial_Interpolation

sessionInfo()
## R version 4.5.1 (2025-06-13 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## 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] leafem_0.2.5       leaflet_2.2.3      RColorBrewer_1.1-3 automap_1.1-20    
##  [5] gstat_2.1-4        stars_0.6-8        abind_1.4-8        sp_2.2-0          
##  [9] sf_1.0-21          terra_1.8-70      
## 
## loaded via a namespace (and not attached):
##  [1] generics_0.1.4          sass_0.4.10             class_7.3-23           
##  [4] KernSmooth_2.23-26      lattice_0.22-7          digest_0.6.37          
##  [7] magrittr_2.0.4          evaluate_1.0.5          grid_4.5.1             
## [10] fastmap_1.2.0           plyr_1.8.9              jsonlite_2.0.0         
## [13] e1071_1.7-16            reshape_0.8.10          DBI_1.2.3              
## [16] crosstalk_1.2.2         scales_1.4.0            codetools_0.2-20       
## [19] jquerylib_0.1.4         cli_3.6.5               rlang_1.1.6            
## [22] units_1.0-0             intervals_0.15.5        base64enc_0.1-3        
## [25] cachem_1.1.0            yaml_2.3.10             FNN_1.1.4.1            
## [28] raster_3.6-32           tools_4.5.1             parallel_4.5.1         
## [31] dplyr_1.1.4             ggplot2_4.0.0           spacetime_1.3-3        
## [34] 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        
## [40] classInt_0.4-11         leaflet.providers_2.0.0 htmlwidgets_1.6.4      
## [43] pkgconfig_2.0.3         pillar_1.11.1           bslib_0.9.0            
## [46] gtable_0.3.6            glue_1.8.0              Rcpp_1.1.0             
## [49] tidyselect_1.2.1        tibble_3.3.0            xfun_0.53              
## [52] rstudioapi_0.17.1       knitr_1.50              farver_2.1.2           
## [55] htmltools_0.5.8.1       rmarkdown_2.30          xts_0.14.1             
## [58] compiler_4.5.1          S7_0.2.0