1. INTRODUCCIÓN

Este cuaderno ilustra dos técnicas de interpolación espacial: Inverse Distance Weighted (IDW) y Ordinary Kriging (OK). El método de IDW asigna pesos a los puntos de datos conocidos en función de su distancia al punto de interpolación. Los puntos más cercanos tienen mayor influencia en el valor interpolado. Si los datos están bien distribuidos y la variabilidad espacial es relativamente baja, el método de IDW puede proporcionar resultados precisos. Por otro lado, el método de Kriging Ordinario (OK) utiliza un modelo estadístico basado en la variograma para estimar los valores desconocidos. El OK tiene en cuenta tanto la distancia como la estructura espacial de los datos. Esto significa que puede proporcionar resultados más precisos cuando hay una variabilidad espacial alta o cuando se presentan patrones espaciales complejos. En este caso, las dos tecnicas serán ilustradas con los datos correspondientes a el carbono orgánico en el suelo del departamento del Meta.

2. CONGIFURACIÓN

Limpiar la memoria de R.

rm(list=ls())

Luego de limpiar la memoria, se cargan las siguientes librerias.

library(sp)
## Warning: package 'sp' was built under R version 4.2.3
library(terra)
## Warning: package 'terra' was built under R version 4.2.3
library(sf)
## Warning: package 'sf' was built under R version 4.2.3
library(stars)
## Warning: package 'stars' was built under R version 4.2.3
library(gstat)
## Warning: package 'gstat' was built under R version 4.2.3
library(automap)
## Warning: package 'automap' was built under R version 4.2.3
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.2.3
library(leafem)
## Warning: package 'leafem' was built under R version 4.2.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.2.3

Instalar la libreria “curl”

library(curl)
## Warning: package 'curl' was built under R version 4.2.3
h <- new_handle()
handle_setopt(h, http_version = 2)

3. LEER LOS DATOS DE ENTRADA

Se debe crear un conjunto de datos con el fin de simular los datos reales del mundo. En base a lo anterior, se lee la capa de carbono orgánico en el suelo (SOC), la cual se obtuvo de ISRIC.

archivo <- ("C:/GB2/Ejercicio1/ejercicio2/Proyecto5/CARBONO_ORGANICO/SOC_META.tif")
(soc <- rast(archivo))
## class       : SpatRaster 
## dimensions  : 1414, 1654, 1  (nrow, ncol, nlyr)
## resolution  : 0.002324104, 0.002324354  (x, y)
## extent      : -74.9196, -71.07554, 1.571639, 4.858276  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : SOC_META.tif 
## name        : SOC_META

Ahora se conviernten los datos de la capa SOC a porcentajes.

soc.perc <-  soc/10

¿Qué es el CRS en los datos reales del mundo?

En esta caso, parece que se necesita una transformación de CRS a el WGS84 CRS:

geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(soc.perc, geog))
## class       : SpatRaster 
## dimensions  : 1414, 1654, 1  (nrow, ncol, nlyr)
## resolution  : 0.00232421, 0.00232421  (x, y)
## extent      : -74.9196, -71.07536, 1.571843, 4.858276  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : SOC_META 
## min value   :   0.0000 
## max value   : 172.6931

El siguiente paso es covertir la capa SpatRaster en un objeto stars:

stars.soc = st_as_stars(geog.soc)

El resultado del mapa es el siguiente:

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

4. MUESTREO DEL MUNDO

En este caso se usará una muestra de aproximadamente 500 lugares de los datos del “mundo real”, utilizando una muestra ubicada aleatoriamente.

set.seed(123456)

MUesta de los 500 sitios aleatorios.

(samples <- spatSample(geog.soc, 500, "random", as.points=TRUE))
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 500, 1  (geometries, attributes)
##  extent      : -74.91844, -71.07885, 1.57533, 4.845493  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : SOC_META
##  type        :    <num>
##  values      :       NA
##                   34.64
##                   42.66

Las coordenadas del departamento del Meta, corresponden a: (x min) = -74.91844 (x max) = -71.07885, (y min) = 1.57533, (y mas) = 4.845493

Ahora, se necesita convertir el spatvector object en un objeto simple característico.

(muestras <- sf::st_as_sf(samples))
nmuestras <- na.omit(muestras)

A continuación, se visualizan las muestras:

longit <- st_coordinates(muestras)[,1]
latit <- st_coordinates(muestras)[,2]
soc <- muestras$SOC_META
id <- seq(1,500,1) 
(sitios <- data.frame(id, longit, latit, soc))

Se remueven los valores de NA:

sitios <- na.omit(sitios)
head(sitios)

Se visualizan las muestras en un mapa intectativo:

m <- 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$soc, clusterOptions = markerClusterOptions())
m  # Print the map

Ahora, se puede continuar con los paso de interpolación.

5. INTERPOLACIÓN

5.1 CREACION DE UN OBJETO GSTAT

