Cómo descargar datos de suelo de SoilGrids
Gabriel Andres Gutierrez Gomez
09/11/2023
1.Introducción.
SoilGrids es un sistema para el mapeo digital global de suelos que utiliza información global del perfil del suelo y datos covariables para modelar la distribución espacial de las propiedades del suelo en todo el mundo. SoilGrids es una colección de mapas de propiedades del suelo para el mundo producidos mediante aprendizaje automático con una resolución de 250 m. Se puede acceder a SoilGrids a través de los siguientes servicios proporcionados por el Centro Internacional de Información y Referencia de Suelos - (ISRIC) [https://www.isric.org/explore/soilgrids]:
Web Mapping Service (WMS): acceso para visualización y resumen de datos.
Servicio de cobertura web (WCS): la mejor manera de obtener un subconjunto de un mapa y utilizar SoilGrids como entrada para otros canales de modelado.
Creación y control de versiones distribuidas en la web (WebDAV): descargue los mapas globales completos en formato VRT.
Usaremos el servicio WebDAV en este cuaderno. En caso de duda, por favor revise en https://www.isric.org/explore/soilgrids/soilgrids-access
2.Configuración
Se limpia el espacio de trabajo
rm(list = ls())Se instalan las librerias desde la consola:
#install.packages('XML')
#install.packages('rgdal')
#install.packages('gdalUtils')
#install.packages('sf')
#install.packages('dplyr')
#install_github("envirometrix/landmap")+
#install.packages(leaflet)
#install.packages(mapview)
##install.packages("devtools")
##devtools:::install_github('gearslaboratory/gdalUtils')Ahora se cargan las respectivas librerias:
library(XML)## Warning: package 'XML' was built under R version 4.3.2library(raster)## Loading required package: sp## Warning: package 'sp' was built under R version 4.3.2library(stars)## Loading required package: abind## Loading required package: sf## Warning: package 'sf' was built under R version 4.3.2## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUElibrary(sf)
library(dplyr)## Warning: package 'dplyr' was built under R version 4.3.2##
## Attaching package: 'dplyr'## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(RColorBrewer)
library(mapview)## Warning: package 'mapview' was built under R version 4.3.2library(tmap)## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')library(devtools)## Warning: package 'devtools' was built under R version 4.3.2## Loading required package: usethis## Warning: package 'usethis' was built under R version 4.3.2library(devtools)
install_github("JoshOBrien/gdalUtilities")## Downloading GitHub repo JoshOBrien/gdalUtilities@HEAD##
## ── R CMD build ─────────────────────────────────────────────────────────────────
##
checking for file 'C:\Users\gagug\AppData\Local\Temp\RtmpCSmkEw\remotes3d847e704fcf\JoshOBrien-gdalUtilities-3d87b46/DESCRIPTION' ...
✔ checking for file 'C:\Users\gagug\AppData\Local\Temp\RtmpCSmkEw\remotes3d847e704fcf\JoshOBrien-gdalUtilities-3d87b46/DESCRIPTION' (734ms)
##
─ preparing 'gdalUtilities': (837ms)
## checking DESCRIPTION meta-information ...
✔ checking DESCRIPTION meta-information
##
─ installing the package to process help pages
##
─ saving partial Rd database (9.8s)
##
─ checking for LF line-endings in source and make files and shell scripts
##
─ checking for empty or unneeded directories
##
─ building 'gdalUtilities_1.2.5.tar.gz'
##
## ## Installing package into 'C:/Users/gagug/AppData/Local/R/win-library/4.3'
## (as 'lib' is unspecified)## Warning in i.p(...): installation of package
## 'C:/Users/gagug/AppData/Local/Temp/RtmpCSmkEw/file3d841c052999/gdalUtilities_1.2.5.tar.gz'
## had non-zero exit statuslibrary(gdalUtilities)## Warning: package 'gdalUtilities' was built under R version 4.3.2##
## Attaching package: 'gdalUtilities'## The following object is masked from 'package:sf':
##
## gdal_rasterizeEstablecer variables de interés:
Primero, definimos la url correspondiente al sitio webDAD de ISRIC:
url = "https://files.isric.org/soilgrids/latest/data/" # Ruta a los datos webDAV.Luego, creamos algunos objetos indicando lo que necesitamos de ISRIC:
# basic strings
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"Ahora definimos la propiedad del suelo que queremos descargar:
#
(layer = paste(variable,depth,quantile, sep="_")) #Capa de interes## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean"Entonces, la especificación VRT está completa:
(vrt_layer = paste(layer, '.vrt', sep=""))## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean.vrt"3. Descargue un Tiff para una región de interés (ROI)
Leamos un shapefile de un departamento previamente descargado del DANE. Estaremos aquí la biblioteca sf:
list.files()## [1] "DescargaDatosSoilGrids.html" "DescargaDatosSoilGrids.nb.html"
## [3] "DescargaDatosSoilGrids.Rmd" "imagen1.png"
## [5] "rsconnect" "Shapefile boyacá.cpg"
## [7] "Shapefile boyacá.dbf" "Shapefile boyacá.prj"
## [9] "Shapefile boyacá.qmd" "Shapefile boyacá.shp"
## [11] "Shapefile boyacá.shx" "soc_igh_15_30.tif"(boy = st_read("Shapefile boyacá.shp"))## Reading layer `Shapefile boyacá' from data source
## `C:\Users\gagug\OneDrive\Escritorio\GEOMATICA-20230824T155303Z-001\GEOMATICA\GB2\Cuadernos SoilGrids\Shapefile boyacá.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 246 features and 8 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## Geodetic CRS: MAGNA-SIRGASTenga en cuenta el CRS de boy. Como las capas ISRIC usan la proyección Homolosine, necesitamos reproyectar nuestra capa usando la biblioteca sf:
igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
boy_igh <- st_transform(boy, igh)Consigamos el cuadro delimitador de nuestra área de interés:
(bbox <- st_bbox(boy_igh))## xmin ymin xmax ymax
## -8325947 518214 -8032661 785421Ahora, usemos los datos de bbox para definir los límites del cuadro de límites tal como los usa la biblioteca GDA. Por cierto, esta es una de las partes más complicadas del uso de GDAL.
## ul significa arriba a la izq.
## lr significa arriba a la drch.
ulx = bbox$xmin
uly = bbox$ymax
lrx= bbox$xmax
lry = bbox$ymin
(bb <- c(ulx, uly, lrx, lry))## xmin ymax xmax ymin
## -8325947 785421 -8032661 518214Ahora podemos usar la función gdal_translate para descargar la capa vrt. Primero, definamos dónde guardar los datos:
sg_url="/vsicurl/https://files.isric.org/soilgrids/latest/data/"
datos = 'soc/soc_15-30cm_mean.vrt'
file = "soc_igh_15_30.tif"Tenga en cuenta que la siguiente tarea puede tardar varios minutos:
gdal_translate(paste0(sg_url,datos), file ,
tr=c(250,250),
projwin=bb,
projwin_srs =igh)Recordemos cuáles son las unidades de la capa SOC: 
Al dividir los valores de las predicciones por los valores de la columna Factor de conversión, el usuario puede obtener las unidades más familiares en la columna Unidades convencionales.
(boy_soc <- raster(file)/10)## class : RasterLayer
## dimensions : 1069, 1173, 1253937 (nrow, ncol, ncell)
## resolution : 250, 250 (x, y)
## extent : -8326000, -8032750, 518250, 785500 (xmin, xmax, ymin, ymax)
## crs : +proj=igh +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## source : memory
## names : soc_igh_15_30
## values : 8.6, 322.9 (min, max)4.Exploración de datos
Un histograma puede ayudar en este momento:
hist (boy_soc)
Un resumen también es útil:
summary(boy_soc)## soc_igh_15_30
## Min. 8.6
## 1st Qu. 26.2
## Median 36.7
## 3rd Qu. 48.8
## Max. 322.9
## NA's 14493.0Finalmente, trazando el tiempo. Primero, una conversión de capa ráster a estrellas.
(nboy_soc <- st_as_stars(boy_soc))## stars object with 2 dimensions and 1 attribute
## attribute(s), summary of first 1e+05 cells:
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## soc_igh_15_30 11 24.9 32.2 37.63789 47.2 127 1298
## dimension(s):
## from to offset delta refsys x/y
## x 1 1173 -8326000 250 +proj=igh +lon_0=0 +x_0=0... [x]
## y 1 1069 785500 -250 +proj=igh +lon_0=0 +x_0=0... [y]names(nboy_soc) <- "soc"Ahora se hace la paleta de colores:
pal <- colorRampPalette(brewer.pal(9, "BrBG"))Luego, hacesmos el mapa:
tm_shape(nboy_soc) +
tm_raster("soc", palette = pal(9), title = "SOC [%]", ) +
tm_shape(boy) + tm_borders() + tm_text("MPIO_CNMBR", size=0.7) +
tm_layout(scale = .5,
legend.position = c("left","bottom"),
legend.bg.color = "white", legend.bg.alpha = .5,
legend.frame = "black")## stars object downsampled to 1047 by 954 cells. See tm_shape manual (argument raster.downsample)
## 5. Referencias. Lizarazo, I. 2023. How to download soil data from
SoilGrids. Available at: https://rpubs.com/ials2un/gettingSoilGrids.
sessionInfo()## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 22621)
##
## 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] gdalUtilities_1.2.5 devtools_2.4.5 usethis_2.2.2
## [4] tmap_3.3-4 mapview_2.11.2 RColorBrewer_1.1-3
## [7] dplyr_1.1.3 stars_0.6-4 sf_1.0-14
## [10] abind_1.4-5 raster_3.6-26 sp_2.1-1
## [13] XML_3.99-0.15
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.0 viridisLite_0.4.2 fastmap_1.1.1 leaflet_2.2.0
## [5] promises_1.2.1 digest_0.6.33 mime_0.12 lifecycle_1.0.4
## [9] ellipsis_0.3.2 processx_3.8.2 terra_1.7-46 magrittr_2.0.3
## [13] compiler_4.3.1 rlang_1.1.1 sass_0.4.7 tools_4.3.1
## [17] utf8_1.2.3 yaml_2.3.7 knitr_1.43 prettyunits_1.1.1
## [21] htmlwidgets_1.6.2 curl_5.0.2 pkgbuild_1.4.2 classInt_0.4-10
## [25] pkgload_1.3.3 KernSmooth_2.23-21 miniUI_0.1.1.1 purrr_1.0.2
## [29] desc_1.4.2 leafsync_0.1.0 grid_4.3.1 stats4_4.3.1
## [33] fansi_1.0.4 urlchecker_1.0.1 profvis_0.3.8 xtable_1.8-4
## [37] e1071_1.7-13 leafem_0.2.3 colorspace_2.1-0 scales_1.2.1
## [41] dichromat_2.0-0.1 cli_3.6.1 rmarkdown_2.24 crayon_1.5.2
## [45] remotes_2.4.2.1 generics_0.1.3 rstudioapi_0.15.0 sessioninfo_1.2.2
## [49] tmaptools_3.1-1 DBI_1.1.3 cachem_1.0.8 proxy_0.4-27
## [53] stringr_1.5.0 parallel_4.3.1 base64enc_0.1-3 vctrs_0.6.3
## [57] jsonlite_1.8.7 callr_3.7.3 crosstalk_1.2.0 jquerylib_0.1.4
## [61] units_0.8-4 glue_1.6.2 lwgeom_0.2-13 codetools_0.2-19
## [65] ps_1.7.5 stringi_1.7.12 later_1.3.1 munsell_0.5.0
## [69] tibble_3.2.1 pillar_1.9.0 htmltools_0.5.6 satellite_1.0.4
## [73] R6_2.5.1 rprojroot_2.0.4 evaluate_0.21 shiny_1.7.5
## [77] lattice_0.21-8 highr_0.10 png_0.1-8 memoise_2.0.1
## [81] httpuv_1.6.11 bslib_0.5.1 class_7.3-22 Rcpp_1.0.11
## [85] xfun_0.40 fs_1.6.3 pkgconfig_2.0.3