Interpolacion del Cauca

1. Introducción:

Este cuaderno de interpolacion espacial nos va aservir a ver las funciones de distintas librerias como lo son: library(sp),library(terra),library(sf),library(stars),library(gstat),library(automap), library(leaflet),library(leafem),library(ggplot2),library(dplyr). Todo esto nos va a ayudar a crear una vision de lo que es una interpolacion espacial del departamento del Cauca. Para lo sigueinte podemos decir que en este cuaderno vamos a ver lo que es la interpolación lo cual es una técnica fundamental en geomática que permite estimar valores desconocidos de una variable geoespacial,como la altitud, la temperatura o la humedad del suelo, en ubicaciones donde no se han realizado mediciones directas, a partir de un conjunto de datos conocidos. Por lo anteiror vamos a ver lo que es el metdod OK (ordinary Kriging) y el metodo IDW.

2. Cofiguración

En esta parte vamos a borrar todos los datos guadados, para que nos quede libre la memoria que vamos a utilizar con los nuevos valores y nuevas variables.

rm(list =ls())

Como lo mensiocnamos anteriormente toca primero utilizar las siguientes librerias.library(sp),library(terra),library(sf),library(stars),library(gstat),library(automap),library(leaflet),library(leafem),library(ggplot2), library(curl),library(dplyr), para luego desarrollar el codigo.

library(sp)
library(terra)

## terra 1.8.29

library(sf)

## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

library(stars)

## Loading required package: abind

library(gstat)
library(automap)
library(leaflet)
library(leafem)
library(ggplot2)
library(dplyr)

## 
## Attaching package: '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)

## Warning: package 'curl' was built under R version 4.3.3

## Using libcurl 8.3.0 with Schannel

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

3. Datos de entrada.

En esta parte lo que vamos a realizar en el ingreso de datos Raster al sistema para luego trabajar con ellos, estos son los datos de carbono organico de suelos (soc).

archivo <- ("soc_igh_15_30.tif")
(soc <- rast(archivo))

## class       : SpatRaster 
## dimensions  : 1056, 969, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8677750, -8435500, 106750, 370750  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine 
## source      : soc_igh_15_30.tif 
## name        : soc_igh_15_30

En este paso vamos a convertir las unidades de SOC dividiendolas sobre 10.

soc.perc <-  soc/10

En esta parte vamos a proyectar a datos geograficos con la funcion “project”, los valores dados se dan en grados.

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

## class       : SpatRaster 
## dimensions  : 1064, 993, 1  (nrow, ncol, nlyr)
## resolution  : 0.002229862, 0.002229862  (x, y)
## extent      : -77.95047, -75.73621, 0.9579311, 3.330504  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : soc_igh_15_30 
## min value   :      14.63106 
## max value   :     269.90146

Aqui lo que hacemos es convertir un objeto de spat raster a stars, que es otra forma de utilizar datos en R.

stars.soc = st_as_stars(geog.soc)

Aqui lo que vamos a hacer es un ploteo utilizando leaflet y leafem. El mapa muestra los valores de carbono organico del suelo en la region del departamento del Cauca

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      stars.soc,
      opacity = 0.8,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 14:270)
    ) 
#
m  # Print the map

4. Toma de muestras

Aqui se toma una muestra aleatoria de 500 puntos del raster de carbono organico del suelo y luego de estos se los convierte a sf para un analiis espacial.

set.seed(123456)

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

##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 500, 1  (geometries, attributes)
##  extent      : -77.94935, -75.74179, 0.9679655, 3.324929  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : soc_igh_15_30
##  type        :         <num>
##  values      :         69.41
##                        51.87
##                           NA

(muestras <- sf::st_as_sf(samples))

Aqui eliminamos los valores que no tienen valores que son NA, osea los valores que no tienen información para los cuales son NA.

nmuestras <- na.omit(muestras)

Aqui vamos obteniendo las coordenadas de longitud y latitud de cada punto. Y se extrae los valores de carbono organico del suelo.

longit <- st_coordinates(muestras)[,1]
latit <- st_coordinates(muestras)[,2]
soc <- muestras$soc_igh_15_30

id <- seq(1,500,1) 

Aqui creamos un data frame, con las siguintes columnas el identificador (id), la longitud, latitud y el valor de carbono organico del suelo.

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

Elimina culquier fila del dataframe que contenga NA en algunas de sus columnas.

sitios <- na.omit(sitios)

Ver las primeras filas

head(sitios)

