Configuración

Limpiar zona de trabajo

rm(list=ls())
# Borrar los mensajes no necesarios
options(showPackageStartupMessages = FALSE)

Instalar librerias

library(terra)
## terra 1.8.5
library(sf)
## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE
library(leaflet)
library(XML)
library(gdalUtilities)
## 
## Adjuntando el paquete: 'gdalUtilities'
## The following object is masked from 'package:sf':
## 
##     gdal_rasterize
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(RColorBrewer)
url = "https://files.isric.org/soilgrids/latest/data/"
voi = "soc"  
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="_"))
## [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"
(Narino <- st_read("C:/Users/danna/Desktop/Geomática Básica 2024-2/Proyecto5/GB_sem_12/Datos/NARINO_.gpkg"))
## Reading layer `municipios_colombia' from data source 
##   `C:\Users\danna\Desktop\Geomática Básica 2024-2\Proyecto5\GB_sem_12\Datos\NARINO_.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 64 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -79.01021 ymin: 0.3613481 xmax: -76.83368 ymax: 2.683898
## Geodetic CRS:  MAGNA-SIRGAS
igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
Narino_igh <- st_transform(Narino, igh)
(bbox <- st_bbox(Narino_igh))
##        xmin        ymin        xmax        ymax 
## -8796326.48    40225.08 -8554103.93   298770.17
ulx = bbox$xmin
uly = bbox$ymax
lrx= bbox$xmax
lry = bbox$ymin
(bb <- c(ulx, uly, lrx, lry))
##        xmin        ymax        xmax        ymin 
## -8796326.48   298770.17 -8554103.93    40225.08
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)
(Narino_soc <- terra::rast(file)/10)
## class       : SpatRaster 
## dimensions  : 1034, 969, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8796500, -8554250, 40500, 299000  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine 
## source(s)   : memory
## varname     : soc_igh_15_30 
## name        : soc_igh_15_30 
## min value   :          13.1 
## max value   :         262.3
terra::hist(Narino_soc)
## Warning: [hist] a sample of 100% of the cells was used (of which 10% was NA)

summary(Narino_soc)
## Warning: [summary] used a sample
##  soc_igh_15_30   
##  Min.   : 14.90  
##  1st Qu.: 41.50  
##  Median : 62.30  
##  Mean   : 79.73  
##  3rd Qu.:110.10  
##  Max.   :253.40  
##  NA's   :10445
(names(Narino_soc) <-  "soc")
## [1] "soc"
valores <- values(Narino_soc, na.rm=TRUE)
orangecyan <- colorNumeric(c("orange","yellow2", "cyan",  "darkseagreen" ), valores,
  na.color = "transparent")
leaflet::leaflet(Narino) %>% 
  addTiles() %>% 
  setView(-78, 1.55, 8) %>% 
   addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.5, fillOpacity = 0.10,
    popup = paste("Municipio: ", Narino$MPIO_CNMBR)) %>%
  addRasterImage(Narino_soc, colors ="Spectral", opacity = 0.8)  %>%
  addLegend(pal = orangecyan, values = valores, title = "Carbono Organico del Suelo (SOC) [g/kg]")
## 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'
set.seed(12345)
geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(Narino_soc, geog))
## class       : SpatRaster 
## dimensions  : 1039, 985, 1  (nrow, ncol, nlyr)
## resolution  : 0.00223491, 0.00223491  (x, y)
## extent      : -79.01988, -76.81849, 0.3638916, 2.685963  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        :       soc 
## min value   :  14.71794 
## max value   : 261.19443
(samples <- terra::spatSample(geog.soc, 2000, "random", as.points=TRUE))
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 2000, 1  (geometries, attributes)
##  extent      : -79.01876, -76.81961, 0.365009, 2.684845  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       :   soc
##  type        : <num>
##  values      : 68.51
##                   NA
##                41.58
(muestras <- sf::st_as_sf(samples))
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. 
##   18.24   41.87   61.58   78.12  103.89  215.61
length(soc)
## [1] 1758
id <- seq(1,1758,1)
(sitios <- data.frame(id, longit, latit, soc))
m <- leaflet() %>%
  addTiles() %>%  
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions())
m 
sf::st_write(nmuestras, "Datos/soc_narino.gpkg", driver = "GPKG",
             append=FALSE)
## Deleting layer `soc_narino' using driver `GPKG'
## Writing layer `soc_narino' to data source `Datos/soc_narino.gpkg' using driver `GPKG'
## Writing 1758 features with 1 fields and geometry type Point.

Referencia

Lizarazo, I. 2025. How to download soil data from SoilGrids. Available at: https://rpubs.com/ials2un/downloadSoilGrids

sessionInfo()
## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## 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] RColorBrewer_1.1-3  dplyr_1.1.4         gdalUtilities_1.2.5
## [4] XML_3.99-0.18       leaflet_2.2.2       sf_1.0-19          
## [7] terra_1.8-5        
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.9         utf8_1.2.4         generics_0.1.3     class_7.3-22      
##  [5] KernSmooth_2.23-24 digest_0.6.37      magrittr_2.0.3     evaluate_1.0.1    
##  [9] grid_4.4.2         fastmap_1.2.0      jsonlite_1.8.9     e1071_1.7-16      
## [13] DBI_1.2.3          fansi_1.0.6        crosstalk_1.2.1    scales_1.3.0      
## [17] codetools_0.2-20   jquerylib_0.1.4    cli_3.6.3          rlang_1.1.4       
## [21] units_0.8-5        munsell_0.5.1      cachem_1.1.0       yaml_2.3.10       
## [25] tools_4.4.2        colorspace_2.1-1   vctrs_0.6.5        R6_2.5.1          
## [29] png_0.1-8          proxy_0.4-27       lifecycle_1.0.4    classInt_0.4-10   
## [33] htmlwidgets_1.6.4  pkgconfig_2.0.3    pillar_1.9.0       bslib_0.8.0       
## [37] glue_1.8.0         Rcpp_1.0.13-1      xfun_0.49          tibble_3.2.1      
## [41] tidyselect_1.2.1   rstudioapi_0.17.1  knitr_1.49         farver_2.1.2      
## [45] htmltools_0.5.8.1  rmarkdown_2.29     compiler_4.4.2