---

title: “cordoba spatial” autor: “Miguel Angel Sosa Camacho

introduccion

Este cuaderno ilustra dos técnicas de interpolación espacial: distancia ponderada inversa (IDW) y Kriging ordinario (OK). IDW es una técnica determinista. OK es probabilístico. Ambas técnicas se utilizan aquí para obtener una superficie continua de carbono orgánico del suelo (COS) a 15-30 cm a partir de muestras obtenidas de SoilGrids 250 m.

Limpiar el espacio de trabajo

{ setup, warning=FALSE} #rm(list=ls())

Para leer y escribir datos geoespaciales, usaremos terra (para ráster) y sf (para vector). Para la interpolación espacial, usaremos las bibliotecas gstat y automap. Para la visualización usaremos ggplot, leaflet y tmap.

Primero, necesitamos instalar las bibliotecas necesarias. Hazlo desde la consola R, no desde esta notebook.

library(terra)

## terra 1.8.29

library(sf)

## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

library(sp)
library(stars)

## Loading required package: abind

library(gstat)
library(automap)
library(RColorBrewer)
library(leaflet)
library(leafem)

Leer los datos de entrada Como se indicó en la introducción, nuestras muestras son valores de SOC a 15-30 cm de profundidad. Se obtuvieron de SoilGrids 250 m usando este cuaderno. Guardamos dichas muestras en formato geopackage, en el directorio de datos, usando soc_{depto}.gpkg como nombre de archivo.

Primero, verifiquemos que los datos estén disponibles para su uso

#list.files(path="data", pattern = "*.gpkg")

Luego, leemos los ejemplos usando sf

samples <- sf::st_read("~/Desktop/soc_cordoba.gpkg")

## Reading layer `soc_cordoba' from data source 
##   `/Users/jsc_89/Desktop/soc_cordoba.gpkg' using driver `GPKG'
## Simple feature collection with 1669 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -76.518 ymin: 7.348765 xmax: -74.71738 ymax: 9.446937
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

Como nuestra área de estudio es cordoba, tiene sentido leer también el geopaquete de municipios correspondiente

munic <- sf::st_read("~/Desktop/Cordoba.gpkg")

## Reading layer `cordoba' from data source `/Users/jsc_89/Desktop/Cordoba.gpkg' using driver `GPKG'
## Simple feature collection with 30 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.51453 ymin: 7.347087 xmax: -74.78095 ymax: 9.447748
## Geodetic CRS:  MAGNA-SIRGAS

munic <- sf::st_read(“~/Desktop/cordoba/soc_cordoba.gpkg”)

summary(samples)

##  soc_igh_15_30               geom     
##  Min.   : 10.45   POINT        :1669  
##  1st Qu.: 18.75   epsg:NA      :   0  
##  Median : 25.96   +proj=long...:   0  
##  Mean   : 30.41                       
##  3rd Qu.: 36.53                       
##  Max.   :140.25

hist(samples$soc)

Redondeemos los valores de SOC a dos dígitos

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

Ahora visualicemos las muestras para conocer su distribución espacial.

pal <- colorNumeric(
  palette = c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"), 
  domain = samples$soc
)

ahora se de asegurar que el objecto Munic contiene datos para el mapa

munic_points <- st_centroid(munic)