Se crea un mapa intrativo que muestra la distribución espacial del SOC como los puntos de meustreo de cada lugar del departamento el cual nos da un valor por cada municipio que se representa con un circulo de color naranja.

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      stars.soc,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("black", "yellow", "cyan", "orange"), 
                                  domain = 14:270)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions())
m  # Print the map

5.Interpolación IDW

En este tipo de interpolación no necesita de un semi variograma.

En g1 se esta creando un objeto (gstat),el soc_igh_15_30 es una variable este esta en raster, nmuestras es el conjunto de puntos muestreados de SOC.

g1 = gstat(formula = soc_igh_15_30 ~ 1, data = nmuestras)

Posteriormente, se define el raster de plantilla para poner los datos. Aqui el 4 nos quiere decir que este coge por cada 4 cuadros que da es una celda que hace, esto hace que redusca la resolución del raster original.

rrr = aggregate(geog.soc, 4)

Se imprime la variable previamente diseñada para ver y comprobar que todos los valores sean correctos para el departamento.

rrr

## class       : SpatRaster 
## dimensions  : 266, 249, 1  (nrow, ncol, nlyr)
## resolution  : 0.008919447, 0.008919447  (x, y)
## extent      : -77.95047, -75.72952, 0.9579311, 3.330504  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : soc_igh_15_30 
## min value   :      15.96767 
## max value   :     258.72571

El sigueinte paso es asignar un valor constante a las celdas este es de 1, lo que va a repercutir en la grilla.

values(rrr) <-1

Aqui renombramos la capa.

names(rrr) <- "valor"

Verificación del raster creado, los valores son 1.

rrr

## class       : SpatRaster 
## dimensions  : 266, 249, 1  (nrow, ncol, nlyr)
## resolution  : 0.008919447, 0.008919447  (x, y)
## extent      : -77.95047, -75.72952, 0.9579311, 3.330504  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : valor 
## min value   :     1 
## max value   :     1

Convertimos el raster rrr a un objeto stars con la funcion de st_as_stars.

stars.rrr = st_as_stars(rrr)

Aqui interpolamos que esta definida en g1 sobre stars.rrr, todos los objeto se guardan en z1 que es un objeto stars.

z1 = predict(g1, stars.rrr)

## [inverse distance weighted interpolation]

Vemos que en z1 en var1.pred tiene los valores interpolados de SOC y en var1.var la vrianza de la predicción.

z1

## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median    Mean  3rd Qu.     Max.  NA's
## var1.pred  18.10718 49.76386 57.72762 67.8056 73.04993 198.2528     0
## var1.var         NA       NA       NA     NaN       NA       NA 66234
## dimension(s):
##   from  to offset     delta                       refsys x/y
## x    1 249 -77.95  0.008919 +proj=longlat +datum=WGS8... [x]
## y    1 266  3.331 -0.008919 +proj=longlat +datum=WGS8... [y]

Aqui podemos ver la primera capa del raster interpolado.

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

Imprimimos la variable crada z1

z1

## stars object with 2 dimensions and 1 attribute
## attribute(s):
##          Min.  1st Qu.   Median    Mean  3rd Qu.     Max.
## soc  18.10718 49.76386 57.72762 67.8056 73.04993 198.2528
## dimension(s):
##   from  to offset     delta                       refsys x/y
## x    1 249 -77.95  0.008919 +proj=longlat +datum=WGS8... [x]
## y    1 266  3.331 -0.008919 +proj=longlat +datum=WGS8... [y]

5.1 Creación de paleta de Color

Aqui creamos una paleta de colores para representar los valores de SOC interpolados.Las celdas NA se presentan transparentes.

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

Creación del mapa. En este podemos ver la superficie interpolada de SOC estimada en IDW,podemos ver los puntos reales muestreados en el mapa, una leyenda para una interpretación de los colores. Los colores se componen de naranja, amarillo, cyan y verde.

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 18:200)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal=paleta, values= z1$soc,
    title = "IDW SOC interpolation [%]"
    )
m  # Print the map

6 OK INTERPOLATION (ordinary krying)

v_emp_ok = variogram(soc_igh_15_30 ~ 1, data=nmuestras)

Creamos un ploteo de la variable v_emp_ok anteriormente creada haciendo variograma de las muestras, para mirar la tendencia de estos datos que nos dan.

plot(v_emp_ok)

