1. Introduction

This is an R Markdown Notebook which illustrates how to access SoilGrids using the webDAV protocol. In particular, it shows how to get soil layers for a region of interest. This notebook aims to help Geomatica Basica students at Universidad Nacional de Colombia to get started with Rgis.

We will use gdalUTils functions to access and download global soil data from the ISRIC website. ISRIC is the International Soil Reference and Information Centre. Its mission is to serve the international community as custodian of global soil information. ISRIC, led by the soil scientist Rik van den Bosch, are striving to increase awareness and understanding of soils in major global issues.

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. Predictions are made for the six standard depth intervals specified in the table below.

SoilGrids uses global models that are calibrated using all available input observations and globally available environmental covariates. This results in globally consistent predictions (no abrupt changes in predicted values at country boundaries, etc).

SoilGrids spatial predictions (layers) are produced using a reproducible soil mapping workflow, and can therefore be regularly updated as new soil data or covariates become available, after quality control and data standardisation/harmonisation.

SoilGrids can be accessed through the following services:

  • 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 here.

In case of any doubt, please review the FAQ SoilGrids

2. Preliminaries

Install packages (if this is your first time with this notebook):

#
#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)
#library(devtools)
#devtools::install_github("https://github.com/be-marc/leaflet.opacity", dependencies=TRUE)

Clean the workspace:

rm(list=ls())

Load packages:

#options("rgdal_show_exportToProj4_warnings"="none")
library(XML)
library(rgdal)
library(gdalUtils)
library(raster)
library(sf)
library(dplyr)
library(RColorBrewer)
library(leaflet.opacity)
library(leaflet)
library(leaflet.opacity)
library(mapview)

Set variables of interest:

We can use the table below for such a definition. It shows the properties currently mapped with SoilGrids, their description and mapped units. All maps produced with SoilGrids store data as integer values to minimise storage space. Therefore, some properties are provided in units that are not so common in soil science. 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.

### Get to know what are soil layers names

We need to define the url for the ISRIC webDAD site:

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

We can visit that url:

Let’s define nitrogen as our variable of interest. We also need to define the depth as well as the quantile.

Let’s note that SoilGrids maps have associated uncertainties as any product derived from a modelling approach. The prediction uncertainty is quantified by probability distributions. For each property and each standard depth interval this distribution is characterised by four parameters as follows:

  • Q0.05 - 5% quantile;
  • Q0.50 - median of the distribution;
  • mean - mean of the distribution;
  • Q0.95 - 95% quantile.

At the url, we can clic on the nitrogen option. We will get:

We need to define several variables:

voi = "nitrogen" # variable of interest
depth = "15-30cm"
quantile = "mean"

(variable = paste(url, voi, sep=""))
## [1] "https://files.isric.org/soilgrids/latest/data/nitrogen"
(layer = paste(variable,depth,quantile, sep="_")) # layer of interest 
## [1] "https://files.isric.org/soilgrids/latest/data/nitrogen_15-30cm_mean"
(vrt_layer = paste(layer, '.vrt', sep=""))
## [1] "https://files.isric.org/soilgrids/latest/data/nitrogen_15-30cm_mean.vrt"

3. Example 1: Describe the layer

We can use the raster library to get layer’s info:

nitro = raster("https://files.isric.org/soilgrids/latest/data/nitrogen/nitrogen_15-30cm_mean.vrt")
#nitro = raster(vrt_layer)

We can check layer’s info:

nitro

Note how large is this layer. It has 58034 rows by 159246 columns: more than 9240 million pixels. Note also the pixel size of the layer.

4. Example 2 :Download a Tiff for a region of interest (ROI)

Let’s read a shapefile of Colombian departments previously downloaded from DANE