Para interpolar, primero se debe crear un objeto “gstat”, usando la funcion “gstat”. Un objeto “gstat” contiene toda la informacion necesaria para poder realizar la interpolacion espacial. - La definición del modelo. - La calibración de los datos.

5.2 INTERPOLACÓN CON EL METODO DE IDW

Para interpolar la propiedad de Carbono organíco del suelo (soc) del Meta usando el metodo de IDW, se crea el siguiente objeto gstat.

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

EL siguiente paso es crear un objeto raster con valores de celda iguales a 1:

# a simple copy
rrr = aggregate(geog.soc, 4)

¿Qué es rrr?

rrr
## class       : SpatRaster 
## dimensions  : 354, 414, 1  (nrow, ncol, nlyr)
## resolution  : 0.009296839, 0.009296839  (x, y)
## extent      : -74.9196, -71.07071, 1.567195, 4.858276  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : SOC_META 
## min value   :   0.0000 
## max value   : 154.5373

Definir nuevos valores:

values(rrr) <-1

Definir nuevos nombres:

names(rrr) <- "valor"

Ahora, ¿qué es rrr?

rrr
## class       : SpatRaster 
## dimensions  : 354, 414, 1  (nrow, ncol, nlyr)
## resolution  : 0.009296839, 0.009296839  (x, y)
## extent      : -74.9196, -71.07071, 1.567195, 4.858276  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : valor 
## min value   :     1 
## max value   :     1
stars.rrr = st_as_stars(rrr)

Por ejemplo, la siguiente expresión interpola los valores “SOC” correpondientes al departamento del Meta según el modelo definido en g1 y la plantilla ráster definida en stars.rrr:

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

¿Qué es z1?

z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean  3rd Qu.     Max.   NA's
## var1.pred  4.039506 34.06714 36.17594 41.23372 42.01884 143.6486      0
## var1.var         NA       NA       NA      NaN       NA       NA 146556
## dimension(s):
##   from  to   offset       delta                       refsys x/y
## x    1 414 -74.9196  0.00929684 +proj=longlat +datum=WGS8... [x]
## y    1 354  4.85828 -0.00929684 +proj=longlat +datum=WGS8... [y]

Se puede crear un subconjunto solo del primer atributo y cambiarle el nombre a “soc”:

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

Se necesita una paleta de colores para el mapa.

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

El ráster SOC interpolado correspondiente al departamento del meta, utilizando IDW, se muestra en la siguiente figura:

m <- 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 = "Interpolación del Meta por IDW [%]"
    )
## 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

5.3 INTERPOLACION OK

Los métodos kriging requieren un modelo de variograma. El modelo de variograma es una forma objetiva de cuantificar el patrón de autocorrelación en los datos y asignar pesos en consecuencia al hacer predicciones.

Como primer paso, se puede calcular y examinar el variograma empírico usando la funcion de variograma.

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

Se traza el variogramaLet’s plot the variogram:

plot(v_emp_ok)

Hay varias formas de ajustar un modelo de variograma a un variograma empírico. En este caso se usará el más simple: ajuste automático usando la función autofitVariogram del paquete automap:

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

El modelo ajustado se puede trazar con la función plot:

plot(v_mod_ok)

El objeto resultante es en realidad una lista con varios componentes, incluido el variograma empírico y el modelo de variograma ajustado. El componente “$var_model” del objeto resultante contiene el modelo real:

v_mod_ok$var_model

Modelo Nugget (Nug): se refiere a la variabilidad espacial aleatoria presente en un conjunto de datos geográficos. Representa la variabilidad que no puede ser explicada por la estructura espacial subyacente. En términos más simples, el modelo Nugget captura la variabilidad que es independiente de la distancia entre los puntos de muestreo, en el caso del departamento del meta, este valor corresponde a 61.98865.

Modelo Structure (Ste): Describe la variabilidad espacial de los datos geográficos en función de la distancia entre los puntos de muestreo. Representa la tendencia o patrón espacial presente en los datos. El modelo Structure se utiliza para estimar la correlación espacial y el alcance de la variabilidad espacial. Para el caso del departamento del Meta, estos valores corresponden a: - Correlacion espacial 152.4789 - El alcance de la variabilidad espacial 562.12296.

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

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

Nuevamente, se subdivide el atributo de valores y se renombra, en este caso se llamará Z2:

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

Las predicciones de Kriging ordinario se muestran en la siguiente figura:

m <- 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 = "Interpolación SOC del Meta por OK [%]"
    )
## 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

6. EVALUACIÓN DE RESULTADOS

6.1 EVALUACIÓN CUALITATIVA

A continuación, se presenta otra vista de las tres interpolaciones (IDW,OK y Real World):

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

EL resultado es:

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 = "Carbono orgánico del suelo [%]"
)
## 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

6.2 Validación cruzada

