Spatial interpolation of soil organic carbon

rm(list=ls())

library(sp)
library(terra)
library(sf)
library(stars)
library(gstat)
library(automap)
library(leaflet)
library(leafem)
library(ggplot2)
library(dplyr)

library("curl")

## Using libcurl 8.14.1 with Schannel

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

archivo <- ("C:/Users/Janus/Desktop/GB2/p7/CASANARE/soc_igh_15_30.tif")
(soc <- rast(archivo))

## class       : SpatRaster 
## size        : 917, 1413, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8145750, -7792500, 477250, 706500  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine 
## source      : soc_igh_15_30.tif 
## name        : soc_igh_15_30

soc.perc  <- soc/10

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

## class       : SpatRaster 
## size        : 934, 1489, 1  (nrow, ncol, nlyr)
## resolution  : 0.002205488, 0.002205488  (x, y)
## extent      : -73.09924, -69.81527, 4.286671, 6.346597  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : soc_igh_15_30 
## min value   :      6.961041 
## max value   :    116.768845

stars.soc = st_as_stars(geog.soc)

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

set.seed(123456)
(samples <- spatSample(geog.soc, 500, "random", as.points=TRUE))

##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 500, 1  (geometries, attributes)
##  extent      : -73.09373, -69.81638, 4.296596, 6.343289  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : soc_igh_15_30
##  type        :         <num>
##  values      :         51.31
##                        21.18
##                        11.89

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

nmuestras <- na.omit(muestras)

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

id <- seq(1,500,1)

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

sitios <- na.omit(sitios)

head(sitios)

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

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

rrr = aggregate(geog.soc, 4)

rrr

## class       : SpatRaster 
## size        : 234, 373, 1  (nrow, ncol, nlyr)
## resolution  : 0.008821954, 0.008821954  (x, y)
## extent      : -73.09924, -69.80866, 4.28226, 6.346597  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : soc_igh_15_30 
## min value   :      8.061533 
## max value   :    103.183873

values(rrr) <-1

names(rrr) <- "valor"

rrr

## class       : SpatRaster 
## size        : 234, 373, 1  (nrow, ncol, nlyr)
## resolution  : 0.008821954, 0.008821954  (x, y)
## extent      : -73.09924, -69.80866, 4.28226, 6.346597  (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)
z1 = predict(g1, stars.rrr)

## [inverse distance weighted interpolation]

z1

## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min. 1st Qu.   Median     Mean  3rd Qu.     Max.  NA's
## var1.pred  9.076501 15.6947 17.46058 21.15953 19.82693 100.4246     0
## var1.var         NA      NA       NA      NaN       NA       NA 87282
## dimension(s):
##   from  to offset     delta                       refsys x/y
## x    1 373  -73.1  0.008822 +proj=longlat +datum=WGS8... [x]
## y    1 234  6.347 -0.008822 +proj=longlat +datum=WGS8... [y]

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

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

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 = "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

m 

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

plot(v_emp_ok)

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

plot(v_mod_ok)

v_mod_ok$var_model

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

## [using ordinary kriging]

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

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 = "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 

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

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 [%]"
)

## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA

m 

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

cv1 = na.omit(cv1)

cv1

## class       : SpatialPointsDataFrame 
## features    : 479 
## extent      : -73.07388, -69.83181, 4.296596, 6.343289  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 6
## names       :        var1.pred, var1.var,         observed,          residual, zscore, fold 
## min values  : 10.0418305481674,       NA, 9.00667572021484, -17.8120498775222,     NA,    1 
## max values  : 53.1702462198962,       NA, 100.604347229004,  56.7423534589487,     NA,  479

cv1 = st_as_sf(cv1)

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

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

## [1] 6.933799

cv2 = st_as_sf(cv2)

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

## [1] 6.026548

sessionInfo()

## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
## 
## 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] curl_6.4.0     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-21     
## [11] terra_1.8-54   sp_2.2-0      
## 
## loaded via a namespace (and not attached):
##  [1] generics_0.1.4     sass_0.4.10        class_7.3-23       KernSmooth_2.23-26
##  [5] lattice_0.22-6     digest_0.6.37      magrittr_2.0.3     evaluate_1.0.3    
##  [9] grid_4.5.0         RColorBrewer_1.1-3 fastmap_1.2.0      plyr_1.8.9        
## [13] jsonlite_2.0.0     e1071_1.7-16       reshape_0.8.10     DBI_1.2.3         
## [17] crosstalk_1.2.1    scales_1.4.0       codetools_0.2-20   jquerylib_0.1.4   
## [21] cli_3.6.5          rlang_1.1.6        units_0.8-7        intervals_0.15.5  
## [25] withr_3.0.2        base64enc_0.1-3    cachem_1.1.0       yaml_2.3.10       
## [29] FNN_1.1.4.1        raster_3.6-32      tools_4.5.0        parallel_4.5.0    
## [33] spacetime_1.3-3    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    classInt_0.4-11   
## [41] htmlwidgets_1.6.4  pkgconfig_2.0.3    pillar_1.10.2      bslib_0.9.0       
## [45] gtable_0.3.6       glue_1.8.0         Rcpp_1.0.14        tidyselect_1.2.1  
## [49] tibble_3.2.1       xfun_0.52          rstudioapi_0.17.1  knitr_1.50        
## [53] farver_2.1.2       htmltools_0.5.8.1  rmarkdown_2.29     xts_0.14.1        
## [57] compiler_4.5.0

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.