## Warning: st_centroid assumes attributes are constant over geometries

  leaflet() %>%
    addCircleMarkers(  
      data = munic_points,
      radius = 2,
      color = "gray",
      fillOpacity = 0.2
    ) %>%
    addCircleMarkers(
      data = samples,
      radius = 1.5, 
      label = ~soc,
      color = ~pal(soc),
      fillOpacity = 1,
      stroke = FALSE
    ) %>%
    addLegend(
      pal = pal,
      values = samples$soc,  # Se usa directamente samples$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'

5. Interpolación 5.1 Creación del objeto gstat Para interpolar, primero necesitamos crear un objeto de clase gstat, usando una función del mismo nombre: gstat.

Un objeto gstat contiene toda la información necesaria para realizar la interpolación espacial, a saber:

La definición del modelo.

Los datos de calibración

Según sus argumentos, la función gstat “entiende” qué tipo de modelo de interpolación queremos usar:

Sin modelo de variograma → Distancia de ponderación inversa (IDW)

Modelo de variograma, sin covariables → Kriging ordinario (OK)

Vamos a utilizar tres parámetros de la función gstat:

fórmula: la “fórmula” de predicción que especifica las variables dependientes e independientes (también conocidas como covariables)

datos: los datos de calibración (también conocidos como datos del tren)

modelo: el modelo de variograma

5.2 interpolación IDW Para interpolar SOC usando el método IDW creamos el siguiente objeto gstat, especificando solo la fórmula y los datos

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

Ahora que nuestro modelo de interpolación g1 está definido, podemos usar la función de predicción para interpolar realmente, es decir, estimar los valores de SOC.

La función de predicción acepta:

Un objeto ráster: estrellas, como dem Un modelo: objeto gstat, como g1 El ráster sirve para dos propósitos:

Especificar las ubicaciones donde queremos hacer predicciones (en todos los métodos) Especificación de valores de covariables (solo en Universal Kriging) Usemos la biblioteca terra para crear un objeto ráster con valores de celda iguales a 1. La extensión de este ráster debe cubrir nuestra área de interés.

(rrr <- terra::rast(xmin=-74.55524, xmax=-72.38985, ymin=5.708308, ymax=8.145474, nrows=370, ncols = 329, vals = 1, crs = "epsg:4326"))

## class       : SpatRaster 
## dimensions  : 370, 329, 1  (nrow, ncol, nlyr)
## resolution  : 0.006581733, 0.006586935  (x, y)
## extent      : -74.55524, -72.38985, 5.708308, 8.145474  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : lyr.1 
## min value   :     1 
## max value   :     1

Ahora necesitamos convertir este SpatRaster en un objeto de stars

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

La siguiente expresión interpola los valores de SOC según el modelo definido en 𝑔1 y el ráster definido en 𝑠𝑡𝑎𝑟𝑠.𝑟𝑟𝑟

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

## [inverse distance weighted interpolation]

¿Qué es z1? acontinuación vemos

z1

## stars object with 2 dimensions and 2 attributes
## attribute(s):
##              Min.  1st Qu.   Median     Mean  3rd Qu.     Max.   NA's
## var1.pred  27.496 29.45154 29.94631 29.91671 30.36229 31.39016      0
## var1.var       NA       NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -74.56  0.006582 WGS 84 [x]
## y    1 370  8.145 -0.006587 WGS 84 [y]

Podemos subconjuntar sólo el primer atributo y cambiarle el nombre, y así podemos manejar los valores

# dealing with NA values
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  27.496 29.45154 29.94631 29.91671 30.36229 31.39016      0
## var1.var       NA       NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -74.56  0.006582 WGS 84 [x]
## y    1 370  8.145 -0.006587 WGS 84 [y]

names(z1) = "soc"

Para trazar la superficie interpolada, necesitamos una paleta de colores

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

Ahora, visualicemos el resultado de la interpolación

combined_domain <- range(c(z1$soc, samples$soc), na.rm = TRUE)

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = combined_domain)
    ) %>%
   addCircleMarkers(
    data = samples,
    radius = 1.5, 
    label = ~soc,
    color = ~paleta(soc),  # Ensure paleta is adjusted for the new range
    fillOpacity = 1,
    stroke = FALSE
  ) %>%
    addLegend("bottomright", pal = paleta, values = c(z1$soc, samples$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

5.3 Interpolación correcta Los métodos de 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 ponderaciones en consecuencia al realizar predicciones.

Como primer paso, podemos calcular y examinar el variograma empírico utilizando la función de variograma.

La función requiere dos argumentos:

fórmula: especifica la variable dependiente y las covariables, como en gstat

datos: la capa de puntos con la variable dependiente y covariables como atributos de puntos

La siguiente expresión calcula el variograma empírico de muestras, sin covariables

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

Tracemos el variograma empírico

plot(v_emp_ok)

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

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

La función elige el tipo de modelo que mejor se adapta y también ajusta sus parámetros. Puedes usar 𝑠ℎ𝑜𝑤.𝑣𝑔𝑚𝑠() para mostrar varios tipos de modelos de variograma.

Tenga en cuenta que la función autofitVariogram no funciona en objetos sf, por lo que convertimos el objeto de muestra en un SpatialPointsDataFrame (paquete sp).

El modelo ajustado se puede trazar con plot

plot(v_mod_ok)

Esta gráfica representa la distancia entre puntos en el eje X y la semivarianza en el eje Y, en este caso que nuestra curva es creciente interpretamos que los puntos separados son valores mas diferentes a medida a que la distancia crece.

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

##   model     psill    range kappa
## 1   Nug  21.35912  0.00000   0.0
## 2   Ste 305.43597 89.14884   0.3

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

#esto toma tiempo según la computadora
g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)
z2= predict(g2, stars.rrr)

## [using ordinary kriging]

Nuevamente, crearemos un subconjunto del atributo de valores predichos y le cambiaremos el nombre

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

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

Evaluacion cualitativa

colores <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), domain = 10:100, na.color = "transparent")
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 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

Validación Cruzada

### 
## Esta parte del código solo debe correrse desde la consola, debido al tiempo de execución, que afecta el knitting process. Dependiendo de la máquina puede tomar hasta 50 minutos.
cv2 = gstat.cv(g2)

Validación cruzada Hemos estimado las superficies de SOC utilizando dos métodos diferentes: IDW y Kriging ordinario. Aunque es útil examinar y comparar los resultados gráficamente, también necesitamos una forma objetiva de evaluar la precisión de la interpolación. De esa manera, podemos elegir objetivamente el método más preciso entre los métodos de interpolación disponibles.

En la validación cruzada Leave-One-Out nosotros:

Saque un punto de los datos de calibración. Haz una predicción para ese punto. Repita para todos los puntos. Al final, lo que obtenemos es una tabla con un valor observado y un valor predicho para todos los puntos. Podemos ejecutar la validación cruzada Leave-One-Out usando la función gstat.cv, que acepta un objeto gstat.

Anote el RSME obtenido para obtener resultados de interpolación correctos.

Luego, repita el proceso con los resultados de la interpolación IDW.

Una vez que calcule ambos valores RMSE, compárelos y analice sus diferencias. ¿Cuáles son los métodos de interpolación más confiables? Explicar.

# ejecutar en la consola
#cv2 = st_as_sf(cv2)

ahora graficar los residuales

#ejecutar en la consola
#sp::bubble(as(cv2[, "residual"], "Spatial"))

# This is the RMSE value for the OK interpolation
# ejecutar en la consola
#sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))

Guardar los resultados de la interpolacion

write_stars(
  z1, dsn = "~/Desktop/IDW_soc_cordoba.tif", layer = 1
)

write_stars(
  z2, dsn = "~/Desktop/OK_soc_stder.tif", layer = 1
)

Conclusiones

Con este ejercicio aprendí sobre nuevos tipos de objetos y funciones utiles para el uso de tecnicas de interpolacion espacial, por un lado La distancia invertida ponderada (IDW) y por otro el Kriging ordinario (OK), esto para generar un superficia continua de SOC. Esto también requirió que buscará mayor información acerca de las técnicas y los retos técnicos de configurar R y Rstudio.