deptos <- st_read("../datos/shapefiles/MGN2018_DPTO_POLITICO/MGN_DPTO_POLITICO.shp")
## Reading layer `MGN_DPTO_POLITICO' from data source `/Users/ivanlizarazo/Documents/ivan/GB/datos/shapefiles/MGN2018_DPTO_POLITICO/MGN_DPTO_POLITICO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 33 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -81.73562 ymin: -4.229406 xmax: -66.84722 ymax: 13.39473
## CRS:            4326
deptos
## Simple feature collection with 33 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -81.73562 ymin: -4.229406 xmax: -66.84722 ymax: 13.39473
## CRS:            4326
## First 10 features:
##    DPTO_CCDGO      DPTO_CNMBR DPTO_NANO_
## 1          18         CAQUETÁ       1981
## 2          19           CAUCA       1857
## 3          86        PUTUMAYO       1991
## 4          76 VALLE DEL CAUCA       1910
## 5          94         GUAINÍA       1991
## 6          99         VICHADA       1991
## 7          85        CASANARE       1991
## 8          91        AMAZONAS       1991
## 9          97          VAUPÉS       1991
## 10         95        GUAVIARE       1991
##                                                      DPTO_CACTO DPTO_NANO
## 1                            Ley 78 del 29 de Diciembre de 1981      2018
## 2                                           15 de junio de 1857      2018
## 3                    Articulo 309 Constitucion Politica de 1991      2018
## 4                         Decreto No 340 de 16 de Abril de 1910      2018
## 5                    Articulo 309 Constitucion Politica de 1991      2018
## 6                      5 de Julio Constitucion Politica de 1991      2018
## 7                      5 de Julio Constitucion Politica de 1991      2018
## 8  Dcto. 2274 del 4 de Octubre de la Constitución Política 1991      2018
## 9                    Articulo 309 Constitucion Politica de 1991      2018
## 10                     5 de Julio Constitucion Politica de 1991      2018
##    Shape_Leng Shape_Area                       geometry
## 1    21.38429   7.318485 MULTIPOLYGON (((-74.89423 2...
## 2    13.95026   2.534419 MULTIPOLYGON (((-76.45922 3...
## 3    12.70792   2.107965 MULTIPOLYGON (((-76.6705 1....
## 4    12.65087   1.679487 MULTIPOLYGON (((-77.2381 4....
## 5    21.17905   5.747937 MULTIPOLYGON (((-67.67638 3...
## 6    17.29261   8.100680 MULTIPOLYGON (((-67.80972 6...
## 7    12.13275   3.615063 MULTIPOLYGON (((-72.33885 6...
## 8    25.35598   8.877480 MULTIPOLYGON (((-71.14469 0...
## 9    20.12983   4.313810 MULTIPOLYGON (((-70.11033 2...
## 10   19.39679   4.511457 MULTIPOLYGON (((-71.31266 2...

Now, let’s filter any department (e.g. Boyaca):

(boyaca <- dplyr::filter(deptos, DPTO_CNMBR =="BOYACÁ"))
## Simple feature collection with 1 feature and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## CRS:            4326
##   DPTO_CCDGO DPTO_CNMBR DPTO_NANO_                    DPTO_CACTO DPTO_NANO
## 1         15     BOYACÁ       1886 Constitucion Politica de 1886      2018
##   Shape_Leng Shape_Area                       geometry
## 1   15.90649   1.888391 MULTIPOLYGON (((-72.17368 7...

We want to download soil raster files for Boyaca. As the ISRIC layers use the Homolosine projection, we need to reproject the Boyaca object:

igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
boyaca_igh <- st_transform(boyaca, igh)
boyaca_igh
## Simple feature collection with 1 feature and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -8325947 ymin: 518214 xmax: -8032661 ymax: 785421
## CRS:            +proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs
##   DPTO_CCDGO DPTO_CNMBR DPTO_NANO_                    DPTO_CACTO DPTO_NANO
## 1         15     BOYACÁ       1886 Constitucion Politica de 1886      2018
##   Shape_Leng Shape_Area                       geometry
## 1   15.90649   1.888391 MULTIPOLYGON (((-8057777 78...

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

(bbox <- st_bbox(boyaca_igh))
##     xmin     ymin     xmax     ymax 
## -8325947   518214 -8032661   785421

Now, let’s use the bbox data to define the boundary box limits as used by GDAL. 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 
## -8325947   785421 -8032661   518214

Now, we can use gdal_translate function to download the nitro layer.

