Descargar datos de suelos en SoilGrids para el departamento de Antioquia

1. Introducción

SoilGrids es un sistema de cartografía digital global del suelo que utiliza información global sobre perfiles de suelo y datos de covariables para modelizar la distribución espacial de las propiedades del suelo en todo el planeta. SoilGrids es una colección de mapas de propiedades del suelo de todo el mundo elaborados 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 Referencia e Información sobre el Suelo - (ISRIC)[https://www.isric.org/explore/soilgrids]:

Web Mapping Service (WMS): acceso para visualización y resumen de datos.

Web Coverage Service (WCS): la mejor manera de obtener un subconjunto de un mapa y utilizar SoilGrids como entrada para otros procesos de modelización.

Web Distributed Authoring and Versioning (WebDAV): descarga del mapa o mapas globales completos en formato VRT.

En este cuaderno utilizaremos el servicio WebDAV. En caso de duda, consulte las instrucciones de SoilGrids WebDAV.

Traducción realizada con la versión gratuita del traductor www.DeepL.com/Translator

2. Configuración

Limpiamos la mermoria:

rm(list = ls())

Instalamos las librerías requridas para este libro:

library(XML)
library(raster)
library(stars)
library(sf)
library(dplyr)
library(RColorBrewer)
library(mapview)
library(tmap)

library(devtools)

library(gdalUtilities)

Establezca las variables de interés

Primero, debemos definir la url correspondiente al sitio webDAD del ISRIC:

url = "https://files.isric.org/soilgrids/latest/data/"

Ahora, vamos a crear algunos objetos que indiquen aquello 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"

Lo siguiente será definir 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"

La espicificación del VTR está completa

vtr_layer = paste(layer, '.vrt', sep = "")

3. Descargar un Tiff para la región de interés (ROI)

Vamos a leer un shapefile de un departamento anteriormente descargado del DANE. Aquí usaremos la biblioteca sf.

Antiq = st_read("Mun_Antiq.shp")

Reading layer `Mun_Antiq' from data source 
  `C:\Users\ASUS\Documents\GB2\SoilGrids\Mun_Antiq.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 125 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -77.12783 ymin: 5.418558 xmax: -73.88128 ymax: 8.873974
Geodetic CRS:  WGS 84

Antiq

Simple feature collection with 125 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -77.12783 ymin: 5.418558 xmax: -73.88128 ymax: 8.873974
Geodetic CRS:  WGS 84
First 10 features:
   MPIO_CCDGO MPIO_CDPMP DPTO_CNMBR  MPIO_CNMBR
1         001      05001  ANTIOQUIA    MEDELLÍN
2         002      05002  ANTIOQUIA   ABEJORRAL
3         004      05004  ANTIOQUIA    ABRIAQUÍ
4         021      05021  ANTIOQUIA  ALEJANDRÍA
5         030      05030  ANTIOQUIA       AMAGÁ
6         031      05031  ANTIOQUIA      AMALFI
7         034      05034  ANTIOQUIA       ANDES
8         036      05036  ANTIOQUIA ANGELÓPOLIS
9         038      05038  ANTIOQUIA   ANGOSTURA
10        040      05040  ANTIOQUIA       ANORÍ
                             MPIO_CRSLC MPIO_NAREA MPIO_TIPO Shape_Leng
1                                  1965  374.83400 MUNICIPIO  1.0353800
2                                  1814  507.14109 MUNICIPIO  1.1585038
3                                  1912  296.89405 MUNICIPIO  0.8121832
4     Decreto departamental 304 de 1907  128.93215 MUNICIPIO  0.7051995
5                                  1912   84.13268 MUNICIPIO  0.4452407
6                                  1843 1209.14802 MUNICIPIO  2.0633596
7                                  1853  402.46597 MUNICIPIO  1.1469836
8  Ordenanza 16 del 12 de julio de 1896   81.87630 MUNICIPIO  0.4274956
9                                  1821  338.50229 MUNICIPIO  0.9164653
10                                 1821 1413.77572 MUNICIPIO  1.8476441
    Shape_Area                       geometry
1  0.030607649 MULTIPOLYGON (((-75.66974 6...
2  0.041383896 MULTIPOLYGON (((-75.46938 5...
3  0.024248255 MULTIPOLYGON (((-76.08351 6...
4  0.010534527 MULTIPOLYGON (((-75.0332 6....
5  0.006866523 MULTIPOLYGON (((-75.67587 6...
6  0.098921310 MULTIPOLYGON (((-74.92268 6...
7  0.032816580 MULTIPOLYGON (((-75.86822 5...
8  0.006683387 MULTIPOLYGON (((-75.69149 6...
9  0.027679560 MULTIPOLYGON (((-75.27173 6...
10 0.115706037 MULTIPOLYGON (((-74.90935 7...

Observe el CRS de Antiq. Como las capas ISRIC utilizan la proyección Homolosine, necesitamos reproyectar nuestra capa utilizando la biblioteca sf:

igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
Antiq_igh <- st_transform(Antiq, igh)

Ahora vamos a mostrar el bounding box de nuestra área de interés:

(bbox <- st_bbox(Antiq_igh))

      xmin       ymin       xmax       ymax 
-8609547.4   603191.2 -8246194.5   987846.2 

Utilizaremos los datos bbox para definir los límites de la caja límite tal y como los utiliza la librería GDAL. Por cierto, esta es una de las partes más complicadas de usar GDAL.

ulx = bbox$xmin
uly = bbox$ymax
lrx= bbox$xmax
lry = bbox$ymin
(bb <- c(ulx, uly, lrx, lry))

      xmin       ymax       xmax       ymin 
-8609547.4   987846.2 -8246194.5   603191.2 

Ahora, podemos utilizar la función gdal_translate para descargar la capa vrt. En primer lugar, vamos a definir 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"

Lo siguiente tardará unos minutos

gdal_translate(paste0(sg_url,datos), file ,
               tr=c(250,250),
               projwin=bb,
               projwin_srs =igh)

Recordemos qué son las unidades de capa SOC:

Dividiendo los valores de las predicciones por los valores de la columna Factor de conversión, el usuario puede obtener las unidades más familiares de la columna Unidades convencionales.

(Antiq_soc <- raster(file)/10)

class      : RasterLayer 
dimensions : 1539, 1453, 2236167  (nrow, ncol, ncell)
resolution : 250, 250  (x, y)
extent     : -8609750, -8246500, 603250, 988000  (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     : 7.8, 350.7  (min, max)

4. Exploración de Datos

Un histograma puede ayudar en este momento:

hist(Antiq_soc)

También es útil un resumen:

summary(Antiq_soc)

        soc_igh_15_30
Min.              7.8
1st Qu.          27.9
Median           42.0
3rd Qu.          70.4
Max.            350.7
NA's         111105.0

Y para finalizar vamos a trazar. Primero una conversión de capa raster a estrellas

(nAntiq_soc <- st_as_stars(Antiq_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   8.3    21.3   29.8 32.59002    40.8 112.7 29172
dimension(s):

names(nAntiq_soc) <-  "soc"

Ahora, el color de la paleta:

pal <- colorRampPalette(brewer.pal(9, "BrBG"))

Y por último el mapa:

tm_shape(nAntiq_soc) +
  tm_raster("soc", palette = pal(10), title = "SOC [%]", )  +  
  tm_shape(Antiq) + 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 = "black")

stars object downsampled to 971 by 1029 cells. See tm_shape manual (argument raster.downsample)

Tenga en cuenta que puede cambiar la paleta de colores y la disposición del mapa para comprender mejor la información de la capa de suelo.

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] 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.2         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] pkgload_1.3.3      viridisLite_0.4.2  shiny_1.7.5.1     
 [4] stats4_4.2.3       remotes_2.4.2.1    sessioninfo_1.2.2 
 [7] pillar_1.9.0       lattice_0.22-5     glue_1.6.2        
[10] digest_0.6.33      promises_1.2.1     colorspace_2.1-0  
[13] htmltools_0.5.7    httpuv_1.6.12      pkgconfig_2.0.3   
[16] purrr_1.0.2        xtable_1.8-4       scales_1.2.1      
[19] processx_3.8.2     terra_1.7-55       satellite_1.0.4   
[22] later_1.3.1        tibble_3.2.1       proxy_0.4-27      
[25] generics_0.1.3     ellipsis_0.3.2     cachem_1.0.8      
[28] leafsync_0.1.0     cli_3.6.1          magrittr_2.0.3    
[31] crayon_1.5.2       mime_0.12          memoise_2.0.1     
[34] 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    
[40] profvis_0.3.8      tools_4.2.3        prettyunits_1.2.0 
[43] lifecycle_1.0.3    stringr_1.5.1      munsell_0.5.0     
[46] callr_3.7.3        compiler_4.2.3     e1071_1.7-13      
[49] rlang_1.1.1        classInt_0.4-10    units_0.8-4       
[52] grid_4.2.3         tmaptools_3.1-1    dichromat_2.0-0.1 
[55] rstudioapi_0.15.0  htmlwidgets_1.6.2  crosstalk_1.2.0   
[58] miniUI_0.1.1.1     leafem_0.2.3       base64enc_0.1-3   
[61] codetools_0.2-19   curl_5.1.0         DBI_1.1.3         
[64] R6_2.5.1           knitr_1.45         fastmap_1.1.1     
[67] utf8_1.2.3         rprojroot_2.0.4    desc_1.4.2        
[70] KernSmooth_2.23-20 stringi_1.8.1      parallel_4.2.3    
[73] Rcpp_1.0.11        vctrs_0.6.3        png_0.1-8         
[76] urlchecker_1.0.1   leaflet_2.2.1      tidyselect_1.2.0  
[79] xfun_0.41