MAPEO DIGITAL DEL DEPARTAMENTO DEL TOLIMA USANDO SOILGRIDS
Santiago Andres Rodriguez Rueda
19/11/2023
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
Iniciaremos por limpiar el entorno de trabajo:
rm (list = ls())Cargamos las librerias:
library(raster)## Loading required package: splibrary(devtools)## Loading required package: usethislibrary(devtools)
library(sf)## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUElibrary(raster)library(gdalUtilities)##
## Attaching package: 'gdalUtilities'## The following object is masked from 'package:sf':
##
## gdal_rasterizelibrary(XML)library(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(devtools)
library(raster)
library(stars)## Loading required package: abindlibrary(RColorBrewer)
library(mapview)library(tmap)## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')Tenga en cuenta que, para utilizar gdalUtils, necesitamos instalar GDAL en su sistema. Más información aquí.
Establecer una Variable
Iniciamos por definir el url del sitio ISRIC webDAD:
url = "https://files.isric.org/soilgrids/latest/data/" Luego creamos algun tipo de objetos indicando lo que necesitamos de ISRIC:
voi = "soc"
depth = "15-30cm"
quantile = "mean"
(variable = paste(url, voi, sep=""))## [1] "https://files.isric.org/soilgrids/latest/data/soc"Ahora definimos la propiedad del suelo que nos interesa observar y pasamos a 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. Descarguar un Tiff para nuestra región de interes (ROI)
Leamos un shapefile de un departamento previamente descargado del DANE. Estaremos aquí la biblioteca sf:
list.files()## [1] "rsconnect" "soc_igh_15_30.tif" "SOILGRIDS.html"
## [4] "SOILGRIDS.nb.html" "SOILGRIDS.Rmd" "Tol_MPIO.shp.CPG"
## [7] "Tol_MPIO.shp.dbf" "Tol_MPIO.shp.prj" "Tol_MPIO.shp.qix"
## [10] "Tol_MPIO.shp.shp" "Tol_MPIO.shp.shx" "Tol_MPIO.shp.xml"(Tol = st_read("./Tol_MPIO.shp.shp"))## Reading layer `Tol_MPIO.shp' from data source
## `C:\Users\SEBAS\Downloads\GB2\SOILGRIDS_2\Tol_MPIO.shp.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 46 features and 11 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -76.10574 ymin: 2.871081 xmax: -74.47482 ymax: 5.313035
## Geodetic CRS: WGS 84Tenga en cuenta el CRS del Tolima. Como las capas ISRIC usan la proyección Homolosine, necesitamos reproyectar nuestra capa usando la biblioteca sf, debemos tener encuenta de estar manejando el mismo sistema de referencia de coordenadas (CRS):
shapefile_path <- "./Tol_MPIO.shp.shp"
original_data <- st_read(shapefile_path)## Reading layer `Tol_MPIO.shp' from data source
## `C:\Users\SEBAS\Downloads\GB2\SOILGRIDS_2\Tol_MPIO.shp.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 46 features and 11 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -76.10574 ymin: 2.871081 xmax: -74.47482 ymax: 5.313035
## Geodetic CRS: WGS 84igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
Tol_igh <- st_transform(Tol, igh)Consigamos el cuadro delimitador para el departamento del Tolima:
(bbox <- st_bbox(Tol_igh))## xmin ymin xmax ymax
## -8475961.7 319607.2 -8297873.1 591444.3Ahora, 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
## -8475961.7 591444.3 -8297873.1 319607.2Ahora 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"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.
(Tol_soc <- raster(file)/10)## class : RasterLayer
## dimensions : 1087, 712, 773944 (nrow, ncol, ncell)
## resolution : 250, 250 (x, y)
## extent : -8476000, -8298000, 319750, 591500 (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)4. 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 5716.0Finalmente, 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.9567 57.6 122.3 5716
## 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 1087 591500 -250 +proj=igh +lon_0=0 +x_0=0... [y](names(nTol_soc) <- ("soc"))## [1] "soc"Luego, la paleta de colores:
pal <- viridisLite::viridis(9, option = "B")
tm_shape(nTol_soc) +
tm_raster("soc", palette = pal, title = "Contenido de SOC (%)") +
tm_shape(Tol) +
tm_borders(lwd = 1, col = "white") +
tm_text("MPIO_CNMBR", size = 0.6, col = "white") +
tm_layout(
title = "Contenido de SOC en el Departamento del Tolima",
title.position = c("center", "top"),
title.size = 1.5,
title.bg.color = "white",
scale = 0.8,
legend.position = c("right", "bottom"),
legend.bg.color = "white",
legend.bg.alpha = 0.2,
legend.frame = "gray50"
)
Al detallar el resltado, es posible identificar patrones geoespaciales que revelen relaciones entre el contenido de carbono orgánico del suelo y características topográficas clave. Por ejemplo, se observa que áreas con pendientes más pronunciadas tienden a tener niveles más alto de carbono orgánico y, tambien se observa que en las áreas urbanas tienden a mostrar niveles diferentes de carbono orgánico en comparación con las zonas rurales.
5. Referencias
Este cuaderno se hace en base al trabajo de Lizarazo, I. 2023. (https://rpubs.com/ials2un/gettingSoilGrids) Este cuaderno lo puede encontrar en (https://rpubs.com/alejandrotrian/SoilGrids)
El cómo descargar datos de suelo desde SoilGrids, Disponible en: https://rpubs.com/ials2un/gettingSoilGrids
sessionInfo()## R version 4.2.3 (2023-03-15 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## 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] stars_0.6-4 abind_1.4-5 dplyr_1.1.4
## [7] XML_3.99-0.15 gdalUtilities_1.2.5 sf_1.0-14
## [10] devtools_2.4.5 usethis_2.2.2 raster_3.6-26
## [13] sp_2.1-1
##
## 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 shiny_1.8.0 highr_0.10 stats4_4.2.3
## [9] yaml_2.3.7 remotes_2.4.2.1 sessioninfo_1.2.2 pillar_1.9.0
## [13] lattice_0.20-45 glue_1.6.2 digest_0.6.33 promises_1.2.1
## [17] colorspace_2.1-0 htmltools_0.5.7 httpuv_1.6.12 pkgconfig_2.0.3
## [21] purrr_1.0.2 xtable_1.8-4 scales_1.2.1 processx_3.8.2
## [25] satellite_1.0.4 terra_1.7-55 later_1.3.1 tibble_3.2.1
## [29] proxy_0.4-27 generics_0.1.3 ellipsis_0.3.2 cachem_1.0.8
## [33] leafsync_0.1.0 cli_3.6.1 magrittr_2.0.3 crayon_1.5.2
## [37] mime_0.12 memoise_2.0.1 evaluate_0.23 ps_1.7.5
## [41] fs_1.6.3 fansi_1.0.5 lwgeom_0.2-13 class_7.3-21
## [45] pkgbuild_1.4.2 profvis_0.3.8 tools_4.2.3 prettyunits_1.2.0
## [49] lifecycle_1.0.4 stringr_1.5.1 munsell_0.5.0 callr_3.7.3
## [53] compiler_4.2.3 jquerylib_0.1.4 e1071_1.7-13 rlang_1.1.2
## [57] tmaptools_3.1-1 classInt_0.4-10 units_0.8-4 grid_4.2.3
## [61] dichromat_2.0-0.1 rstudioapi_0.15.0 htmlwidgets_1.6.2 crosstalk_1.2.0
## [65] miniUI_0.1.1.1 leafem_0.2.3 base64enc_0.1-3 rmarkdown_2.25
## [69] codetools_0.2-19 DBI_1.1.3 R6_2.5.1 knitr_1.45
## [73] fastmap_1.1.1 utf8_1.2.4 KernSmooth_2.23-20 stringi_1.8.1
## [77] parallel_4.2.3 Rcpp_1.0.11 png_0.1-8 vctrs_0.6.4
## [81] leaflet_2.2.1 tidyselect_1.2.0 xfun_0.41 urlchecker_1.0.1