##1. Introducción En este cuaderno veremos dos formas de interpolación espacial: Inverse Distance Weighted (IDW), y Ordinary Kriging (OK), utilizaremos estas técinas para obtener una superficie contínua de Arcilla (clay) en una profundidad de 15-30 cm a partir de las muestras obtenidas en SoilGrids250m.

rm(list=ls())
library(terra)
## terra 1.9.27
library(sf)
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
library(sp)
## Warning: package 'sp' was built under R version 4.6.1
library(stars)
## Warning: package 'stars' was built under R version 4.6.1
## Cargando paquete requerido: abind
library(gstat)
## Warning: package 'gstat' was built under R version 4.6.1
library(automap)
## Warning: package 'automap' was built under R version 4.6.1
library(RColorBrewer)
library(leaflet)
library(leafem)
## Warning: package 'leafem' was built under R version 4.6.1
list.files()
##  [1] "clay_igh_15_30.tif"            "clay_stder.gpkg"              
##  [3] "IDW_clay_stder.tif"            "interpolacionespacial.html"   
##  [5] "interpolacionespacial.nb.html" "interpolacionespacial.Rmd"    
##  [7] "interpolacionespacial_files"   "muestreo.gpkg"                
##  [9] "municipios_col.gpkg"           "MunicipiosSantander.gpkg"     
## [11] "OK_clay_stder.tif"             "puntos.gpkg"                  
## [13] "Soilgrids.nb.html"             "Soilgrids.Rmd"
samples <- sf::st_read("clay_stder.gpkg")
## Reading layer `clay_stder' from data source 
##   `C:\Users\isabe\OneDrive\Documentos\GB2\Proyecto7\clay_stder.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 2000 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.99655 ymin: 6.006384 xmax: -72.00211 ymax: 7.998851
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs
munic <- sf::st_read("MunicipiosSantander.gpkg")
## Reading layer `col_adm2' from data source 
##   `C:\Users\isabe\OneDrive\Documentos\GB2\Proyecto7\MunicipiosSantander.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 87 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.55179 ymin: 5.705701 xmax: -72.5161 ymax: 8.120555
## Geodetic CRS:  WGS 84
summary(samples)
##  clay_igh_15_30             geom     
##  Min.   :  0.0   POINT        :2000  
##  1st Qu.:343.5   epsg:NA      :   0  
##  Median :370.7   +proj=long...:   0  
##  Mean   :358.7                       
##  3rd Qu.:393.5                       
##  Max.   :488.0
hist(samples$clay)

