¿Como descargar datos de suelo de SoilGrids
Edna Valentina Acosta Romero
2023-11-22
1. INTRODUCCION
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. CONFIGURACION
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')#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")## 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)## Warning: package 'XML' was built under R version 4.3.2library(raster)## Warning: package 'raster' was built under R version 4.3.2## Loading required package: splibrary(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 INTERES:
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 REGION DE INTERES(ROI)
Leamos un shapefile de un departamento previamente descargado del DANE. Estaremos aquí la biblioteca sf
(cundi = st_read("C:/GB/MGN_ANM_MPIOS.shp"))## Reading layer `MGN_ANM_MPIOS' from data source `C:\GB\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## 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
## First 10 features:
## DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR MPIO_CDPMP VERSION AREA
## 1 18 001 FLORENCIA 18001 2018 2547637532
## 2 18 029 ALBANIA 18029 2018 414122070
## 3 18 094 BELÉN DE LOS ANDAQUÍES 18094 2018 1191618572
## 4 18 247 EL DONCELLO 18247 2018 1106076151
## 5 18 256 EL PAUJÍL 18256 2018 1234734145
## 6 18 410 LA MONTAÑITA 18410 2018 1701061430
## 7 18 460 MILÁN 18460 2018 1220575999
## 8 18 479 MORELIA 18479 2018 462479591
## 9 18 610 SAN JOSÉ DEL FRAGUA 18610 2018 1304769006
## 10 18 860 VALPARAÍSO 18860 2018 1330210015
## LATITUD LONGITUD STCTNENCUE STP3_1_SI STP3_2_NO STP3A_RI STP3B_TCN
## 1 1.749139 -75.55824 71877 32 71845 32 0
## 2 1.227865 -75.88233 2825 24 2801 24 0
## 3 1.500923 -75.87565 4243 54 4189 54 0
## 4 1.791386 -75.19394 8809 0 8809 0 0
## 5 1.617746 -75.23404 5795 0 5795 0 0
## 6 1.302860 -75.23573 5113 15 5098 15 0
## 7 1.146693 -75.38665 3462 363 3099 363 0
## 8 1.382996 -75.67381 1907 0 1907 0 0
## 9 1.304702 -76.11001 5756 83 5673 83 0
## 10 1.061292 -75.62705 3958 0 3958 0 0
## STP4_1_SI STP4_2_NO STP9_1_USO STP9_2_USO STP9_3_USO STP9_4_USO STP9_2_1_M
## 1 0 71877 61176 2178 8436 87 39
## 2 0 2825 1826 49 948 2 3
## 3 1 4242 3223 109 900 11 4
## 4 0 8809 6598 357 1850 4 11
## 5 0 5795 4891 204 695 5 4
## 6 0 5113 4077 241 786 9 2
## 7 0 3462 2594 118 745 5 1
## 8 0 1907 1221 36 650 0 0
## 9 27 5729 4493 302 955 6 1
## 10 0 3958 2308 40 1602 8 0
## STP9_2_2_M STP9_2_3_M STP9_2_4_M STP9_2_9_M STP9_3_1_N STP9_3_2_N STP9_3_3_N
## 1 1550 566 18 5 54 2591 1061
## 2 34 12 0 0 3 21 32
## 3 99 6 0 0 8 88 6
## 4 259 87 0 0 21 334 124
## 5 161 38 1 0 5 239 104
## 6 205 32 2 0 3 103 123
## 7 104 12 1 0 0 75 41
## 8 36 0 0 0 1 27 16
## 9 273 28 0 0 4 242 82
## 10 29 11 0 0 4 35 42
## STP9_3_4_N STP9_3_5_N STP9_3_6_N STP9_3_7_N STP9_3_8_N STP9_3_9_N STP9_3_10
## 1 535 368 3172 233 7 19 371
## 2 728 53 92 5 0 0 14
## 3 2 61 626 8 0 8 93
## 4 807 89 361 42 0 27 39
## 5 4 70 211 16 0 0 45
## 6 17 84 362 30 0 0 58
## 7 39 24 514 8 0 8 35
## 8 442 31 114 8 0 2 9
## 9 0 101 307 55 0 1 159
## 10 441 56 982 19 0 3 15
## STP9_3_99 STVIVIENDA STP14_1_TI STP14_2_TI STP14_3_TI STP14_4_TI STP14_5_TI
## 1 25 63354 47817 13764 1624 21 8
## 2 0 1875 1793 40 22 17 0
## 3 0 3332 3189 113 24 2 2
## 4 6 6955 6006 775 160 1 1
## 5 1 5095 4700 145 188 4 1
## 6 6 4318 3890 224 161 6 2
## 7 1 2712 2204 85 33 386 0
## 8 0 1257 1160 81 15 0 0
## 9 4 4795 4372 262 135 9 0
## 10 5 2348 2197 52 82 1 2
## STP14_6_TI STP15_1_OC STP15_2_OC STP15_3_OC STP15_4_OC TSP16_HOG STP19_EC_1
## 1 120 49809 2681 2150 8714 51430 48638
## 2 3 1409 13 55 398 1559 1300
## 3 2 2883 2 107 340 3161 2595
## 4 12 5767 304 388 496 6129 5375
## 5 57 4568 5 323 199 5848 4195
## 6 35 3553 151 308 306 3748 2159
## 7 4 1985 4 399 324 2182 1152
## 8 1 1013 0 70 174 1233 923
## 9 17 3444 1 875 475 3873 2794
## 10 14 1813 3 199 333 1881 1697
## STP19_ES_2 STP19_EE_1 STP19_EE_2 STP19_EE_3 STP19_EE_4 STP19_EE_5 STP19_EE_6
## 1 1171 34851 10343 2169 509 13 3
## 2 109 1184 106 1 0 0 0
## 3 288 2118 366 17 2 1 0
## 4 392 3548 962 793 1 1 1
## 5 373 3330 770 84 0 1 1
## 6 1394 1964 144 9 1 0 0
## 7 833 1098 31 0 1 0 1
## 8 90 719 187 9 0 0 0
## 9 650 2486 283 7 0 1 0
## 10 116 1624 68 2 0 0 0
## STP19_EE_9 STP19_ACU1 STP19_ACU2 STP19_ALC1 STP19_ALC2 STP19_GAS1 STP19_GAS2
## 1 750 45179 4630 41138 8671 37028 12074
## 2 9 808 601 703 706 26 1371
## 3 91 2017 866 1806 1077 52 2796
## 4 69 4175 1592 4323 1444 57 5549
## 5 9 2505 2063 2359 2209 1463 3041
## 6 41 1441 2112 1329 2224 67 3454
## 7 21 359 1626 595 1390 20 1950
## 8 8 588 425 573 440 15 982
## 9 17 2576 868 2121 1323 1329 2078
## 10 3 775 1038 842 971 22 1772
## STP19_GAS9 STP19_REC1 STP19_REC2 STP19_INT1 STP19_INT2 STP19_INT9 STP27_PERS
## 1 707 45491 4318 13362 35727 720 156789
## 2 12 727 682 27 1370 12 4514
## 3 35 1905 978 73 2775 35 9075
## 4 161 4348 1419 211 5395 161 17775
## 5 64 2414 2154 125 4379 64 13014
## 6 32 1273 2280 64 3457 32 12128
## 7 15 661 1324 39 1931 15 7507
## 8 16 588 425 6 991 16 3350
## 9 37 2559 885 131 3276 37 11364
## 10 19 816 997 48 1739 26 6082
## STPERSON_L STPERSON_S STP32_1_SE STP32_2_SE STP34_1_ED STP34_2_ED STP34_3_ED
## 1 4315 152474 77620 79169 25503 30249 29951
## 2 151 4363 2323 2191 725 1016 717
## 3 346 8729 4551 4524 1592 2254 1388
## 4 203 17572 8790 8985 3047 3811 2601
## 5 192 12822 6601 6413 2346 2882 2170
## 6 604 11524 6437 5691 2229 3022 1836
## 7 527 6980 3953 3554 1550 2018 1125
## 8 0 3350 1709 1641 562 678 570
## 9 399 10965 5812 5552 2129 2822 1757
## 10 404 5678 3204 2878 1013 1592 836
## STP34_4_ED STP34_5_ED STP34_6_ED STP34_7_ED STP34_8_ED STP34_9_ED STP51_PRIM
## 1 23602 17235 14349 8969 4687 2244 48848
## 2 568 536 445 253 162 92 1940
## 3 1121 986 816 487 286 145 3541
## 4 2302 2032 1792 1135 707 348 7571
## 5 1587 1460 1188 703 430 248 6072
## 6 1563 1441 1010 578 323 126 5932
## 7 912 736 562 356 172 76 3344
## 8 428 431 304 212 114 51 1433
## 9 1494 1267 933 542 275 145 4952
## 10 696 742 608 308 187 100 2601
## STP51_SECU STP51_SUPE STP51_POST STP51_13_E STP51_99_E Shape_Leng Shape_Area
## 1 59610 21898 4592 5892 3799 2.942508 0.20692777
## 2 1712 231 41 215 46 1.112829 0.03361758
## 3 3340 490 119 720 123 2.234657 0.09674460
## 4 6287 1029 228 1095 171 3.154370 0.08986744
## 5 4066 639 108 916 99 3.529316 0.10030928
## 6 3870 407 94 724 182 3.402939 0.13817351
## 7 2669 239 38 465 86 1.863197 0.09912782
## 8 1194 215 16 181 32 1.518688 0.03755356
## 9 4155 516 87 579 94 2.040837 0.10589313
## 10 2292 241 31 368 108 2.313848 0.10800551
## geometry
## 1 MULTIPOLYGON (((-75.42074 2...
## 2 MULTIPOLYGON (((-75.89506 1...
## 3 MULTIPOLYGON (((-75.78705 1...
## 4 MULTIPOLYGON (((-75.36167 2...
## 5 MULTIPOLYGON (((-75.36638 2...
## 6 MULTIPOLYGON (((-75.40346 1...
## 7 MULTIPOLYGON (((-75.39362 1...
## 8 MULTIPOLYGON (((-75.77185 1...
## 9 MULTIPOLYGON (((-76.16722 1...
## 10 MULTIPOLYGON (((-75.73128 1...Filtremos el departamento:
cundi %>% filter(DPTO_CCDGO=="25") %>% select(MPIO_CDPMP, MPIO_CNMBR,AREA) -> Cundi1Cundi1## Simple feature collection with 116 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.89063 ymin: 3.730129 xmax: -73.05256 ymax: 5.837258
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## MPIO_CDPMP MPIO_CNMBR AREA geometry
## 1 25019 ALBÁN 50835100 MULTIPOLYGON (((-74.47852 4...
## 2 25035 ANAPOIMA 123955490 MULTIPOLYGON (((-74.55201 4...
## 3 25040 ANOLAIMA 120952347 MULTIPOLYGON (((-74.43043 4...
## 4 25053 ARBELÁEZ 142459967 MULTIPOLYGON (((-74.42848 4...
## 5 25086 BELTRÁN 176845987 MULTIPOLYGON (((-74.76139 4...
## 6 25095 BITUIMA 61046194 MULTIPOLYGON (((-74.51062 4...
## 7 25099 BOJACÁ 102319881 MULTIPOLYGON (((-74.32819 4...
## 8 25120 CABRERA 421839593 MULTIPOLYGON (((-74.46463 4...
## 9 25123 CACHIPAY 53452616 MULTIPOLYGON (((-74.40513 4...
## 10 25126 CAJICÁ 51275631 MULTIPOLYGON (((-73.99959 4...Tenga en cuenta el CRS de stder. 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'
cundi_igh <- st_transform(Cundi1, igh)Consigamos el cuadro delimitador de nuestra área de interés:
(bbox <- st_bbox(cundi_igh))## xmin ymin xmax ymax
## -8344613.5 415236.1 -8142417.0 649800.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.
## 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
## -8344613.5 649800.5 -8142417.0 415236.1Ahora 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.
(cundi1_soc <- raster(file)/10)## class : RasterLayer
## dimensions : 938, 809, 758842 (nrow, ncol, ncell)
## resolution : 250, 250 (x, y)
## extent : -8344750, -8142500, 415500, 650000 (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.9, 322.9 (min, max)5. EXPLORACION DE LOS DATOS
Un histograma puede ayudar en este momento:
hist(cundi1_soc)
El SOC estara en los intervalos de 0 a 120 en el mapa, predominando el
de 0 a 50%, tienen muy poco carbono organico
Un resumen también es útil:
summary(cundi1_soc)## soc_igh_15_30
## Min. 8.9
## 1st Qu. 25.9
## Median 38.3
## 3rd Qu. 51.4
## Max. 322.9
## NA's 13496.0Finalmente, trazando el tiempo. Primero, una conversión de capa ráster a estrellas.
(ncundi_soc <- st_as_stars(cundi1_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.9 25.9 38.3 39.68411 51.4 322.9 13496
## dimension(s):
## from to offset delta refsys x/y
## x 1 809 -8344750 250 +proj=igh +lon_0=0 +x_0=0... [x]
## y 1 938 650000 -250 +proj=igh +lon_0=0 +x_0=0... [y]names(ncundi_soc) <- "soc"Luego, la paleta de colores:
pal <- colorRampPalette(brewer.pal(5, "YlOrRd"))Entonces, el mapa:
tm_shape(ncundi_soc) +
tm_raster("soc", palette = pal(5), title = "SOC [%]", ) +
tm_shape(Cundi1) + tm_borders() + tm_text("MPIO_CNMBR", size=0.7, col = "black") +
tm_layout(scale = .7,
legend.position = c("right","bottom"),
legend.bg.color = "white", legend.bg.alpha = .10,
legend.frame = "black")
Viendo el mapa los valores de 0 a 50% son paratebueno, medina, puerto
salgar, tocaiman, guataqui, agua de dios, entre otros, luego ene l rango
de 50 a 100% estaria gutierrez, venecia, junin, une entre otros y en el
rangp de 100 a 150% esta tequendama, teniendo un numero mayor de
porcentaje de carbono.
## 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 10 x64 (build 19045)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=Spanish_Mexico.utf8 LC_CTYPE=Spanish_Mexico.utf8
## [3] LC_MONETARY=Spanish_Mexico.utf8 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Mexico.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.4 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.3
## [9] ellipsis_0.3.2 processx_3.8.2 terra_1.7-55 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.4 yaml_2.3.7 knitr_1.44 prettyunits_1.2.0
## [21] htmlwidgets_1.6.2 pkgbuild_1.4.2 classInt_0.4-10 curl_5.1.0
## [25] pkgload_1.3.3 KernSmooth_2.23-21 miniUI_0.1.1.1 withr_2.5.1
## [29] purrr_1.0.2 leafsync_0.1.0 grid_4.3.1 stats4_4.3.1
## [33] fansi_1.0.5 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.25 crayon_1.5.2
## [45] generics_0.1.3 remotes_2.4.2.1 rstudioapi_0.15.0 tmaptools_3.1-1
## [49] sessioninfo_1.2.2 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.4
## [57] jsonlite_1.8.7 callr_3.7.3 crosstalk_1.2.0 jquerylib_0.1.4
## [61] units_0.8-4 lwgeom_0.2-13 glue_1.6.2 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.1 satellite_1.0.4
## [73] R6_2.5.1 evaluate_0.22 shiny_1.8.0 lattice_0.22-5
## [77] png_0.1-8 memoise_2.0.1 httpuv_1.6.12 bslib_0.5.1
## [81] class_7.3-22 Rcpp_1.0.11 xfun_0.40 fs_1.6.3
## [85] pkgconfig_2.0.3