download_and_aggregate_soilgrids.R

#This script downloads soilgrids (https://www.isric.org/explore/soilgrids) data 
#at a given resolution (variable "res", in degrees) and the depth intervals 0-5 and 5-15
#to compute the mean for the depth interval 0-15.
#It also applies the conversion factor shown in the FAQ (https://www.isric.org/explore/soilgrids/faq-soilgrids#How_can_I_access_SoilGrids)
#to reconstitute the original data units.

library(gdalUtils)
library(raster)

#variables to download
vars <- c(
  "clay",
  "sand",
  "silt",
  "nitrogen",
  "phh2o",
  "soc"
)

#resolution
res <- 0.05
#res <- 0.008 #~1km takes around 24 hours

#iterating through variables
for(var.i in vars){
  
  #downloading two depths
  
  gdalwarp(
    t_srs = "EPSG:4326", 
    multi = TRUE,
    wm = 200, 
    co = c("BIGTIFF=YES", "COMPRESS=DEFLATE", "TILED=TRUE"),
    tr = c(res,res), # Desired output resolution
    verbose = T,
    paste0(
      "/vsicurl?max_retry=30&retry_delay=60&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/",
      var.i,
      "/",
      var.i,
      "_0-5cm_mean.vrt"
    ),
    paste0(
      var.i,
      "_0_5cm.tif"
    )
  )
  
  gdalwarp(
    t_srs = "EPSG:4326", 
    multi = TRUE,
    wm = 200, 
    co = c("BIGTIFF=YES", "COMPRESS=DEFLATE", "TILED=TRUE"),
    tr = c(res,res), # Desired output resolution
    verbose = T,
    paste0(
      "/vsicurl?max_retry=60&retry_delay=60&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/",
      var.i,
      "/",
      var.i,
      "_5-15cm_mean.vrt"
    ),
    paste0(
      var.i,
      "_5_15cm.tif"
    )
  )
  
  #mean
  
  temp.stack <- raster::stack(
    paste0(
      var.i,
      "_0_5cm.tif"
    ),
    paste0(
      var.i,
      "_5_15cm.tif"
    )
  )
  
  temp.raster <- raster::stackApply(
    x = temp.stack,
    indices = c(1, 1),
    fun = mean
  ) / ifelse(
    var.i == "nitrogen",
    100,
    10
  )
  
  raster::writeRaster(
    x = round(
      temp.raster,
      1
    ),
    filename = paste0(
      var.i,
      "_0_30cm.tif"
    )
  )
  
  #removing depth files
  
  file.remove(
    paste0(
      var.i,
      "_0_5cm.tif"
    ),
    paste0(
      var.i,
      "_5_15cm.tif"
    )
  )
  
}