How to download soil data from SoilGrids

1. Introduction
2. Setup
- Set variables of interest:
3. Download a Tiff for a region of interest (ROI)
4. Data exploration
5. Reference

1. Introduction

SoilGrids is a system for global digital soil mapping that makes use of global soil profile information and covariate data to model the spatial distribution of soil properties across the globe. SoilGrids is a collections of soil property maps for the world produced using machine learning at 250 m resolution.

SoilGrids can be accessed through the following services provided by the International Soil Reference and Information Centre - (ISRIC)[https://www.isric.org/explore/soilgrids]:

Web Mapping Service (WMS): access for visualisation and data overview.
Web Coverage Service (WCS): best way to obtain a subset of a map and use SoilGrids as input to other modelling pipelines.
Web Distributed Authoring and Versioning (WebDAV): download the complete global map(s) in VRT format.

We will use the WebDAV service in this notebook. In case of any doubt, please review the SoilGrids WebDAV instructions.

2. Setup

Clean the workspace:

rm(list=ls())

Install the libraries from the console:

#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')

Load the libraries:

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

Note that, in order to use gdalUtils, we need to install GDAL in your system. More information here.

Set variables of interest:

First, we define the url corresponding to the ISRIC webDAD site:

url = "https://files.isric.org/soilgrids/latest/data/" # Path to the webDAV data.

Then, we create some objects indicating what we need from 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"

Now, we define the soil property we want to download>

#
(layer = paste(variable,depth,quantile, sep="_")) # layer of interest

## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean"

Then, the VRT specification is complete:

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

## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean.vrt"

3. Download a Tiff for a region of interest (ROI)

Let’s read a shapefile of a department previously downloaded from DANE. We will be here the sf library:

(stder = st_read("68_SANTANDER/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

## Reading layer `MGN_MPIO_POLITICO' from data source 
##   `/Users/ials/Documents/unal/GB2023/68_SANTANDER/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 87 features and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.52895 ymin: 5.707536 xmax: -72.47706 ymax: 8.14501
## Geodetic CRS:  WGS 84

Note the CRS of stder. As the ISRIC layers use the Homolosine projection, we need to reproject our layer using the sf library:

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

Let’s get the bounding box of our area of interest:

(bbox <- st_bbox(stder_igh))

##       xmin       ymin       xmax       ymax 
## -8313449.6   635360.0 -8089703.8   906698.4

Now, let’s use the bbox data to define the boundary box limits as used by the GDA libraryL. By the way, this is one of the trickiest parts of using 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 
## -8313449.6   906698.4 -8089703.8   635360.0

Now, we can use gdal_translate function to download the vrt layer. First, let’s define where to save the data:

sg_url="/vsicurl/https://files.isric.org/soilgrids/latest/data/"
datos = 'soc/soc_15-30cm_mean.vrt'
file = "soc_igh_15_30.tif"

Note that the following task may be take several minutes:

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

## Checking gdal_installation...

## Scanning for GDAL installations...

## Checking Sys.which...

## GDAL version 3.6.4

## GDAL command being used: "/usr/local/Cellar/gdal/3.6.4_4/bin/gdal_translate" -tr 250 250 -projwin -8313449.59600213 906698.380845578 -8089703.76012057 635360.004766508 -of "GTiff" -projwin_srs "+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs" "/vsicurl/https://files.isric.org/soilgrids/latest/data/soc/soc_15-30cm_mean.vrt" "soc_igh_15_30.tif"

## Input file size is 159243, 580340...10...20...30...40...50...60...70...80...90...100 - done.

## NULL

Let’s remind us what are the SOC layer units:

By dividing the predictions values by the values in the Conversion factor column, the user can obtain the more familiar units in the Conventional units column.

(stder_soc <- raster(file)/10)

## class      : RasterLayer 
## dimensions : 1085, 895, 971075  (nrow, ncol, ncell)
## resolution : 250, 250  (x, y)
## extent     : -8313500, -8089750, 635500, 906750  (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. Data exploration

A histogram can help at this point in time:

hist(stder_soc)

A summary is also helpful:

summary(stder_soc)

##         soc_igh_15_30
## Min.              8.7
## 1st Qu.          25.0
## Median           33.2
## 3rd Qu.          44.5
## Max.            130.3
## NA's          15388.0

Finally, plotting time. First a raster layer to stars conversion

(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   33.2 35.78297    44.5 130.3 15388
## dimension(s):
##   from   to   offset delta                       refsys x/y
## x    1  895 -8313500   250 +proj=igh +lon_0=0 +x_0=0... [x]
## y    1 1085   906750  -250 +proj=igh +lon_0=0 +x_0=0... [y]

names(nstder_soc) <-  "soc"

Then, the color palette:

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

Then, the map:

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")

Note that you can change the color palette and the map layout to better understand the soil layer’s information.

Now it’s your turn to try these functionalities. Good luck!!!

5. Reference

Cite this work as follows: Lizarazo, I. 2023. How to download soil data from SoilGrids. Available at: https://rpubs.com/ials2un/gettingSoilGrids

sessionInfo()

## R version 4.3.0 (2023-04-21)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.3.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## 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-3         mapview_2.11.0     RColorBrewer_1.1-3 dplyr_1.1.2       
##  [5] stars_0.6-1        sf_1.0-12          abind_1.4-5        raster_3.6-20     
##  [9] gdalUtils_2.0.3.2  rgdal_1.6-6        sp_1.6-0           XML_3.99-0.14     
## 
## loaded via a namespace (and not attached):
##  [1] xfun_0.39          bslib_0.4.2        htmlwidgets_1.6.2  lattice_0.21-8    
##  [5] vctrs_0.6.2        tools_4.3.0        crosstalk_1.2.0    generics_0.1.3    
##  [9] stats4_4.3.0       parallel_4.3.0     tibble_3.2.1       proxy_0.4-27      
## [13] fansi_1.0.4        highr_0.10         pkgconfig_2.0.3    R.oo_1.25.0       
## [17] KernSmooth_2.23-20 satellite_1.0.4    webshot_0.5.4      leaflet_2.1.2     
## [21] lifecycle_1.0.3    compiler_4.3.0     munsell_0.5.0      terra_1.7-29      
## [25] leafsync_0.1.0     codetools_0.2-19   htmltools_0.5.5    class_7.3-21      
## [29] sass_0.4.6         yaml_2.3.7         pillar_1.9.0       jquerylib_0.1.4   
## [33] R.utils_2.12.2     classInt_0.4-9     cachem_1.0.8       lwgeom_0.2-11     
## [37] iterators_1.0.14   foreach_1.5.2      tidyselect_1.2.0   digest_0.6.31     
## [41] fastmap_1.1.1      grid_4.3.0         colorspace_2.1-0   cli_3.6.1         
## [45] magrittr_2.0.3     base64enc_0.1-3    dichromat_2.0-0.1  utf8_1.2.3        
## [49] leafem_0.2.0       e1071_1.7-13       scales_1.2.1       rmarkdown_2.21    
## [53] png_0.1-8          R.methodsS3_1.8.2  evaluate_0.21      knitr_1.42        
## [57] tmaptools_3.1-1    viridisLite_0.4.2  rlang_1.1.1        Rcpp_1.0.10       
## [61] glue_1.6.2         DBI_1.1.3          rstudioapi_0.14    jsonlite_1.8.4    
## [65] R6_2.5.1           units_0.8-2