spatial interpolation

Indice

Objetivo

Obtener una superficie continua de carbono orgánico del suelo (SOC) a 15-30 cm a partir de muestras obtenidas de SoilGrids 250 m.

Preparación del Entorno

Limpiar espacio de trabajo
Cargar librerías

Carga de Datos

Cargar capas de puntos y polígonos del anterior notebook

Análisis Exploratorio de Datos

Gráficos y estadísticos descriptivos
- Rampa de colores para mapear puntos de muestreo de Carbono orgánico
- Mapeo de puntos de muestreo

Interpolación

Interpolación Inversa a la Distancia (IDW)
- Crear objeto gstat para modelar SOC
- Crear raster para el departamento de Tolima
- Convertir raster a stars y realizar cálculo de pronósticos
- Visualización de resultados IDW
Kriging Ordinario (OK)
- Variograma empírico
- Ajuste del variograma
- Ejecución del Kriging Ordinario
- Visualización de resultados OK

Evaluación de Resultados

Comparación visual de IDW y OK
Cross-Validation
Validación cruzada para IDW
Validación cruzada para OK
Cálculo del RMSE para IDW y OK

Exportación de Resultados

Guardar resultados de IDW y OK en archivos

8.Conclusiones

Limpiar espacio de trabajo

rm(list = ls())

Cargar librerias

library(terra)

## terra 1.8.15

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)

cargar capas de puntos y polígonos del anterior notebook

samples <- sf::st_read("SOC_tolima.gpkg")

## Reading layer `SOC_tolima' from data source 
##   `C:\Users\yusel\Downloads\SPATIAL\SPATIAL\SOC_tolima.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1898 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -76.10565 ymin: 2.876745 xmax: -74.4456 ymax: 5.319164
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

munic <- sf::st_read("TOLIMA_MAGNA_SIRGAS.shp")

## Reading layer `TOLIMA_MAGNA_SIRGAS' from data source 
##   `C:\Users\yusel\Downloads\SPATIAL\SPATIAL\TOLIMA_MAGNA_SIRGAS.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 47 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.10574 ymin: 2.871081 xmax: -74.47482 ymax: 5.319342
## Geodetic CRS:  WGS 84

Graficos y estadísticos descriptivos

summary(samples)

##       soc                   geom     
##  Min.   :12.44   POINT        :1898  
##  1st Qu.:25.27   epsg:NA      :   0  
##  Median :38.39   +proj=long...:   0  
##  Mean   :42.07                       
##  3rd Qu.:57.74                       
##  Max.   :96.99

hist(samples$soc)

Rampa de colores para mapear puntos de muestreo de Carbono orgánico

samples$soc = round(samples$soc,2)
pal <- colorNumeric(
  c("#253a40", "#425751", "#be9d63", "#b88750", "#9d7751"),
  domain = samples$soc,reverse = T
  )

Mapeo de puntos de muestreo

leaflet() %>%
  addPolygons(
    data = munic,
    color = "green",
    opacity = 0.8,
    weight = 0.5,
    fillOpacity = 0.1) %>%
 addCircleMarkers(
    data = samples,
    radius= 1.3, 
    label = ~soc,
    color = ~pal(soc),
    fillOpacity = 0.7,
    stroke = F
  ) %>%
  addLegend(
    data = samples,
    pal = pal,
    values = ~soc,
    position = "bottomright",
    title = "SOC:",
    opacity = 0.9) %>%
  addProviderTiles("OpenStreetMap")

5. Interpolation

Crear el objeto gstat que se quiere modelar el Carbono organico en función de una medida constante sin covairables, de los valores cercanos.

IDW

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

raster para el departamento de Tolima

(rrr <- terra::rast(xmin=-76.10565, xmax=-74.4456, ymin=2.876745, ymax=5.319164, nrows=270, ncols = 329, vals = 1, crs = "epsg:4326"))

## class       : SpatRaster 
## dimensions  : 270, 329, 1  (nrow, ncol, nlyr)
## resolution  : 0.005045745, 0.009045996  (x, y)
## extent      : -76.10565, -74.4456, 2.876745, 5.319164  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : lyr.1 
## min value   :     1 
## max value   :     1