En esta sección del análisis se realiza el cálculo y ajuste automático de un variograma experimental, el cual es una herramienta fundamental en geoestadística para analizar la estructura espacial de una variable. En este caso, la variable estudiada es el carbono orgánico del suelo en el horizonte de 15 a 30 cm de profundidad, representada por soc_igh_15_30. El código utiliza la función autofitVariogram() del paquete automap, la cual permite generar un variograma experimental a partir de los datos de muestreo y, de forma automatizada, ajustar un modelo teórico (como esférico, exponencial o gaussiano) que mejor represente la relación espacial entre los puntos. El objeto nmuestras, que contiene los datos georreferenciados, es transformado a la clase “Spatial” para asegurar compatibilidad con la función. La fórmula soc_igh_15_30 ~ 1 indica que se ajusta un modelo simple, sin variables explicativas adicionales.

v_mod_ok = autofitVariogram(soc_igh_15_30 ~ 1, as(nmuestras, "Spatial"))

Una vez ejecutada la función, el resultado (v_mod_ok) contiene tanto el variograma experimental como el modelo ajustado. Este resultado se visualiza con plot(v_mod_ok), lo que permite evaluar gráficamente la calidad del ajuste. El modelo resultante será clave para aplicar técnicas de interpolación espacial como el kriging, permitiendo estimar valores de carbono del suelo en ubicaciones no muestreadas, basándose en la estructura espacial observada en los datos originales.

plot(v_mod_ok)

El modelo de variograma muestra cómo varía el carbono orgánico del suelo (SOC) según la distancia entre muestras. Tiene un rango de unos 101 metros, lo que indica que los datos están correlacionados espacialmente dentro de esa distancia. El efecto nugget es bajo, lo que sugiere buena calidad en los datos. Este modelo sirve para aplicar kriging y estimar valores de SOC en áreas donde no se tomaron muestras.

v_mod_ok$var_model

Se aplica kriging ordinario para estimar valores de carbono orgánico del suelo (SOC) en zonas no muestreadas, usando la estructura espacial capturada por el variograma.

g2 = gstat(formula = soc_igh_15_30 ~ 1, model = v_mod_ok$var_model, data = nmuestras)
z2= predict(g2, stars.rrr)

## [using ordinary kriging]

Se extrae la capa de var1.pred y se le cambia el nombre por soc.

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

Se genera un mapa interactivo utilizando leaflet, integrando los valores interpolados de carbono orgánico del suelo (SOC) y los puntos originales de muestreo. La capa raster muestra la distribución espacial del SOC mediante una paleta de colores personalizada, donde los tonos van desde naranja (valores bajos) hasta verde (valores altos). Los puntos de muestreo se presentan como marcadores agrupados, cada uno con información detallada en popups. Esta visualización facilita la exploración espacial del SOC y aporta una herramienta efectiva para comunicar los resultados del análisis geoestadístico.

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 18:200)
    ) %>%
  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

m  # Print the map

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

El mapa interactivo generado permite contrastar visualmente tres representaciones del carbono orgánico del suelo (SOC): los valores reales obtenidos en campo, los estimados mediante interpolación IDW y los obtenidos por kriging ordinario. Esta herramienta facilita el análisis espacial, revelando diferencias en cómo cada método distribuye los valores estimados a lo largo del territorio.Al alternar entre capas, se observa cómo IDW suaviza los valores basándose solo en la distancia entre puntos, mientras que el krigeado aprovecha la estructura espacial identificada en el variograma, generando una estimación más ajustada a la variabilidad real del suelo. La visualización permite identificar zonas con alta concentración de SOC, transiciones entre tipos de suelo y posibles áreas críticas para intervención agrícola o conservación. Este enfoque comparativo potencia la toma de decisiones, al mostrar claramente qué método refleja mejor la distribución espacial del SOC en el área de estudio.

m <- 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 [%]"
)
m  # Print the map

7. Validacion cruzada para sacar los errores.

Se aplica una validación cruzada con gstat.cv() para evaluar el desempeño del modelo de kriging. Este procedimiento predice los valores de carbono orgánico del suelo (SOC) dejando fuera cada muestra, uno a uno.

### Este código se corre desde la consola, linea por línea.
#cv1 = gstat.cv(g1)
#cv2 = gstat.cv(g2)

#cv1 = na.omit(cv1)

#cv1 = gstat.cv(g1)

Aqui se omiten los valores que o tienen información en cv1.

#cv1 = na.omit(cv1)

#cv1

cv1(cross validation)

Para evaluar la precisión del modelo de interpolación por kriging, se transforman los resultados de la validación cruzada (cv1) en un objeto espacial sf, lo que permite una visualización geográfica directa de los residuos (errores de predicción). Luego, se genera un gráfico de burbujas con bubble(), donde cada punto representa un sitio de muestreo y el tamaño de la burbuja indica la magnitud del error. Este gráfico facilita la identificación de áreas donde el modelo tuvo menor desempeño, permitiendo ajustar futuras estrategias de interpolación o reconocer zonas con mayor variabilidad del suelo. Es una herramienta clave para validar y comunicar la calidad del análisis espacial realizado.

