Clean workspace

rm(list = ls())
library(XML)
library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
library(gdalUtilities)
## 
## Adjuntando el paquete: 'gdalUtilities'
## The following object is masked from 'package:sf':
## 
##     gdal_rasterize
library(terra)
## terra 1.8.50
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:terra':
## 
##     intersect, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(leaflet)
library(RColorBrewer)
url = "https://files.isric.org/soilgrids/latest/data/"
voi = "soc" # soil organic carbon 
depth = "15-30cm"
quantile = "mean"
# concatenate the strings
(variable = paste(url, voi, sep=""))
## [1] "https://files.isric.org/soilgrids/latest/data/soc"
(layer = paste(variable,depth,quantile, sep="_")) # layer of interest
## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean"
(vrt_layer = paste(layer, '.vrt', sep=""))
## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean.vrt"
munic = st_read("Division_muncipios_completos.shp")
## Reading layer `Division_muncipios_completos' from data source 
##   `/home/mateo/Downloads/Geom_Liza/Division_muncipios_completos.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1122 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -81.73562 ymin: -4.225938 xmax: -66.85689 ymax: 13.39473
## Geodetic CRS:  WGS 84
guainia = munic %>% 
  filter(Depto == "Guainía")
igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
guainia_igh <- st_transform(guainia, igh)
(bbox <- st_bbox(guainia_igh)) ## bounding box
##       xmin       ymin       xmax       ymax 
## -7900629.1   129741.4 -7443322.2   449947.7
ulx = bbox$xmin
uly = bbox$ymax
lrx= bbox$xmax
lry = bbox$ymin
(bb <- c(ulx, uly, lrx, lry))
##       xmin       ymax       xmax       ymin 
## -7900629.1   449947.7 -7443322.2   129741.4
sg_url="/vsicurl/https://files.isric.org/soilgrids/latest/data/"
datos = 'soc/soc_15-30cm_mean.vrt'
file = "soc_igh_15_30.tif"

# gdal_translate(paste0(sg_url,datos), file ,
#               tr=c(250,250),
#              projwin=bb,
#              projwin_srs =igh)

(guainia_soc <- terra::rast(file)/10)
## class       : SpatRaster 
## dimensions  : 1281, 1829, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -7900750, -7443500, 129750, 450000  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine 
## source(s)   : memory
## varname     : soc_igh_15_30 
## name        : soc_igh_15_30 
## min value   :           6.8 
## max value   :         175.4
#Descriptivo
terra::hist(guainia_soc, main = 'Distribución del Contenido de Carbono Orgánico del Suelo en la Guainia', xlab = 'g/kg')
## Warning: [hist] a sample of 43% of the cells was used (of which 1% was NA)

summary(guainia_soc)
## Warning: [summary] used a sample
##  soc_igh_15_30   
##  Min.   :  7.20  
##  1st Qu.: 11.70  
##  Median : 14.10  
##  Mean   : 16.03  
##  3rd Qu.: 17.90  
##  Max.   :167.80  
##  NA's   :1058
(names(guainia_soc) <-  "soc")
## [1] "soc"
valores <- values(guainia_soc, na.rm=TRUE)
orangecyan <- colorNumeric(c("orange",
                             "yellow2", 
                             "darkseagreen",
                             "cyan"),
                           valores,
                           na.color = "transparent")

El Carbono Orgánico del Suelo en el Guainia mayormente está menor a 50 g/kg

leaflet::leaflet(guainia) %>% 
  addTiles() %>% 
  setView(-68.52293954604566, 2.7601747799884104, 7.6) %>% 
   addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.5, fillOpacity = 0.10,
    popup = paste("Municipio: ", guainia$MpNombre)) %>%
  addRasterImage(guainia_soc, colors =orangecyan, opacity = 0.8)  %>%
  addLegend(pal = orangecyan, values = valores, title = "Soil Organic Carbon (SOC) [g/kg]")

Se observa que el SOC en todo el espacio dela Guainia es muy similar entre sí con valores bajos, menores a 50.

set.seed(123456)
# Conversion of SOC coordinate reference system
geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(guainia_soc, geog))
## class       : SpatRaster 
## dimensions  : 1295, 1883, 1  (nrow, ncol, nlyr)
## resolution  : 0.002221812, 0.002221812  (x, y)
## extent      : -70.96764, -66.78397, 1.165172, 4.042419  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        :        soc 
## min value   :   7.133821 
## max value   : 172.682816
# Random sampling of 2000 points
(samples <- terra::spatSample(geog.soc, 2000, "random", as.points=TRUE))
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 2000, 1  (geometries, attributes)
##  extent      : -70.9643, -66.7873, 1.168504, 4.041308  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       :   soc
##  type        : <num>
##  values      : 19.63
##                13.82
##                29.85
## Toma muestras como tipo puntos
(muestras <- sf::st_as_sf(samples))
##Omitir NA
nmuestras <- na.omit(muestras)
longit <- st_coordinates(nmuestras)[,1]
latit <- st_coordinates(nmuestras)[,2]
soc <- nmuestras$soc

summary(soc)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   8.221  11.907  14.295  16.274  18.143 135.738
length(soc)
## [1] 1952
id <- 1:1952
(sitios <- data.frame(id, longit, latit, soc))
m <- leaflet() %>%
  addTiles() %>%  
  addMarkers(lng=sitios$longit,
             lat=sitios$latit,
             popup=sitios$soc,
             clusterOptions = markerClusterOptions())
m  # Print the map

link del notebook : https://rpubs.com/ials2un/downloadSoilGrids