convertir el raster a stars y realizamos el calculo de los pronósticos mediante la interpolación IDW

stars.rrr <- stars::st_as_stars(rrr)
z1 = predict(g1, stars.rrr)

## [inverse distance weighted interpolation]

z1 contiene las predicciones y la varianza del SOC

z1

## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean  3rd Qu.     Max.  NA's
## var1.pred  13.23168 30.72785 40.51173 42.00109 52.85222 95.50745     0
## var1.var         NA       NA       NA      NaN       NA       NA 88830
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -76.11  0.005046 WGS 84 [x]
## y    1 270  5.319 -0.009046 WGS 84 [y]

Tomar solo los valores ajustados

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

## [1] 13.23168

max(z1[[1]])

## [1] 95.50745

paleta <- colorNumeric(palette = c("#253a40", "#425751", "#be9d63", "#b88750", "#9d7751"), domain = 10:100, na.color = "transparent", reverse = T)

Se observan los valores del carbono orgánico en espacios bien definidos

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("#9d7751","#b88750","#be9d63","#425751", "#253a40"), domain = 10:100)
    ) %>%
   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 [%]"
    )
m

Kriging Ordinal

Variograma empírico para ver la autocorrelación espacial del carbono orgánico

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

Ajustar el variograma

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

plot(v_mod_ok)

El modelo ajustado es un Spherical que asume que la correlación espacial del carbono orgánico disminuye con la distancia. El nugget 29 sugiere que hay una cantidad significativa de variabilidad a pequeña escala. La variabilidad total capturada por el modelo es 431, puede ser clasificada como “moderada”.

v_mod_ok$var_model

##   model     psill    range
## 1   Nug  29.05621  0.00000
## 2   Sph 401.78290 96.33817

Ejecutando el Kriging Ordinal

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

## [using ordinary kriging]

Tomar los valores predichos

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("#253a40", "#425751", "#be9d63", "#b88750", "#9d7751"),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 write_stars.stars(x, dsn = fl): all but first attribute are ignored

Se observa un comportamiento un poco diferente al del IDW, las franjas donde el SOC es mayor visualmente ocupa una menor área a la presentada en la anterior salida.

6. Evaluar los resultados

colores <- colorOptions(palette = c("#9d7751","#b88750","#be9d63","#425751", "#253a40"), 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

Se logra distinguir un tipo de suavizado y mayor detalle de la interpolación OK, y debe ser porque este método considera las distancias entre los puntos y la correlación espacial entre los mismos, mientras que la IDW solo la distancia.

6.2 Cross-Validation

Eliminando y estimando dato a dato Primero con el IDW

#cv1 = gstat.cv(g1)

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

# run from the console
#sp::bubble(as(cv1[, "OK residuals"], "Spatial"))

Luego con el OK

#cv2 = gstat.cv(g2)

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

# run from the console
#sp::bubble(as(cv2[, "OK residuals"], "Spatial"))

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

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

# write_stars(
#   z1, dsn = "IDW_soc_stder.tif", layer = 1
# )

# write_stars(
#   z2, dsn = "OK_soc_stder.tif", layer = 1
# )

8. Conclusions

El IDW con un RMSE de 9.56 indica que en promedio las predicciones de SOC realizadas con IDW se desvían esa 9.56 unidades de los valores observados. En el OK con un RMSE de 7.98 indica que en promedio las predicciones de SOC realizadas con el método de Kriging Ordinario se desvían 7.98 unidades de los valores observados. En conclusión el método OK muestra que las predicciones son más precisas. Sin embargo, también los recursos computacionales requeridos para este desarrollo pueden ser muy altos y en algún experimento muy grande habria que definir una metodología que me permita tomar decisiones con los recursos disponibles.

El notebook presentó herramientas muy útiles para conocer el comportamiento de variables del suelo a lo largo del departamento de Tolima que pueden ser muy útiles para la toma de decisiones agrícolas como el establecimiento de cultivos que presenten propiedades que incrementen rendimientos en presencia de altos contenidos de Carbono Orgánico en el Suelo.