descarga de los datos de suelo de SoilGrids para el departamento de Santander
Leidy Stefanny Contreras Gomez
16.11.2022
#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, revise las instrucciones de SoilGrids WebDAV.
#2. Configuraciones
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.packages(leaflet)
#install.packages(mapview)
##install.packages("devtools")
##devtools:::install_github('gearslaboratory/gdalUtils')Cargue las bibliotecas:
library(devtools)## Loading required package: usethisinstall_github("JoshOBrien/gdalUtilities")## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install Rtools 4.2 from https://cran.r-project.org/bin/windows/Rtools/.## Skipping install of 'gdalUtilities' from a github remote, the SHA1 (3d87b464) has not changed since last install.
## Use `force = TRUE` to force installationlibrary(gdalUtilities)library(XML)library(raster)## Loading required package: splibrary(stars)## Loading required package: abind## Loading required package: sf## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE##
## Attaching package: 'sf'## The following object is masked from 'package:gdalUtilities':
##
## gdal_rasterizelibrary(sf)
library(dplyr)##
## 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)library(tmap)## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')library(devtools)devtools::install_github("r-lib/devtools")## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install Rtools 4.2 from https://cran.r-project.org/bin/windows/Rtools/.## Skipping install of 'devtools' from a github remote, the SHA1 (bfb42091) has not changed since last install.
## Use `force = TRUE` to force installationTenga en cuenta que, para utilizar gdalUtils, necesitamos instalar GDAL en su sistema. Mas informacion en https://gdal.org/download.html
Establecer variables de interés:
Primero, definimos la url correspondiente al sitio webDAD de ISRIC:
url = "https://files.isric.org/soilgrids/latest/data/"Luego, creamos algunos objetos indicando lo que necesitamos de ISRIC:
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="_")) ## [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:
install.packages("sf") ## Warning: package 'sf' is in use and will not be installedlibrary(sf)(stder = st_read('datos/santander-mun.shp'))## Reading layer `santander-mun' from data source
## `E:\CuadernoAguacate\datos\santander-mun.shp' using driver `ESRI Shapefile'
## Simple feature collection with 87 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.52895 ymin: 5.719626 xmax: -72.5161 ymax: 8.119248
## Geodetic CRS: MAGNA-SIRGASTenga en cuenta el CRS de santander. 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'
stder_igh <- st_transform(stder, igh)Consigamos el cuadro delimitador de nuestra área de interés:
(bbox <- st_bbox(stder_igh))## xmin ymin xmax ymax
## -8313449.6 636705.8 -8093778.6 903830.5Ahora, 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.
ulx = bbox$xmin
uly = bbox$ymax
lrx= bbox$xmax
lry = bbox$ymin
(bb <- c(ulx, uly, lrx, lry))## xmin ymax xmax ymin
## -8313449.6 903830.5 -8093778.6 636705.8Ahora 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"devtools::install_github("r-lib/devtools")## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install Rtools 4.2 from https://cran.r-project.org/bin/windows/Rtools/.## Skipping install of 'devtools' from a github remote, the SHA1 (bfb42091) has not changed since last install.
## Use `force = TRUE` to force installationTenga 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:
[soil.png]
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.
(stder_soc <- raster(file)/10)## class : RasterLayer
## dimensions : 1068, 879, 938772 (nrow, ncol, ncell)
## resolution : 250, 250 (x, y)
## extent : -8313500, -8093750, 637000, 904000 (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.7, 130.3 (min, max)4. Exploración de datos.
Un histograma puede ayudar en este momento:
hist(stder_soc)
Un resumen también es útil:
summary(stder_soc)## soc_igh_15_30
## Min. 8.7
## 1st Qu. 25.1
## Median 33.1
## 3rd Qu. 44.2
## Max. 130.3
## NA's 14971.0Finalmente, trazando el tiempo. Primero, una conversión de capa ráster a estrellas.
(nstder_soc <- st_as_stars(stder_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 8.7 25.1 33.1 35.62564 44.2 130.3 14971
## dimension(s):
## from to offset delta refsys x/y
## x 1 879 -8313500 250 +proj=igh +lon_0=0 +x_0=0... [x]
## y 1 1068 904000 -250 +proj=igh +lon_0=0 +x_0=0... [y]Luego, la paleta de colores:
names(nstder_soc) <- "soc"pal <- colorRampPalette(brewer.pal(9, "BrBG"))Entonces, el mapa:
library(tmap)tm_shape(nstder_soc) +
tm_raster("soc", palette = pal(10), title = "SOC [%]", ) +
tm_shape(stder) + tm_borders() + tm_text("MPIO_CNMBR", size=0.5) +
tm_layout(scale = .8,
legend.position = c("left","bottom"),
legend.bg.color = "white", legend.bg.alpha = .2,
legend.frame = "gray50")
#5. Referencia
sessionInfo()## R version 4.2.3 (2023-03-15 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 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
##
## 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.2 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.9000
## [13] usethis_2.2.2
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.7 pkgload_1.3.3 viridisLite_0.4.2 jsonlite_1.8.7
## [5] bslib_0.5.1 highr_0.10 stats4_4.2.3 yaml_2.3.7
## [9] remotes_2.4.2.1 sessioninfo_1.2.2 pillar_1.9.0 lattice_0.22-5
## [13] glue_1.6.2 digest_0.6.33 colorspace_2.1-0 htmltools_0.5.7
## [17] pkgconfig_2.0.3 purrr_1.0.2 scales_1.2.1 processx_3.8.2
## [21] terra_1.7-55 satellite_1.0.4 tibble_3.2.1 proxy_0.4-27
## [25] generics_0.1.3 ellipsis_0.3.2 cachem_1.0.8 leafsync_0.1.0
## [29] cli_3.6.1 magrittr_2.0.3 crayon_1.5.2 memoise_2.0.1
## [33] evaluate_0.23 ps_1.7.5 fs_1.6.3 fansi_1.0.4
## [37] lwgeom_0.2-13 class_7.3-21 pkgbuild_1.4.2 tools_4.2.3
## [41] prettyunits_1.2.0 lifecycle_1.0.4 munsell_0.5.0 callr_3.7.3
## [45] compiler_4.2.3 jquerylib_0.1.4 e1071_1.7-13 rlang_1.1.1
## [49] classInt_0.4-10 units_0.8-4 grid_4.2.3 tmaptools_3.1-1
## [53] dichromat_2.0-0.1 rstudioapi_0.15.0 htmlwidgets_1.6.2 crosstalk_1.2.0
## [57] leafem_0.2.3 base64enc_0.1-3 rmarkdown_2.25 codetools_0.2-19
## [61] DBI_1.1.3 curl_5.1.0 R6_2.5.1 knitr_1.45
## [65] fastmap_1.1.1 utf8_1.2.3 KernSmooth_2.23-20 parallel_4.2.3
## [69] Rcpp_1.0.11 vctrs_0.6.3 png_0.1-8 leaflet_2.2.1
## [73] tidyselect_1.2.0 xfun_0.41