La validación cruzada se puede realizar con el siguiente código.

cv1 = gstat.cv(g1)
cv2 = gstat.cv(g2)
cv1 = na.omit(cv1)
cv1
## class       : SpatialPointsDataFrame 
## features    : 275 
## extent      : -74.62094, -71.08117, 1.610193, 4.696743  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 6
## names       :        var1.pred, var1.var,         observed,          residual, zscore, fold 
## min values  : 11.8164702466386,       NA, 3.80667662620544, -36.9964308054443,     NA,    1 
## max values  : 98.8344054610177,       NA, 143.717758178711,  82.4097071435825,     NA,  275

Se convierte e objeto “cv1”.

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

Ahora, se calculan los indices de precisión de predicción como el error cuadrático medio (RMSE).

Este es el valor del (RMSE) para la interpolación por el método IDW:

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

Ahora, se repite el proceso con los resultados del método OK:

Se realiza la conversión:

cv2 = st_as_sf(cv2)
sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))
## [1] 9.582831

Para el caso del departamento del Meta y de acuerdo a los datos obtenidos del ISRIC, el método de interpolación más adecuado a usar entre el kriging Ordinario (OK) y el Inverse Distance Weighting (IDW), es el primero, es decir, el Kriging Ordinario (OK) dado que es un método más avanzado que se basa en el análisis de la estructura espacial de los datos. Adicional a esto, utiliza modelos estadísticos para estimar la variabilidad espacial y las correlaciones entre los puntos de datos. Al hacerlo, puede capturar patrones más complejos y generar superficies más realistas (Mapa #5), el cual permite observar tanto el mapa del Carbono Orgánico del suelo correspondiente al departamento del Meta, comparando los métodos IDW y OK . El OK también proporciona estimaciones de incertidumbre, lo que permite evaluar la confiabilidad de las interpolaciones. Sin embargo, el OK requiere la estimación de parámetros adicionales y puede ser más computacionalmente intensivo. En resumen, si los datos muestran una estructura espacial compleja y se requiere una mayor precisión, el Kriging Ordinario (OK) es generalmente considerado como un método más preciso. Sin embargo, si los datos son más simples y se desea una implementación más sencilla, el Inverse Distance Weighting (IDW) puede ser suficiente. En última instancia, la elección del método de interpolación depende de las características de los datos y los requisitos específicos del proyecto.

7. REFERENCIAS

Para la realizacion de este cuaderno, se tuvo en cuenta el siguiente cuaderno: Lizarazo, I., 2023. Spatial interpolation of soil organic carbon. Available at https://rpubs.com/ials2un/soc_interp.

sessionInfo()
## R version 4.2.1 (2022-06-23 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.1252  LC_CTYPE=Spanish_Colombia.1252   
## [3] LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] curl_5.0.1    dplyr_1.1.1   ggplot2_3.4.2 leafem_0.2.0  leaflet_2.1.2
##  [6] automap_1.1-9 gstat_2.1-1   stars_0.6-1   abind_1.4-5   sf_1.0-12    
## [11] terra_1.7-23  sp_1.6-0     
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.10        lattice_0.20-45    FNN_1.1.3.2        png_0.1-8         
##  [5] class_7.3-20       zoo_1.8-11         digest_0.6.31      utf8_1.2.3        
##  [9] R6_2.5.1           plyr_1.8.8         evaluate_0.20      e1071_1.7-13      
## [13] highr_0.10         pillar_1.9.0       rlang_1.1.0        rstudioapi_0.14   
## [17] raster_3.6-20      jquerylib_0.1.4    rmarkdown_2.21     rgdal_1.6-5       
## [21] htmlwidgets_1.6.2  munsell_0.5.0      proxy_0.4-27       compiler_4.2.1    
## [25] xfun_0.38          pkgconfig_2.0.3    base64enc_0.1-3    htmltools_0.5.5   
## [29] tidyselect_1.2.0   tibble_3.2.1       intervals_0.15.3   codetools_0.2-18  
## [33] reshape_0.8.9      fansi_1.0.4        spacetime_1.3-0    withr_2.5.0       
## [37] grid_4.2.1         jsonlite_1.8.4     lwgeom_0.2-11      gtable_0.3.3      
## [41] lifecycle_1.0.3    DBI_1.1.3          magrittr_2.0.3     units_0.8-1       
## [45] scales_1.2.1       KernSmooth_2.23-20 cli_3.6.1          cachem_1.0.7      
## [49] farver_2.1.1       bslib_0.4.2        ellipsis_0.3.2     xts_0.13.1        
## [53] generics_0.1.3     vctrs_0.6.1        tools_4.2.1        glue_1.6.2        
## [57] crosstalk_1.2.0    parallel_4.2.1     fastmap_1.1.1      yaml_2.3.7        
## [61] colorspace_2.1-0   classInt_0.4-9     knitr_1.42         sass_0.4.5