samples$clay = round(samples$clay,2)
pal <- colorNumeric(c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"), domain = samples$clay)
leaflet() %>%
  addPolygons(
    data = munic,
    color = "gray",
    opacity = 1,
    weight = 1,
    fillOpacity = 0.2) %>%
 addCircleMarkers(
    data = samples,
    radius= 1.5, 
    label = ~clay,
    color = ~pal(clay),
    fillOpacity = 1,
    stroke = F
  ) %>%
  addLegend(
    data = samples,
    pal = pal,
    values = ~clay,
    position = "bottomleft",
    title = "CLAY:",
    opacity = 0.9) %>%
  addProviderTiles("OpenStreetMap")
g1 = gstat(formula = clay ~1, data = samples)
(rrr<- terra::rast(xmin = -74.55179, xmax = -72.00211, ymin = 5.705701, ymax = 8.120555, nrows = 370, ncols = 329, vals = 1, crs = "epsg:4326"))
## class       : SpatRaster
## size        : 370, 329, 1  (nrow, ncol, nlyr)
## resolution  : 0.007749787, 0.006526632  (x, y)
## extent      : -74.55179, -72.00211, 5.705701, 8.120555  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s)   : memory
## name        : lyr.1
## min value   :     1
## max value   :     1
stars.rrr <- stars::st_as_stars(rrr)
z1 = predict(g1, stars.rrr)
## [inverse distance weighted interpolation]
z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
## Warning in (function (..., deparse.level = 1) : number of columns of result is
## not a multiple of vector length (arg 1)
##                Min.  1st Qu.   Median    Mean  3rd Qu.     Max.    NAs     NA's
## var1.pred  1.334276 341.5757 365.7062 356.876 376.1191 484.6013      0 1.334276
## var1.var         NA       NA       NA     NaN       NA       NA 121730 0.000000
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -74.55   0.00775 WGS 84 [x]
## y    1 370  8.121 -0.006527 WGS 84 [y]
library("terra")
library("leaflet")
r_final <- terra::rast(z1[1])

pal_raster <- colorNumeric(
palette = "YlGn",
domain = terra::values(r_final),
na.color = "transparent"
)

leaflet() %>%
addProviderTiles("OpenStreetMap") %>%
addRasterImage(
r_final,
colors = pal_raster, opacity = 0.7 ) %>% addPolygons( data = munic, color = "gray", weight = 1, fillOpacity = 0.2 ) %>% addCircleMarkers( data = samples, radius = 1.5, label = ~clay, color = ~pal(clay), fillOpacity = 1, stroke = FALSE ) %>% addLegend( pal = pal_raster, values = terra::values(r_final), position = "bottomleft", title = "Soil Clay [%]" )
a <- which(is.na(z1[[1]]))
z1[[1]][a] = 0.0
z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
## Warning in (function (..., deparse.level = 1) : number of columns of result is
## not a multiple of vector length (arg 1)
##                Min.  1st Qu.   Median    Mean  3rd Qu.     Max.    NAs     NA's
## var1.pred  1.334276 341.5757 365.7062 356.876 376.1191 484.6013      0 1.334276
## var1.var         NA       NA       NA     NaN       NA       NA 121730 0.000000
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -74.55   0.00775 WGS 84 [x]
## y    1 370  8.121 -0.006527 WGS 84 [y]
names(z1) = "clay"
v_emp_ok = variogram(clay ~1, data= samples)
plot(v_emp_ok)

v_mod_ok = automap::autofitVariogram(clay ~1, as (samples, "Spatial"))
## Warning in automap::autofitVariogram(clay ~ 1, as(samples, "Spatial")): Some models where removed for being either NULL or having a negative sill/range/nugget, 
##  set verbose == TRUE for more information
plot(v_mod_ok)

v_mod_ok$var_model
g2 = gstat(formula = clay ~1, model = v_mod_ok$var_model, data = samples, nmax = 30)
z2 = predict(g2, stars.rrr)
## [using ordinary kriging]
a<- which (is.na(z2[[1]]))
z2[[1]][a]=0.0
names(z2) = "clay"
paleta <- colorNumeric(
  palette = c("orange", "yellow", "cyan", "green"),
  domain = c(0, 100)
)
m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
    z2,
    opacity = 0.7,
    colorOptions = colorOptions(
      palette = c("orange", "yellow", "cyan", "green"),
      domain = range(z2$clay, na.rm = TRUE)
    )
  ) %>%
  addCircleMarkers(
    data = samples,
    radius = 1.5,
    label = ~clay,
    color = ~paleta(clay),
    fillOpacity = 1,
    stroke = FALSE
  ) %>%
  addLegend(
    "bottomright",
    pal = paleta,
    values = c(0,100),
    title = "OK CLAY interpolation [%]"
  )
## Warning in paleta(clay): Some values were outside the color scale and will be
## treated as NA
## Warning in paleta(clay): 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
colores <- colorOptions(
  palette = c("orange", "yellow", "cyan", "green"),
  domain = c(0,500),
  na.color = "transparent"
)
m <- leaflet() %>%
  addTiles()%>%
  addGeoRaster(z1, 
               colorOptions = colores,
               opacity = 0.8,
               group = "IDW"
               )%>%
  addGeoRaster(z2,
               colorOptions = colores,
               opacity = 0.8,
               group = "OK"
               )%>%
  addLayersControl(
    overlayGroups = c("IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  addLegend(
    "bottomright",
    pal = paleta,
    values =c(0,500),
    title = "Soil Clay [%]"
  )
## Warning in pal(c(r[1], cuts, r[2])): 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
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m