¿Cómo descargar datos de suelo de SoilGrids?
Valentina Díaz Díaz
16.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.
2. Configuración
-Limpiar el espacio de trabajo:
rm(list=ls())-Instale las bibliotecas 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')-Cargue las bibliotecas:
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.2install_github("JoshOBrien/gdalUtilities")## Downloading GitHub repo JoshOBrien/gdalUtilities@HEAD##
## ── R CMD build ─────────────────────────────────────────────────────────────────
##
checking for file 'C:\Users\gagug\AppData\Local\Temp\RtmpaivzY5\remotes4c7c712a584f\JoshOBrien-gdalUtilities-3d87b46/DESCRIPTION' ...
checking for file 'C:\Users\gagug\AppData\Local\Temp\RtmpaivzY5\remotes4c7c712a584f\JoshOBrien-gdalUtilities-3d87b46/DESCRIPTION' ...
✔ checking for file 'C:\Users\gagug\AppData\Local\Temp\RtmpaivzY5\remotes4c7c712a584f\JoshOBrien-gdalUtilities-3d87b46/DESCRIPTION' (429ms)
##
─ preparing 'gdalUtilities':
## checking DESCRIPTION meta-information ...
✔ checking DESCRIPTION meta-information
##
─ installing the package to process help pages
##
─ saving partial Rd database (9s)
##
─ 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/RtmpaivzY5/file4c7c8e51682/gdalUtilities_1.2.5.tar.gz'
## had non-zero exit statuslibrary(gdalUtilities)## Warning: package 'gdalUtilities' was built under R version 4.3.2library(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 TRUE##
## Attaching package: 'sf'## The following object is masked from 'package:gdalUtilities':
##
## gdal_rasterizelibrary(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)3.Establecer variables de interés:
-Primero, definimos la url correspondiente al sitio webDAD de ISRIC:
url = "https://files.isric.org/soilgrids/latest/data/" # Path to the webDAV data.-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="_")) # layer of interest ## [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"4. Descargue un TIFF para una región de interés (ROI)
-Leamos un shapefile de un departamento previamente descargado del DANE. Encontramos aquí la biblioteca sf:
(tol= st_read("C:/Users/gagug/OneDrive/Escritorio/gb_vale/Soil_datos/SHP_MGN2018_INTGRD_MPIO (1)/MGN_ANM_MPIOS.shp"))## Reading layer `MGN_ANM_MPIOS' from data source
## `C:\Users\gagug\OneDrive\Escritorio\gb_vale\Soil_datos\SHP_MGN2018_INTGRD_MPIO (1)\MGN_ANM_MPIOS.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1122 features and 90 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -81.73562 ymin: -4.229406 xmax: -66.84722 ymax: 13.39473
## Geodetic CRS: MAGNA-SIRGAS-Filtremos el departamento
tol %>% filter(DPTO_CCDGO=="73") %>% select(MPIO_CDPMP , MPIO_CNMBR,AREA) -> tolimtolim-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'
tol_igh <- st_transform(tolim,igh)-Consigamos el cuadro delimitador de nuestra área de interés:
(bbox <- st_bbox(tol_igh))## xmin ymin xmax ymax
## -8475961.7 319607.2 -8297873.1 592146.4-Ahora, usemos los datos de bbox para definir los límites del cuadro de límites tal como los usa la biblioteca GDA:
## ul means upper left
## lr means lower right
ulx = bbox$xmin
uly = bbox$ymax
lrx= bbox$xmax
lry = bbox$ymin
(bb <- c(ulx, uly, lrx, lry))## xmin ymax xmax ymin
## -8475961.7 592146.4 -8297873.1 319607.2-Ahora 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)-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.
(tol_soc <- raster(file)/10)## class : RasterLayer
## dimensions : 1090, 712, 776080 (nrow, ncol, ncell)
## resolution : 250, 250 (x, y)
## extent : -8476000, -8298000, 319750, 592250 (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 : 10.4, 122.3 (min, max)5. Exploración de datos
-Un histograma puede ayudar en este momento:
hist(tol_soc)
-Un resumen también es útil:
summary(tol_soc)## soc_igh_15_30
## Min. 10.4
## 1st Qu. 25.3
## Median 38.4
## 3rd Qu. 57.6
## Max. 122.3
## NA's 5735.0-Finalmente, trazando el tiempo. Primero, una conversión de capa ráster a estrellas.
(ntol_soc <- st_as_stars(tol_soc))## stars object with 2 dimensions and 1 attribute
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## soc_igh_15_30 10.4 25.3 38.4 41.97771 57.6 122.3 5735
## dimension(s):
## from to offset delta refsys x/y
## x 1 712 -8476000 250 +proj=igh +lon_0=0 +x_0=0... [x]
## y 1 1090 592250 -250 +proj=igh +lon_0=0 +x_0=0... [y]names(ntol_soc) <- "soc"-Luego, la paleta de colores:
pal <- colorRampPalette(brewer.pal(11, "BrBG"))-Entonces, el mapa:
tm_shape(ntol_soc) +
tm_raster("soc", palette = pal(9), title = "SOC [%]", ) +
tm_shape(tolim) + 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 = .10,
legend.frame = "black")
6. 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] tmap_3.3-4 mapview_2.11.2 RColorBrewer_1.1-3
## [4] dplyr_1.1.3 stars_0.6-4 sf_1.0-14
## [7] abind_1.4-5 raster_3.6-26 sp_2.1-1
## [10] XML_3.99-0.15 gdalUtilities_1.2.5 devtools_2.4.5
## [13] usethis_2.2.2
##
## 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 pkgbuild_1.4.2 classInt_0.4-10 curl_5.0.2
## [25] pkgload_1.3.3 KernSmooth_2.23-21 miniUI_0.1.1.1 withr_2.5.2
## [29] purrr_1.0.2 leafsync_0.1.0 desc_1.4.2 grid_4.3.1
## [33] stats4_4.3.1 fansi_1.0.4 urlchecker_1.0.1 profvis_0.3.8
## [37] xtable_1.8-4 e1071_1.7-13 leafem_0.2.3 colorspace_2.1-0
## [41] scales_1.2.1 dichromat_2.0-0.1 cli_3.6.1 rmarkdown_2.24
## [45] crayon_1.5.2 generics_0.1.3 remotes_2.4.2.1 rstudioapi_0.15.0
## [49] tmaptools_3.1-1 sessioninfo_1.2.2 DBI_1.1.3 cachem_1.0.8
## [53] proxy_0.4-27 stringr_1.5.0 parallel_4.3.1 base64enc_0.1-3
## [57] vctrs_0.6.3 jsonlite_1.8.7 callr_3.7.3 crosstalk_1.2.0
## [61] jquerylib_0.1.4 units_0.8-4 lwgeom_0.2-13 glue_1.6.2
## [65] codetools_0.2-19 ps_1.7.5 stringi_1.7.12 later_1.3.1
## [69] munsell_0.5.0 tibble_3.2.1 pillar_1.9.0 htmltools_0.5.6
## [73] satellite_1.0.4 R6_2.5.1 rprojroot_2.0.4 evaluate_0.21
## [77] shiny_1.7.5 lattice_0.21-8 highr_0.10 png_0.1-8
## [81] memoise_2.0.1 httpuv_1.6.11 bslib_0.5.1 class_7.3-22
## [85] Rcpp_1.0.11 xfun_0.40 fs_1.6.3 pkgconfig_2.0.3