sg_url="/vsicurl/https://files.isric.org/soilgrids/latest/data/"
datos = 'nitrogen/nitrogen_15-30cm_mean.vrt'
lfile = "./datos/nitro_ntest_igh_15_30.tif"

Note that the following task may be take several minutes:

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

If everything worked, we receive a message similar to this:

Let’s read the downloaded layer:

(boyaca_nitro <- raster(lfile)/100)
## class      : RasterLayer 
## dimensions : 1069, 1173, 1253937  (nrow, ncol, ncell)
## resolution : 250, 250  (x, y)
## extent     : -8326000, -8032750, 518250, 785500  (xmin, xmax, ymin, ymax)
## crs        : +proj=igh +ellps=WGS84 +units=m +no_defs 
## source     : memory
## names      : nitro_ntest_igh_15_30 
## values     : 0.91, 13.54  (min, max)

Now, let’s reproject the layer using WGS84 geographic coordinates. This is needed to plot it using mapview:

#crs_wgs84 = "+proj=longlat +datum=WGS84 +no_defs" 
crs_wgs84 = CRS('+init=EPSG:4326')
(boyaca_nitro_wgs84 = projectRaster(boyaca_nitro, crs=crs_wgs84))
## class      : RasterLayer 
## dimensions : 1077, 1233, 1327941  (nrow, ncol, ncell)
## resolution : 0.00226, 0.00225  (x, y)
## extent     : -74.72159, -71.93501, 4.644267, 7.067517  (xmin, xmax, ymin, ymax)
## crs        : +init=EPSG:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /private/var/folders/l_/1nv5l0dx78g_d3k3_t85lrjr0000gn/T/RtmpeGB1vm/raster/r_tmp_2020-10-19_165706_11534_46032.grd 
## names      : nitro_ntest_igh_15_30 
## values     : 0.8871728, 12.47475  (min, max)

A histogram can help at this point in time:

hist(boyaca_nitro_wgs84)

Finally, plotting time:

pal <- colorRampPalette(brewer.pal(9, "BrBG"))
mapview(boyaca_nitro_wgs84, col.regions = pal(100), at = seq(0, 6, 1), maxpixels =  1327941,legend = TRUE)

Note that you can change the background layers in the map to better understand the soil layer’s information.

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

5. Reproducibility

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] mapview_2.9.0         leaflet_2.0.3.9000    leaflet.opacity_0.1.0
##  [4] RColorBrewer_1.1-2    dplyr_1.0.2           sf_0.9-6             
##  [7] raster_3.3-13         gdalUtils_2.0.3.2     rgdal_1.5-18         
## [10] sp_1.4-4              XML_3.99-0.3         
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.1.0   xfun_0.18          purrr_0.3.4        lattice_0.20-41   
##  [5] colorspace_1.4-1   vctrs_0.3.4        generics_0.0.2     htmltools_0.5.0   
##  [9] stats4_3.6.3       yaml_2.2.1         base64enc_0.1-3    rlang_0.4.8       
## [13] R.oo_1.24.0        e1071_1.7-4        pillar_1.4.6       glue_1.4.2        
## [17] DBI_1.1.0          R.utils_2.10.1     foreach_1.5.1      lifecycle_0.2.0   
## [21] stringr_1.4.0      munsell_0.5.0      R.methodsS3_1.8.1  htmlwidgets_1.5.2 
## [25] codetools_0.2-16   evaluate_0.14      knitr_1.30         crosstalk_1.1.0.1 
## [29] class_7.3-17       leafem_0.1.1       Rcpp_1.0.5         KernSmooth_2.23-17
## [33] scales_1.1.1       classInt_0.4-3     satellite_1.0.2    webshot_0.5.2     
## [37] png_0.1-7          digest_0.6.26      stringi_1.5.3      grid_3.6.3        
## [41] tools_3.6.3        magrittr_1.5       tibble_3.0.4       crayon_1.3.4      
## [45] pkgconfig_2.0.3    ellipsis_0.3.1     rmarkdown_2.4.6    iterators_1.0.13  
## [49] R6_2.4.1           units_0.6-7        compiler_3.6.3