#cv1 = st_as_sf(cv1)

Lo que podemos interpretar de la imagen es que tenemos unos datos mas impresisos a medida que nos alejsmos del centro del departamento, son datos mas impresisos que los datos que estan en el centro. Y aun mas se encuantran datos mas impresisos en los lugares como los de la costa pacifica. Esto se dice que se espera que los datos no sean muy grandes ni muy pequeños, osea que no se alejen del 10% de sobreestimación.

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

Este valor indica que, en promedio, las estimaciones del modelo se desvían ±21 unidades del valor real de carbono orgánico del suelo (SOC). Un RMSE bajo sugiere que el modelo tiene buena precisión; en este caso, dependerá del rango total de valores de SOC para juzgar si ese error es aceptable.

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

#cv2 = gstat.cv(g2)

Analisis de burbujas de OK

En este analizamos con cun diagrama de burbujas la presicion de los datos que se dieron en el mapa y se mira que tan acertados son y por lo cual que tan confiable es hacerlo.

#cv2 = na.omit(cv2)

#cv2

#cv2 = st_as_sf(cv2)

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

El gráfico muestra la distribución espacial de los residuos obtenidos tras aplicar un modelo de interpolación en el departamento del Cauca. Los residuos positivos, representados en color verde, indican que el modelo subestimó los valores reales, mientras que los residuos negativos (en rosa) muestran sobreestimación. El tamaño de los puntos refleja la magnitud del error en cada ubicación: los puntos más grandes señalan áreas donde la discrepancia entre el valor observado y el estimado es mayor. Este tipo de visualización es útil para identificar zonas donde el modelo presenta menor precisión, permitiendo ajustes como la incorporación de más datos de entrada, la modificación de parámetros o el uso de métodos alternativos de interpolación. Evaluar los residuos es esencial para mejorar la calidad del análisis espacial y asegurar que las decisiones basadas en estos datos sean confiables.

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

El modelo de interpolación aplicado en el departamento del Cauca presenta una desviación promedio de 17.54 unidades respecto a los valores observados, lo cual indica una diferencia moderada entre los datos reales y los estimados. Si bien este valor no es excesivamente alto, los residuos espaciales muestran que hay zonas específicas donde el modelo subestima o sobreestima con mayor intensidad, lo que sugiere que el ajuste no es uniforme en todo el territorio.

sessionInfo()

## R version 4.3.2 (2023-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## 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_5.2.1     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-20     
## [11] terra_1.8-29   sp_2.2-0      
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.4       xfun_0.42          bslib_0.6.1        raster_3.6-32     
##  [5] htmlwidgets_1.6.4  lattice_0.21-9     vctrs_0.6.5        tools_4.3.2       
##  [9] crosstalk_1.2.1    generics_0.1.3     parallel_4.3.2     tibble_3.2.1      
## [13] proxy_0.4-27       spacetime_1.3-3    fansi_1.0.6        highr_0.10        
## [17] xts_0.14.1         pkgconfig_2.0.3    KernSmooth_2.23-22 lifecycle_1.0.4   
## [21] farver_2.1.1       compiler_4.3.2     FNN_1.1.4.1        munsell_0.5.1     
## [25] codetools_0.2-19   htmltools_0.5.7    class_7.3-22       sass_0.4.8        
## [29] yaml_2.3.8         pillar_1.9.0       jquerylib_0.1.4    ellipsis_0.3.2    
## [33] classInt_0.4-11    cachem_1.1.0       tidyselect_1.2.1   digest_0.6.34     
## [37] fastmap_1.2.0      grid_4.3.2         colorspace_2.1-0   cli_3.6.2         
## [41] magrittr_2.0.3     base64enc_0.1-3    utf8_1.2.4         e1071_1.7-16      
## [45] withr_3.0.0        scales_1.3.0       rmarkdown_2.25     zoo_1.8-12        
## [49] png_0.1-8          evaluate_0.23      knitr_1.45         rlang_1.1.3       
## [53] Rcpp_1.0.12        glue_1.7.0         DBI_1.2.3          rstudioapi_0.17.1 
## [57] reshape_0.8.10     jsonlite_1.8.8     R6_2.5.1           plyr_1.8.9        
## [61] intervals_0.15.5   units_0.8-7